{**************************************************}
{ }
{ MathPack 87 Lite Library for Virtual Pascal/2 }
{ }
{ Copyright (C) 1996 NITEK Corporation }
{ Permission is granted to use this library }
{ with Virtual Pascal/2. }
{ }
{**************************************************}
unit MP87LITE;
interface
{$PureInt+}
Type
NLINEQFUNC = procedure(Var x, f : Double; n: LongInt);
JACFUNC = procedure(Var x, f, h : Double; d2h, n : LongInt);
MALLOCFUNC = function( n : LongInt ): pointer;
FREEFUNC = procedure(p: pointer);
REPORTFUNC = procedure(msg: LongInt; x : Double);
{$OrgName+}
{---------------------------------------}
{ Vector/scalar manipulation procedures }
{---------------------------------------}
{ kfillv -- fill n element of x with constant k }
procedure kfillv (Var x: Double; k: Double; n: LongInt);
{ lfillv -- fill n elements with linear sequence }
procedure lfillv (Var x: Double; x0,dx : Double; n: LongInt);
{ zerov -- zero a vector }
procedure zerov(Var x: Double; n: LongInt);
{ copyv -- copy a vector }
procedure copyv(Var y, x: Double; n: LongInt);
{ sumv -- sum n elements of vector x }
function sumv (Var x: Double; n : LongInt) : Extended;
{ dotprod -- dot product of x and y }
function dotprod(Var x: Double; incx: LongInt;
Var y: Double; incy, n: LongInt) : Extended;
{ crossprod -- cross product of x, y and z }
procedure crossprod(Var x, y, z: Double);
{ vaddv -- x = y + z for vectors of length n }
procedure vaddv(Var x, y, z: Double; n: LongInt);
{ vsubv -- x = y - z for vectors of length n }
procedure vsubv(Var x, y, z: Double; n: LongInt);
{ vxv -- x = y * z for vectors of length n }
procedure vxv(Var x, y, z: Double; n: LongInt);
{ kxv -- y = k*x for vectors y and x }
procedure kxv(Var y: Double; k: Double; Var x: Double; n: LongInt);
{ magv -- Euclidean norm of vector x }
function magv(Var x: Double; n: LongInt) : Extended;
{ minval -- find minimum value in x }
function minval(Var x: Double; n: LongInt) : LongInt;
{ maxval -- find maximum value in x }
function maxval(Var x: Double; n: LongInt) : LongInt;
{ minaval -- find minimum absolute value in x }
function minaval(Var x: Double; n: LongInt) : LongInt;
{ maxaval -- find maximum absolute value in x }
function maxaval(Var x: Double; n: LongInt) : LongInt;
{ normalize -- normalize a vector }
procedure normalize(Var x: Double; n: LongInt);
{-------------------------}
{ Complex number routines }
{-------------------------}
Type
Complex = record
r,i : Single;
end;
LongComplex = record
r,i : Double;
end;
{ cmplx -- make LongComplex }
function cmplx(r,i: Double) : LongComplex;
{ cadd -- add complex numbers }
function cadd( a, b : LongComplex ) : LongComplex;
{ csub -- subtract complex numbers }
function csub( a, b : LongComplex ) : LongComplex;
{ cmul -- multiply complex numbers }
function cmul( a, b : LongComplex ) : LongComplex;
{ csqr -- square complex number }
function csqr( a : LongComplex ) : LongComplex;
{ cdiv -- divide complex number }
function cdiv( a, b : LongComplex ) : LongComplex;
{ conj -- conjugate complex number }
function conj( a : LongComplex ) : LongComplex;
{ csqrt -- square root of complex number }
function csqrt( a : LongComplex ) : LongComplex;
{----------------------------------}
{ Polynomial manipulation routines }
{----------------------------------}
{ poly87 -- evaluate polynomial of order n }
function poly87(Var a: Double; x: Double; n: LongInt): Extended;
{ polyder -- evaluate polynomial derivative of order n }
function polyder(Var a: Double; x: Double; n: LongInt): Extended;
{ laguerre -- find polynomial root by Laguerre's method }
procedure laguerre(Var a, x: LongComplex; m: LongInt);
{ findroots -- find all polynomia roots }
procedure findroots(Var a, ad, roots: LongComplex; m: LongInt);
{ poldiv -- divide polynomial by a polynomial }
procedure poldiv(Var q, u, v, r: LongComplex; nu, nv: LongInt);
{ polmul -- multiply polynomials }
procedure polmul(Var p, u, v: LongComplex; nu, nv: LongInt);
{-------------------------------------}
{ Simple matrix manipulation routines }
{-------------------------------------}
{ kfillm -- fill a matrix with a constant }
procedure kfillm(Var a: Double; k: Double; nr, nc, d2a: LongInt);
{ copym -- copy a matrix }
procedure copym(Var a, b: Double; nr, nc, d2a, d2b: LongInt);
{ kxm -- constant times matrix }
procedure kxm(Var a: Double; k: Double; Var b: Double;
nr, nc, d2a, d2b: LongInt);
{ mxv -- matrix times vector }
procedure mxv(Var y, a, x: Double; n, m, d2a: LongInt);
{ mxm -- matrix times matrix }
procedure mxm(Var a, b, c: Double; rb, cb, cc, d2a, d2b, d2c: LongInt);
{ maddm -- add matrix plus matrix }
procedure maddm(Var a, b, c: Double; nr, nc, d2a, d2b, d2c: LongInt);
{ msubm -- add matrix minus matrix }
procedure msubm(Var a, b, c: Double; nr, nc, d2a, d2b, d2c: LongInt);
{ transpm -- transpose matrix }
procedure transpm(Var a, b: Double; nrb, ncb, d2a, d2b: LongInt);
{ mxmt -- matrix x matrix transpose }
procedure mxmt(Var a, b, c: Double; nrb, ncb, nrc, d2a, d2b, d2c: LongInt);
{-------------------------------------}
{ Solving systems of linear equations }
{-------------------------------------}
{ decomp -- LU decomposition of a real matrix }
procedure decomp (Var a: Double; Var ip: LongInt;
n, d2a: LongInt; Var error: LongInt);
{ solve -- LU backsolving of a real matrix }
procedure solve (Var a, b: Double; Var ip: LongInt; n, d2a: LongInt);
{ invmatr -- Gauss-Jordan matrix inversion }
procedure invmatr (Var a: Double; n, d2a: LongInt; Var ip: LongInt;
Var det: Double);
{ tmdecomp -- LU decomposition of tridiagonal matrix }
procedure tmdecomp(Var a: Double; d2a, n: LongInt; Var error: LongInt);
{ tmsolve -- LU backsolving of a tridiagonal matrix }
procedure tmsolve (Var a, b: Double; d2a, n: LongInt);
{----------------------------------------}
{ Solving systems of nonlinear equations }
{----------------------------------------}
{ newton -- solution of system of nonlinear equations }
procedure newton(Var x, h: Double;
func: NLINEQFUNC;
jac: JACFUNC;
n, maxit, d2h: LongInt; Var iter: LongInt;
dxmax, tol: Double; goodh, verbose: LongInt;
malloc: MALLOCFUNC;
free: FREEFUNC;
report: REPORTFUNC
);
{ fdjacobian -- finite difference jacobian computation }
procedure fdjacobian(Var x, f, h: Double;
func: NLINEQFUNC;
d2h, n: LongInt;
malloc: MALLOCFUNC;
free: FREEFUNC
);
{---------------------------------------}
{ Statistical and data-fitting routines }
{---------------------------------------}
{ meansdev -- calculate mean and standard deviation }
procedure meansdev(Var x, mean, sdev: Double; n : LongInt);
{ linfit -- least-squares linear regression }
procedure linfit(Var x, y, m, b, r, sm, sb: Double; n : LongInt);
{--------------------------------}
{ Fast Fourier Transforms (FFTs) }
{--------------------------------}
{ fftc -- FFT of a complex vector }
procedure fftc(Var x: LongComplex; m, fwd: LongInt);
{ fft2d -- 2-dimensional FFT }
procedure fft2d(Var a, w: LongComplex; m1, m2, d2a, fwd: LongInt);
{-----------------------}
{ Numerical integration }
{-----------------------}
Type
ROMBFUNC = function ( x : Double ): Extended;
{ romberg -- Romberg numerical integration routine }
function romberg(a, b, tol: Double;
Var error: LongInt; func: ROMBFUNC) : Extended;
{----------------------------------------}
{ Ordinary Differential Equations (ODEs) }
{----------------------------------------}
Type
FYTFUNC = function ( y, t : Double ) : Extended;
{ hammde -- Hamming predictor-corrector ODE routine }
procedure hammde(n: LongInt; Var y: Double; y0, t0, dt: Double; fyt: FYTFUNC);
{ rk4step -- Runge Kutta fourth order method }
function rk4step(y0, t0, dt: Double; fyt: FYTFUNC): Extended;
implementation
begin
end.
