NumLib Documentation
UNDER CONSTRUCTION
Introduction
NumLib is a colletion of routines for numerical mathematics. It was written by C.J.J.M. van Ginneken, W.J.P.M. Kortsmit, and L. Reij at the Technical University of Einhoven (Netherlands) and donated to the Free Pascal project. Its beginnings date back to 1965. Originally written in Algol 60, ported to poorly formatted Pascal with strange procedure names, it is hard stuff, in particular because an official documentation is not available.
But after all, this library contains a series of very useful routines:
- Calculation of special functions
- Solution of systems of linear equations
- Finding roots of equations
- Fitting routines
- Matrix inversion
- Calculation of Eigen values
Therefore, it is the intention of this wiki article to create some useful documentation. Each routine described here was tested and its description will be accompanied by a working code example.
Data types and declarations (unit "typ")
Basic floating point type, ArbFloat
In NumLib, a new type ArbFloat
is defined which is used as floating point type throughout the entire library. This allows to switch between IEEE Double (64bit) and Extended(80bit) (though big endian state is unknown at present).
At the time of this writing the Define ArbExtended
is enabled, and therefore ArbFloat
is the type extended
.
In order to change the floating point type define or undefine ArbExtended
and add the machineconstants change to the type selected.
Basic integer type, ArbInt
Because in plain FPC mode the type Integer is 16-bits (for TP compatibility), and in Delphi 32-bits, all integers were changed to ArbInt
. The basic idea is the same as with ArbFloat
.
Vectors and matrices
NumLib often passes vectors as one ArbFloat
plus some integer value. This ArbFloat
value is the first value of an array containing the vector, the integer is the count of elements used from this array. Internally an array type is mapped over the first value in order to access the other array elements. This way both 0- and 1-based arrays as well as static and dynamic arrays can be passed to the functions and procedures of NumLib.
procedure DoSomething(var AVector: ArbFloat; n: Integer);
type
TLocalVector = array[0..MaxElements] of ArbFloat;
var
i: Integer;
x: TLocalVector absolute AVector;
begin
for i:=0 to n-1 do
// do something with x[i]
end;
var
x: array[1..100]; // Note: This array begins with index 1!
...
DoSomething(x[1], 100);
Similarly, a matrix is passed in the same way as the value of the first ArbFloat
value of the first row/column, plus two integers specifying the number of rows and columns. Data are considered to be arranged in rows and usually must be contiguous, i.e. the last element of a row must be followed immediately by the first element of the next row.
Complex numbers
NumLib defines the type complex
as a TP object
. Real and imaginary parts are denoted as xreal
and imag
, respectively. Methods for basic operations are provided:
type
complex = object
xreal, imag : ArbFloat;
procedure Init(r, i: ArbFloat); // Initializes the number with real and imaginary parts r and i, resp
procedure Add(c: complex); // Adds another complex number
procedure Sub(c: complex); // Subtracts a complex number
function Inp(z:complex):ArbFloat; // Calculates the inner product (a.Re * b.Re + a.Im * b.Im)
procedure Conjugate; // Conjugates the value (i.e. negates the imaginary part)
procedure Scale(s: ArbFloat); // Multiplies real and imaginary part by the real value s
Function Norm: ArbFloat; // Calculates the "norm", i.e. sqr(Re) + sqr(Im)
Function Size: ArbFloat; // Size = sqrt(Norm)
Function Re: ArbFloat; // Returns the real part of the complex number
procedure Unary; // Negates real and imaginary parts
Function Im: ArbFloat; // Returns the imaginary part of the complex number
Function Arg: ArbFloat; // Calculates the phase of the complex number
procedure MinC(c: complex); // Converts itself to the minimum of real and imaginary parts
procedure MaxC(c: complex); // Converts itself to the maximum of real and imaginary ports
Procedure TransF(var t: complex);
end;
The real and imaginary parts of a complex number can also be retrieved by these stand-alone functions:
function Re(c: complex): ArbFloat;
function Im(c: complex): ArbFloat;
Finding the roots of a function (unit "roo")
The roots are the x values at which a function f(x) is zero.
Roots of the quadratic equation
The quadratic equation
always has two, not necessarily different, complex roots. These solutions are determined by the procedure rooqua
.
procedure rooqua(p, q: ArbFloat; var z1, z2: complex);
p
andq
are the (real) coefficients of the quadratic equationz1
andz2
return the two complex roots. See unit typ for the declaration of the typecomplex
Example
Determine the roots of the equation z2 + 2 z + 5 = 0.
program quadratic_equ;
uses
SysUtils, typ, roo;
var
z1, z2: complex;
const
SIGN: array[boolean] of string = ('+', '-');
begin
rooqua(2, 5, z1, z2);
WriteLn(Format('1st solution: %g %s %g i', [z1.re, SIGN[z1.im < 0], abs(z1.im)]));
WriteLn(Format('2nd solution: %g %s %g i', [z2.re, SIGN[z2.im < 0], abs(z2.im)]));
end.
There are always two, not necessarily different, complex solutions.
Bisection
In the bisection method two x values a and b are estimated to be around the expected root such that the function values have opposite signs at a and b. The center point of the interval is determined, and the subinterval for which the function has opposite signs at its endpoints is selected for a new iteration. The process ends when a given precision, i.e. interval length, is achieved.
In NumLib, this approach is supported by the procedure roof1r
:
procedure roof1r(f: rfunc1r; a, b, ae, re: ArbFloat; var x: ArbFloat; var term: ArbInt);
f
is the function for which the root is to be determined. It must be a function of one floating point argument (typeArbFloat
). The type of the function,rfunc1r
, is declared in unittyp
.a
andb
are the endpoints of the test interval. The root must be located between these two values, i. e. the function values f(a) and f(b) must have different signs.ae
andre
determine the absolute and relative precision, respectively, with which the root will be determined.re
is relative to the maximum of abs(a) and abs(b). Note that precision and speed are conflicting issues. Highest accuracy is achieved ifae
is given asMachEps
(see unittyp
). Both parameters must not be negative.x
returns the value of the found root.term
returns whether the process has been successful:- 1 - successful termination, a zero point has been found with absolute accuracy
ae
or relative accuracyre
- 2 - the required accuracy of the root could not be reached; However, the value of x is called the "best achievable" approach
- 3 - error in input parameters:
ae
< 0 orre
< 0, or f(a)*f(b) > 0
- 1 - successful termination, a zero point has been found with absolute accuracy
Example
The following program determines the square root of 2. This is the x value at which the function f(x) = x^2 - 2 is zero. Since f(1) = 1^2 - 2 = -1 < 0 and f(2) = 2^2 - 2 = 2 > 0 we can assume a and b to be 1 and 2, respectively.
program bisection_demo;
uses
typ, roo;
function f(x: ArbFloat): ArbFloat;
begin
Result := x*x - 2;
end;
var
x: ArbFloat = 0.0;
term : ArbInt;
begin
roof1r(@f, 1.0, 2.0, 1e-9, 0, x, term);
WriteLn('Bisection result ', x);
WriteLn('sqrt(2) ', sqrt(2.0));
end.
Special functions (unit "spe")
Error function
The Gauss error function erf(x) is the integral over the Gauss function, normalized to unity:
Sometimes also its complement, the complemenatry error function erfc(x), is needed:
Both functions can be calculated by the NumLib functions speerf
and speefc
, respectively
function speerf(x: ArbFloat): ArbFloat; // --> erf(x)
function speefc(x: ArbFloat): ArbFloat; // --> erfc(x)
Example
A table of some values of the error function and the complementary error function is calculated in the following demo project:
program erf_Table;
uses
SysUtils, StrUtils, typ, spe;
const
Wx = 7;
Wy = 20;
D = 6;
var
i: Integer;
x: ArbFloat;
fs: TFormatSettings;
begin
fs := DefaultFormatSettings;
fs.DecimalSeparator := '.';
WriteLn('x':Wx, 'erf(x)':Wy+1, 'erfc(x)':Wy+1);
WriteLn;
for i:= 0 to 20 do begin
x := 0.1*i;
WriteLn(Format('%*.1f %*.*f %*.*f', [
Wx, x,
Wy, D, speerf(x),
Wy, D, speefc(x)
], fs));
end;
end.
Gamma function
The gamma function is needed by many probability functions. It is defined by the integral
NumLib provides two functions for its calculation:
function spegam(x: ArbFloat): ArbFloat;
function spelga(x: ArbFloat): ArbFloat;
The first one, spegam
, calculates the function directly. But since the gamma function grows rapidly for even not-too large arguments this calculation very easily overflows.
The second function, spelga
calculates the natural logarithm of the gamma function which is more suitable to combinatorial calculations where multiplying and dividing the large gamma values can be avoided by adding or subtracting their logarithms.
Example
The following demo project prints a table of some values of the gamma function and of its logarithm. Note that spegam(x)
overflows above about x = 170 for data type extended
.
program Gamma_Table;
uses
SysUtils, StrUtils, typ, spe;
const
VALUES: array[0..2] of ArbFloat = (1, 2, 5);
Wx = 7;
Wy = 30;
Wln = 20;
var
i: Integer;
x: ArbFloat;
magnitude: ArbFloat;
begin
WriteLn('x':Wx, 'Gamma(x)':Wy, 'ln(Gamma(x))':Wln);
WriteLn;
magnitude := 1E-3;
while magnitude <= 1000 do begin
for i := 0 to High(VALUES) do begin
x := VALUES[i] * magnitude;
Write(PadLeft(FloatToStr(x), Wx));
if x <= 170 then // Extended overflow above 170
Write(FormatFloat('0.000', spegam(x)):Wy)
else
Write('overflow':Wy);
WriteLn(spelga(x):Wln:3);
if abs(x-1000) < 1E-6 then
break;
end;
magnitude := magnitude * 10;
end;
end.
Incomplete gamma function
The lower and upper incomplete gamma functions are defined by
In fpc 3.2, these functions were added to NumLib for their calculation
function gammap(s, x: ArbFloat): ArbFloat;
function gammaq(s, x: ArbFloat): ArbFloat;
Example
The following project prints out a table of gammap
and gammaq
values:
program IncompleteGamma_Table;
uses
SysUtils, StrUtils, typ, spe;
const
s = 1.25;
Wx = 7;
Wy = 20;
D = 6;
var
i: Integer;
x: ArbFloat;
fs: TFormatSettings;
begin
fs := DefaultFormatSettings;
fs.DecimalSeparator := '.';
WriteLn('x':Wx, 'GammaP('+FloatToStr(s,fs)+', x)':Wy+1, 'GammaQ('+FloatToStr(s,fs)+', x)':Wy+1);
WriteLn;
for i := 0 to 20 do begin
x := 0.5*i;
WriteLn(Format('%*.1f %*.*f %*.*f', [
Wx, x,
Wy, D, gammap(s, x),
Wy, D, gammaq(s, x)
], fs));
end;
end.