NumLib Documentation: what changed

Revision 2 to 3. Removed lines have a minus sign; added lines have a plus sign.

Old	New		Line
1	1		## Introduction
2	2
3	3		NumLib is a collection 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 Eindhoven (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 readily available (well... -- Dutch documentation files in TeX can be found at [<http://www.stack.nl/~marcov/numlib/>](http://www.stack.nl/~marcov/numlib/)).
4	4
5	5		But after all, this library contains a series of very useful routines grouped into several units:
6	6
7	7		- Calculation of special functions (unit spe)
8	8		- Solution of systems of linear equations (unit sle)
9	9		- Finding roots of equations (unit roo)
10	10		- Fitting routines (unit ipf)
11	11		- Basic operations with matrices and vectors (unit omv)
12	12		- Matrix inversion (unit inv)
13	13		- Calculation of determinants (unit det)
14	14		- Calculation of eigenvalues (unit eig)
15	15
16	16		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.
17	17
18	18		> Note: This page heavily makes use of mathematical formulas written in LaTeX. In order to see these formulas in the browser JavaScript must be enabled.
19	19
20	20		### Naming convention of routines
21	21
22	22		Would you imagine that a procedure named `spears` calculates the arcsin? The original NumLib routines limit the procedure/function names to 6 characters (maybe a limitation of Algol?). The first three characters always are identical with the name of the unit in which they are implemented. This leaves only three characters to differentiate between the routines...
23	23
24	24		Routines later added to the original library are given more conventional names, such as `normaldist` for calculation of the normal distribution.
25	25
26	26		## Data types and declarations (unit typ)
27	27
28	28		### Basic floating point type, ArbFloat
29	29
30	30		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).
31	31
32	32		At the time of this writing the `$DEFINE ArbExtended` is enabled, and therefore `ArbFloat` is identical to the standard type `extended`. In order to change the floating point type define or undefine `ArbExtended` and adapt the machineconstants to the type selected.
33	33
34	34		### Basic integer type, ArbInt
35	35
36	36		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`.
37	37
38	38		### Vectors
39	39
40	40		NumLib often passes vectors to procedures 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.
41	41
42	42		``` pascal
43	43		procedure DoSomething(var AVector: ArbFloat; n: Integer);
44	44		type
45	45		TLocalVector = array[0..MaxElements] of ArbFloat;
46	46		var
47	47		i: Integer;
48	48		x: TLocalVector absolute AVector;
49	49		begin
50	50		for i:=0 to n-1 do
51	51		// do something with x[i]
52	52		end;
53	53
54	54		var
55	55		x: array[1..100]; // Note: This array begins with index 1!
56	56		...
57	57		DoSomething(x[1], 100);
58	58		```
59	59
60	60		### Matrices
61	61
62	62		A matrix is a two-dimensional arrangement of numbers in m rows and n columns.
63	63
64	64
65	65		$`{A} = \left[
66	66		\begin{array}{cccc}
67	67		a_{11} & a_{12} & a_{13} & \ldots & a_{1n}\\
68	68		a_{21} & a_{22} & a_{23} & \ldots & a_{2n}\\
69	69		\vdots & \vdots & \vdots & \ddots & \vdots\\
70	70		a_{m1} & a_{m2} & a_{m3} & \ldots & a_{mn}
71	71		\end{array}
72	72		\right]`$
73	73
74	74		The matrix elements a<sub>ij</sub> are addressed by row and column indexes, i and j, respectively. In mathematics, the index values usually start with 1 (i = 1..m, j = 1..n) while in computing the first index is preferred to be 0 (i = 0...m-1, j = 0...n-1).
75	75
76	76		#### Generic matrices
77	77
78	78		In this type of matrix, no assumption is made on symmetry and shape of the data arrangement.
79	79
80	80		A constant (real) matrix can be declared in Pascal as a two-dimensional array of floating point numbers: `A: array[1..m, 1..n] of ArbFloat` (or `A: array[0..m-1, 0..n-1] of ArbFloat` in 0-based index notation). Similarly, a constant complex matrix can be declared as an `array[1..m, 1..n] of complex`.
81	81
82	82		When a NumLib procedure needs a matrix the first (top left) matrix element must be passed as a var parameter, along with the number of rows and columns. It is even possible to allocate larger arrays as needed, and in this case the allocated row length (or number of columns) must be specified.
83	83
84	84		Alternatively a matrix also can be specified as a one-dimensional array. This works because all data are found within the same memory block. The array index, `k` must be converted to the 2D matrix indexes, `i, j`:
85	85
86	86		``` pascal
87	87		k := i * n + j <==> i = k div n, j := k mod n {for 0-based array and matrix indexes}
88	88		k := (i-1) * n + j <==> i := (k-1) div n + 1, j := (k-1) mod n + 1 {for 1-based array and matrix indexes}
89	89		```
90	90
91	91		Using a one-dimensional array is also the correct way to dynamically allocate memory for a matrix. The first idea of using an `array of array of ArbFloat` will not work because this does not allocate a single memory block. Instead, individual memory blocks will be allocated for the inner array (rows), and these pointers will be stored in the outer array - this way of storing 2D data is not compatible with the requirements of NumLib.
92	92
93	93		Here is a simple example how to create, populate and read a matrix dynamically
94	94
95	95		``` math
96	96		A = \left[
97	97		\begin{array}{ccc}
98	98		1 & 2 & 3 \\
99	99		4 & 5 & 6
100	100		\end{array}
101	101		\right]
102	102		```
103	103
104	104		``` pascal
105	105		var
106	106		A: array of Arbfloat; // Note: the index in a dynamic array is 0-based.
107	107		ij: Integer;
108	108		i, j: Integer;
109	109		p: ^ArbFloat;
110	110		m, n: Integer;
111	111		begin
112	112		m := 2;
113	113		n := 3;
114	114		SetLength(A, m, n);
115	115		A[0] := 1; A[1] := 2; A[2] := 3;
116	116		A[3] := 4; A[4] := 5; A[5] := 6;
117	117
118	118		for i := 0 to m-1 do begin
119	119		for j := 0 to n-1 do begin
120	120		ij := i * n + j;
121	121		Write(A[ij]);
122	122		end;
123	123		WriteLn;
124	124		end;
125	125
126	126		// or
127	127		p := @A[0]; // pointer to first array element
128	128		for i:=0 to m - 1 do begin
129	129		for j := 0 to n - 1 do begin
130	130		WriteLn(p^);
131	131		inc(p); // advance pointer to next array element
132	132		end;
133	133		WriteLn;
134	134		end;
135	135		```
136	136
137	137		#### Symmetric matrices
138	138
139	139		The elements of a symmetric matrix are mirrored at the main diagonal, i.e. a<sub>ij</sub> = a<sub>ji</sub>. In other words: a symmetric matrix is equal to its transpose: A = A<sup>T</sup>.
140	140
141	141		NumLib contains some routine optimized for such matrices. They require the matrix being stored like a generic matrix in a rectangular m x n array or a linear m*n array. However, only the half of the matrix below the diagonal must be populated (but it can be, of course).
142	142
143	143		#### Symmetric positive definite matrices
144	144
145	145		A matrix A is positive definite if the product d<sup>T</sup> A d is a positive number for any non-zero vector d. A matrix is symmetric positive definite if, additionally, A<sup>T</sup> = A.
146	146
147	147		This type of matrices has some interesting properties (determinant is never 0, eigenvalues are real and positive). Therefore NumLib offers specialized routines for many matrix problems with symmetric positive matrices which are superior to those for general matrices.
148	148
149	149		The matrix is stored in the conventional way like a rectangular 2D array of m rows and n columns or a linear 1D array with m*n elements.
150	150
151	151		#### Band matrices
152	152
153	153		A n x n band matrix has zero elements everywhere except for a band along the main diagonal characterized by a left and right bandwidth (l and r, respectively).
154	154
155	155		The NumLib procedures for band matrices assume that the matrix stores only the elements of the band. It is expected the the matrix is passed as a 1D array. Starting at the top/left corner point of the matrix (element A<sub>11</sub>) the diagonal elements are collected by running along the rows of the matrix, from left to right and from top to bottom. In total, there will be n (l+1+r) - (r(r+1) + l(l+1)) div 2 elements. The calling procedure must allocate a 1D array of at least this size for the diagonal elements. Note that a band matrix cannot be passed as a 2D array.
156	156
157	157		As an example, suppose this 7 x 7 matrix with left bandwidth 2 and right bandwidth 1:
158	158
159	159		$`A = \begin{bmatrix}
160	160		A_{11} & A_{12} & 0 & 0 & 0 & 0 & 0 \\
161	161		A_{21} & A_{22} & A_{23} & 0 & 0 & 0 & 0 \\
162	162		A_{31} & A_{32} & A_{33} & A_{34} & 0 & 0 & 0 \\
163	163		0 & A_{42} & A_{43} & A_{44} & A_{45} & 0 & 0 \\
164	164		0 & 0 & A_{53} & A_{54} & A_{55} & A_{56} & 0 \\
165	165		0 & 0 & 0 & A_{64} & A_{65} & A_{66} & A_{67} \\
166	166		0 & 0 & 0 & 0 & A_{75} & A_{76} & A_{77}
167	167		\end{bmatrix}`$
168	168
169	169		The 1D data array to be passed to the NumLib routines consists of these elements
170	170
171	171		``` pascal
172	172		const
173	173		arr: array[1..24] of ArbFloat = (
174	174		A11, A22,
175	175		A21, A22, A23,
176	176		A31, A32, A33, A34,
177	177		A42, A43, A44, A45,
178	178		A53, A54, A55, A56,
179	179		A64, A65, A66, A67,
180	180		A75, A76, A77
181	181		);
182	182		```
183	183
184	184		#### Symmetric band matrices
185	185
186	186		In a symmetric matrix, the matrix elements can be mirrored at the main diagonal of the matrix, i.e. A<sub>ij</sub> = A<sub>ji</sub>. All elements out side a band of width w to the left and right of the main diagonal are zero.
187	187
188	188		Like with a generic band matrix, NumLib assumes that the only the elements of the diagonals are passed as a 1D array to the procedures dedicated to this matrix type. Since in a symmetric band matrix the left and right band widths are equal it is assumed that only the left band and the main diagonal are stored. Again the relevant diagonal elements are collected by running across the rows of the matrix from left to right and from top to bottom, stopping in each row at the main diagonal. In total, there will be n (w+1) - w(w+1) div 2 elements in the data array.
189	189
190	190		As an example consider this 7 x 7 symmetric band matrix with band width 2:
191	191
192	192		$`A = \begin{bmatrix}
193	193		A_{11} & A_{21} & A_{31} & 0 & 0 & 0 & 0 \\
194	194		A_{21} & A_{22} & A_{32} & A_{42} & 0 & 0 & 0 \\
195	195		A_{31} & A_{32} & A_{33} & A_{43} & A_{53} & 0 & 0 \\
196	196		0 & A_{42} & A_{43} & A_{44} & A_{54} & A_{64} & 0 \\
197	197		0 & 0 & A_{53} & A_{54} & A_{55} & A_{65} & A_{75} \\
198	198		0 & 0 & 0 & A_{64} & A_{65} & A_{66} & A_{76} \\
199	199		0 & 0 & 0 & 0 & A_{75} & A_{76} & A_{77}
200	200		\end{bmatrix}`$
201	201
202	202		The 1D array for the NumLib procedures is
203	203
204	204		``` pascal
205	205		const
206	206		arr: array[1..18] of ArbFloat = (
207	207		A11,
208	208		A21, A22,
209	209		A31, A32, A33,
210	210		A42, A43, A44,
211	211		A53, A54, A55,
212	212		A64, A65, A66,
213	213		A75, A76, A77
214	214		);
215	215		```
216	216
217	217		#### Tridiagonal matrices
218	218
219	219		A tridiagonal matrix has zero elements everywhere except for the main diagonal and the diagonal immediately above and below it. This way the matrix can be stored more efficiently than a generic matrix because all the empty elements can be dropped. The elements of an n x n square tridiagonal matrix can be stored as three vectors
220	220
221	221		- an array `l` with `n-1` elements for the subdiagonal
222	222		- an array `d` with `n` elements for the main diagonal, and
223	223		- an array `u` with `n-1` elements for the superdiagonal
224	224
225	225
226	226		$`A=
227	227		\begin{bmatrix}
228	228		d[0] & u[0] & & \dots & 0 \\
229	229		l[1] & d[1] & u[1] & & \\
230	230		& l[2] & d[2] & u[2] & \\
231	231		\vdots & & \ddots & \ddots & \ddots \\
232	232		& & l[n-1] & d[n-1] & u[n-1] \\
233	233		0 & & & l[n] & d[n]
234	234		\end{bmatrix}`$
235	235
236	236		#### Symmetric tridiagonal matrices
237	237
238	238		In these tridiagonal matrices, the subdiagonal and superdiagonal elements are equal. Consequently, NumLib requests only two, instead of three, one-dimensional arrays:
239	239
240	240		- an array `d` with `n` elements for the main diagonal, and
241	241		- an array `c` with `n-1` elements for the superdiagonal and subdiagonal
242	242
243	243
244	244		$`A=
245	245		\begin{bmatrix}
246	246		d[0] & c[1] & & \dots & 0 \\
247	247		c[1] & d[1] & c[2] & & \\
248	248		& c[2] & d[2] & c[3] & \\
249	249		\vdots & & \ddots & \ddots & \ddots \\
250	250		& & c[n-1] & d[n-1] & c[n] \\
251	251		0 & & & c[n] & d[n]
252	252		\end{bmatrix}`$
253	253
254	254		### Complex numbers
255	255
256	256		NumLib defines the type `complex` as a Turbo Pascal `object`. Real and imaginary parts are denoted as `xreal` and `imag`, respectively. Methods for basic operations are provided:
257	257
258	258		``` pascal
259	259		type
260	260		complex = object
261	261		xreal, imag : ArbFloat;
262	262		procedure Init(r, i: ArbFloat); // Initializes the number with real and imaginary parts r and i, resp
263	263		procedure Add(c: complex); // Adds another complex number
264	264		procedure Sub(c: complex); // Subtracts a complex number
265	265		function Inp(z:complex):ArbFloat; // Calculates the inner product (a.Re * b.Re + a.Im * b.Im)
266	266		procedure Conjugate; // Conjugates the value (i.e. negates the imaginary part)
267	267		procedure Scale(s: ArbFloat); // Multiplies real and imaginary part by the real value s
268	268		function Norm: ArbFloat; // Calculates the "norm", i.e. sqr(Re) + sqr(Im)
269	269		function Size: ArbFloat; // Size = sqrt(Norm)
270	270		function Re: ArbFloat; // Returns the real part of the complex number
271	271		procedure Unary; // Negates real and imaginary parts
272	272		function Im: ArbFloat; // Returns the imaginary part of the complex number
273	273		function Arg: ArbFloat; // Calculates the phase phi of the complex number z = x + y i = r exp(i phi)
274	274		procedure MinC(c: complex); // Converts itself to the minimum of real and imaginary parts of itself and another value, c
275	275		procedure MaxC(c: complex); // Converts itself to the maximum of real and imaginary ports of itself and another value, c
276	276		procedure TransF(var t: complex);
277	277		end;
278	278		```
279	279
280	280		The real and imaginary parts of a complex number can also be retrieved by these stand-alone functions:
281	281
282	282		``` pascal
283	283		function Re(c: complex): ArbFloat;
284	284		function Im(c: complex): ArbFloat;
285	285		```
286	286
287	287		## Basic operations with matrices and vectors (unit omv)
288	288
289	289		### Inner product of two vectors
290	290
291	291		The result of the inner product of two vectors a = \[a<sub>1</sub>, a<sub>2</sub>, ... a<sub>n</sub>\] and b = \[b<sub>1</sub>, b<sub>2</sub>, ... b<sub>n</sub>\] (also called "dot product" or "scalar product") is a scalar, i.e. a single number. It is calculated as the sum of the element-by-element products of the vectors. Both vectors must have the same lengths.
292	292
293	293
294	294		$`\mathbf{a} \cdot \mathbf{b} = \sum_{i=1}^n a_i b_i = a_1 b_1 + a_2 b_2 + \dots + a_n b_n`$
295	295
296	296		NumLib provides the function `omvinp` for calculation of the inner product:
297	297
298	298		``` pascal
299	299		function omvinp(var a, b: ArbFloat; n: ArbInt): ArbFloat;
300	300		```
301	301
302	302		- `a` and `b` are the first elements of 1-dimensional arrays representing the vectors a and b, respectively.
303	303		- `n` defines the dimension of the vectors (count of array elements). Both vectors must have the same number of elements.
304	304
305	305		Example
306	306
307	307		Calculate the dot product of the vectors a = \[0, 1, 2, 2, -1\] and b = \[3, -1, -2, 2, -1\].
308	308
309	309		``` pascal
310	310		program dot_product;
311	311
312	312		uses
313	313		typ, omv;
314	314
315	315		const
316	316		n = 5;
317	317		a: array[0..n-1] of ArbFloat = (0, 1, 2, 2, -1);
318	318		b: array[0..n-1] of Arbfloat = (3, -1, -2, 2, -1);
319	319
320	320		var
321	321		ab: ArbFloat;
322	322		i: Integer;
323	323
324	324		begin
325	325		ab := omvinp(a[0], b[0], n);
326	326
327	327		Write('Vector a = [');
328	328		for i:= Low(a) to High(a) do
329	329		Write(a[i]:5:0);
330	330		WriteLn(' ]');
331	331
332	332		Write('Vector b = [');
333	333		for i:= Low(b) to High(b) do
334	334		Write(b[i]:5:0);
335	335		WriteLn(' ]');
336	336
337	337		Write('Dot product a b = ');
338	338		WriteLn(ab:0:0);
339	339		end.
340	340		```
341	341
342	342		### Product of two matrices
343	343
344	344		The product of a m x n matrix A and a n x p matrix B is a m x p matrix C in which the elements are calculated as the inner product of the row vectors of A and the column vectors of B (see [here](https://en.wikipedia.org/wiki/Matrix_multiplication) for more details):
345	345
346	346
347	347		$`C_{ij} = \sum_{k=0}^n A_{ik} B_{kj}`$
348	348
349	349		The NumLib function to perform this multipication is `omvmmm` ("mmm" = multiply matrix by matrix):
350	350
351	351		``` pascal
352	352		procedure omvmmm(var a: ArbFloat; m, n, rwa: ArbInt;
353	353		var b: ArbFloat; p, rwb: ArbInt;
354	354		var c: ArbFloat; rwc: ArbInt);
355	355		```
356	356
357	357		- `a` is the first element of the first input matrix A. The array will not be changed.
358	358		- `m` specifies the number of rows (or: column length) of matrix A.
359	359		- `n` specifies the number or columns (or: row length) of matrix A.
360	360		- `rwa` is the allocated row length of matrix A. This means that the array can be allocated larger than actually needed. Of course, `rwa` cannot be smaller than `n`.
361	361
362	362		- `b` is the first element of the second input matrix B. Again, this array will not be changed.
363	363		- `p` specifies the number or columns (or: row length) of matrix B. The number of rows is already defined by the number of columns of matrix A, `n`.
364	364		- `rwb` is the allocated row length of matrix B. `rwb` cannot be smaller than `p`.
365	365
366	366		- `c` is the first element of the output matrix C. After the calculation, this and the following array elements will contain the elements of the product matrix.
367	367		- The size of the output matrix is `m` rows x `p` columns as set up by the input matrices A and B.
368	368		- `rwc` is the allocated row length of matrix C. `rwc` cannot be smaller than `p`.
369	369
370	370		Example
371	371
372	372		Calculate the product of the matrices
373	373
374	374
375	375		$`A = \begin{bmatrix}
376	376		3 & \ 5 & 4 & 1 & -4 \\
377	377		-2 & \ 3 & 4 & 1 & 0 \\
378	378		0 & \ 1 & -1 & -2 & 5
379	379		\end{bmatrix}
380	380		\ \ \ \ \
381	381		B = \begin{bmatrix}
382	382		3 & 4 & -2 & 0 \\
383	383		0 & -2 & 2 & 1 \\
384	384		-1 & 0 & 3 & 4 \\
385	385		-2 & -3 & 0 & 1 \\
386	386		1 & 3 & 0 & -1
387	387		\end{bmatrix}`$
388	388
389	389		``` pascal
390	390		program matrix_multiplication;
391	391
392	392		uses
393	393		typ, omv;
394	394
395	395		const
396	396		m = 3;
397	397		n = 5;
398	398		p = 4;
399	399
400	400		A: array[1..m, 1..n] of ArbFloat = (
401	401		( 3, 5, 4, 1, -4),
402	402		(-2, 3, 4, 1, 0),
403	403		( 0, 1, -1, -2, 5)
404	404		);
405	405		B: array[1..n, 1..p] of ArbFloat = (
406	406		( 3, 4, -2, 0),
407	407		( 0, -2, 2, 1),
408	408		(-1, 0, 3, 4),
409	409		(-2, -3, 0, 1),
410	410		( 1, 3, 0, -1)
411	411		);
412	412
413	413		var
414	414		C: array[1..m, 1..p] of ArbFloat;
415	415		i, j: Integer;
416	416
417	417		begin
418	418		// Calculate product
419	419		omvmmm(
420	420		A[1,1], m, n, n,
421	421		B[1,1], p, p,
422	422		C[1,1], p
423	423		);
424	424
425	425		// Print A
426	426		WriteLn('A = ');
427	427		for i:= 1 to m do begin
428	428		for j := 1 to n do
429	429		Write(A[i, j]:10:0);
430	430		WriteLn;
431	431		end;
432	432		WriteLn;
433	433
434	434		// Print B
435	435		WriteLn('B = ');
436	436		for i:= 1 to n do begin
437	437		for j := 1 to p do
438	438		Write(B[i, j]:10:0);
439	439		WriteLn;
440	440		end;
441	441		WriteLn;
442	442
443	443		// Print C
444	444		WriteLn('C = A B = ');
445	445		for i:= 1 to m do begin
446	446		for j := 1 to p do
447	447		Write(C[i, j]:10:0);
448	448		WriteLn;
449	449		end;
450	450		end.
451	451		```
452	452
453	453		### Product of a matrix with a vector
454	454
455	455		A matrix A of size m x n multiplied with a vector b of length n yields another vector of length m. The elements of the result vector c are the dot products of the matrix rows with b. Note that the length of vector b must be equal to the number of columns of matrix A.
456	456
457	457
458	458		$`\begin{align*}
459	459		\mathbf{c} = A\ \mathbf{b}=
460	460		\left[
461	461		\begin{array}{cccc}
462	462		a_{11} & a_{12} & a_{13} & \ldots & a_{1n}\\
463	463		a_{21} & a_{22} & a_{23} & \ldots & a_{2n}\\
464	464		\vdots & \vdots & \vdots & \ddots & \vdots\\
465	465		a_{m1} & a_{m2} & a_{m3} & \ldots & a_{mn}
466	466		\end{array}
467	467		\right]
468	468		\ \
469	469		\left[
470	470		\begin{array}{c}
471	471		b_1\\
472	472		b_2\\
473	473		b_3\\
474	474		\vdots\\
475	475		b_n
476	476		\end{array}
477	477		\right]
478	478		=
479	479		\left[
480	480		\begin{array}{c}
481	481		a_{11}b_1+a_{12}b_2 + a_{13}b_3 + \cdots + a_{1n} b_n\\
482	482		a_{21}b_1+a_{22}b_2 + a_{23}b_3 + \cdots + a_{2n} b_n\\
483	483		\vdots\\
484	484		a_{m1}b_1+a_{m2}b_2 + a_{m3}b_3 + \cdots + a_{mn} b_n\\
485	485		\end{array}
486	486		\right].
487	487		\end{align*}`$
488	488
489	489		NumLib's procedure for this task is called `omvmmv` ("mmv" = multiply a matrix with a vector):
490	490
491	491		``` pascal
492	492		procedure omvmmv(var a: ArbFloat; m, n, rwidth: ArbInt; var b, c: ArbFloat);
493	493		```
494	494
495	495		Example
496	496
497	497		Calculate the product
498	498
499	499
500	500		$`\begin{align*}
501	501		\mathbf{c} = A\ \mathbf{b}=
502	502		\left[
503	503		\begin{array}{cccc}
504	504		3 & 5 & 4 & 1 & -4\\
505	505		-2 & 3 & 4 & 1 & 0\\
506	506		0 & 1 & -1 & -2 & 5
507	507		\end{array}
508	508		\right]
509	509		\ \
510	510		\left[
511	511		\begin{array}{c}
512	512		3\\
513	513		0\\
514	514		-1\\
515	515		-2\\
516	516		1
517	517		\end{array}
518	518		\right]
519	519		\end{align*}`$
520	520
521	521		``` pascal
522	522		program matrix_vector_multiplication;
523	523
524	524		uses
525	525		typ, omv;
526	526
527	527		const
528	528		m = 3;
529	529		n = 5;
530	530		p = 4;
531	531
532	532		A: array[1..m, 1..n] of ArbFloat = (
533	533		( 3, 5, 4, 1, -4),
534	534		(-2, 3, 4, 1, 0),
535	535		( 0, 1, -1, -2, 5)
536	536		);
537	537		b: array[1..n] of ArbFloat = (
538	538		3, 0, -1, -2, 1
539	539		);
540	540
541	541		var
542	542		c: array[1..m] of ArbFloat;
543	543		i, j: Integer;
544	544
545	545		begin
546	546		// Calculate product c = A b
547	547		omvmmv(
548	548		A[1,1], m, n, n,
549	549		b[1],
550	550		c[1]
551	551		);
552	552
553	553		// Print A
554	554		WriteLn('A = ');
555	555		for i:= 1 to m do begin
556	556		for j := 1 to n do
557	557		Write(A[i, j]:10:0);
558	558		WriteLn;
559	559		end;
560	560		WriteLn;
561	561
562	562		// Print vector b
563	563		WriteLn('b = ');
564	564		for i:= 1 to n do
565	565		Write(b[i]:10:0);
566	566		WriteLn;
567	567		WriteLn;
568	568
569	569		// Print result vector c
570	570		WriteLn('c = A b = ');
571	571		for i:= 1 to m do
572	572		Write(c[i]:10:0);
573	573		WriteLn;
574	574		end.
575	575		```
576	576
577	577		### Transpose matrix
578	578
579	579		The [transpose matrix](https://en.wikipedia.org/wiki/Transpose) A<sup>T</sup> of matrix A is obtained by flipping rows and columns:
580	580
581	581
582	582		$`\begin{align*}
583	583		{A} = \left[
584	584		\begin{array}{cccc}
585	585		a_{11} & a_{12} & a_{13} & \ldots & a_{1n}\\
586	586		a_{21} & a_{22} & a_{23} & \ldots & a_{2n}\\
587	587		\vdots & \vdots & \vdots & \ddots & \vdots\\
588	588		a_{m1} & a_{m2} & a_{m3} & \ldots & a_{mn}
589	589		\end{array}
590	590		\right]
591	591		\ \ \
592	592		{A^T} =
593	593		\left[
594	594		\begin{array}{cccc}
595	595		a_{11} & a_{21} & \ldots & a_{m1}\\
596	596		a_{12} & a_{22} & \ldots & a_{m2}\\
597	597		a_{13} & a_{23} & \ldots & a_{m3} \\
598	598		\vdots & \vdots & \ddots & \vdots \\
599	599		a_{1n} & a_{2n} & \ldots & a{mn}
600	600		\end{array}
601	601		\right]
602	602		\end{align*}`$
603	603
604	604		Use the procedure `omvtrm` to perform this operation with NumLib:
605	605
606	606		``` pascal
607	607		procedure omvtrm(
608	608		var a: ArbFloat; m, n, rwa: ArbInt;
609	609		var c: ArbFloat; rwc: ArbInt
610	610		);
611	611		```
612	612
613	613		- `a` denotes the first element of the input matrix A. The elements of this array are not changed by the procedure.
614	614		- `m` is the number of rows of matrix A.
615	615		- `n` is the number of columns of matrix A.
616	616		- `rwa` is the number of allocated columns of A. This way the array of A can be larger than acutally needed.
617	617		- `c` is the first element of the transposed matrix A<sup>T</sup>. It has `n` rows and `m` columns.
618	618		- `rwc` is the number of allocated columns for the transposed matrix.
619	619
620	620		Example
621	621
622	622		Obtain the transpose of the 2x4 matrix
623	623
624	624
625	625		$`{A} = \left[
626	626		\begin{array}{cccc}
627	627		1 & 2 & 3 & 4 \\
628	628		11 & 22 & 33 & 44
629	629		\end{array}
630	630		\right]`$
631	631
632	632		``` pascal
633	633		program transpose;
634	634
635	635		uses
636	636		typ, omv;
637	637
638	638		const
639	639		m = 2;
640	640		n = 4;
641	641
642	642		A: array[1..m, 1..n] of ArbFloat = (
643	643		( 1, 2, 3, 4),
644	644		(11, 22, 33, 44)
645	645		);
646	646
647	647		var
648	648		C: array[1..n, 1..m] of ArbFloat;
649	649		i, j: Integer;
650	650
651	651		begin
652	652		// Transpose A
653	653		omvtrm(
654	654		A[1,1], m, n, n,
655	655		C[1,1], m
656	656		);
657	657
658	658		// Print A
659	659		WriteLn('A = ');
660	660		for i:= 1 to m do begin
661	661		for j := 1 to n do
662	662		Write(A[i, j]:10:0);
663	663		WriteLn;
664	664		end;
665	665		WriteLn;
666	666
667	667		// Print C
668	668		WriteLn('AT = ');
669	669		for i:= 1 to n do begin
670	670		for j := 1 to m do
671	671		Write(C[i, j]:10:0);
672	672		WriteLn;
673	673		end;
674	674		end.
675	675		```
676	676
677	677		### Norm of vectors and matrices
678	678
679	679		A norm is a function which assigns a positive "length" to a [vector a](https://en.wikipedia.org/wiki/Norm_(mathematics)) or a [matrix M](https://en.wikipedia.org/wiki/Matrix_norm).
680	680
681	681		In case of vectors, NumLib supports these norms:
682	682
683	683		- The 1-norm of a vector is defined as the sum of the absolute values of the vector components. It is also called "Taxicab" or "Manhattan" norm because it corresponds to the distance that a taxis has to drive in a rectangular grid of streets.
684	684
685	685
686	686		$`\\|a\\|_1 = \sum_{i=1}^n \|{a_i}\|`$
687	687
688	688		- The 2-norm (also: Euclidian norm) corresponds to the distance of point a from the origin.
689	689
690	690
691	691		$`\\|a\\|_2 = \sqrt{\sum_{i=1}^n {a_i}^2}`$
692	692
693	693		- The maximum infinite norm of a vector is the maximum absolute value of its components.
694	694
695	695
696	696		$`\\|a\\|_\infty = \max({a_1}, {a_2}, ... {a_n})`$
697	697
698	698		The norms of matrices, as implemented by NumLib, are usually defined by row or column vectors:
699	699
700	700		- The 1-norm of a matrix is the maximum of the absolute column sums:
701	701
702	702
703	703		$`\\|M\\|_1 = \max_{1 \le j \le {n}} \sum_{i=1}^m\|M_{ij}\|`$
704	704
705	705		- The maximum infinite norm of a matrix is the maximum of the absolute row sum.
706	706
707	707
708	708		$`\\|M\\|_\infty = \max_{1 \le i \le\ m} \sum_{j=1}^n \|M_{ij}\|`$
709	709
710	710		- Only the Frobenius norm is calculated by individual matrix elements. It corresponds to the 2-norm of vectors, only the sum runs over all matrix elements:
711	711
712	712
713	713		$`\\|M\\|_F = \sqrt{\sum_{i=1}^m \sum_{j=1}^n {M_{ij}}^2}`$
714	714
715	715		These are the NumLib procedure for calculation of norms:
716	716
717	717		``` pascal
718	718		function omvn1v(var a: ArbFloat; n: ArbInt): ArbFloat; // 1-norm of a vector
719	719		function omvn1m(var a: ArbFloat; m, n, rwidth: ArbInt): ArbFloat; // 1-norm of a matrix
720	720		function omvn2v(var a: ArbFloat; n: ArbInt): ArbFloat; // 2-norm of a vector
721	721		function omvnfm(Var a: ArbFloat; m, n, rwidth: ArbInt): ArbFloat; // Frobenius norm of a matrix
722	722		function omvnmv(var a: ArbFloat; n: ArbInt): ArbFloat; // Maximum infinite norm of a vector
723	723		function omvnmm(var a: ArbFloat; m, n, rwidth: ArbInt): ArbFloat; // Maximum infinite norm of a matrix
724	724		```
725	725
726	726		- `a` is the first element of the vector or matrix of which the norm is to be calculated
727	727		- `m`, in case of a matrix norm, is the number of matrix rows.
728	728		- `n` is the number of vector elements in case of a vector norm, or the number of columns in case of a matrix norm.
729	729		- `rwidth`, in case of a matrix norm, is the allocated count of columns. This way a larger matrix can be allocated than actually needed.
730	730
731	731		Example
732	732
733	733		``` pascal
734	734		program norm;
735	735
736	736		uses
737	737		typ, omv;
738	738
739	739		const
740	740		n = 2;
741	741		m = 4;
742	742		a: array[0..m-1] of ArbFloat = (0, 1, 2, -3);
743	743		b: array[0..m-1, 0..n-1] of Arbfloat = ((3, -1), (-2, 2), (0, -1), (2, 1));
744	744
745	745		var
746	746		i, j: Integer;
747	747
748	748		begin
749	749		Write('Vector a = [');
750	750		for i:= Low(a) to High(a) do
751	751		Write(a[i]:5:0);
752	752		WriteLn(' ]');
753	753		WriteLn(' 1-norm of vector a: ', omvn1v(a[0], m):0:3);
754	754		WriteLn(' 2-norm of vector a: ', omvn2v(a[0], m):0:3);
755	755		WriteLn(' maximum norm of vector a: ', omvnmv(a[0], m):0:3);
756	756		WriteLn;
757	757
758	758		WriteLn('Matrix b = ');
759	759		for i:= 0 to m-1 do begin
760	760		for j := 0 to n-1 do
761	761		Write(b[i,j]:5:0);
762	762		WriteLn;
763	763		end;
764	764		WriteLn(' 1-norm of matrix b: ', omvn1m(b[0, 0], m, n, n):0:3);
765	765		WriteLn('Forbenius norm of matrix b: ', omvnfm(b[0, 0], m, n, n):0:3);
766	766		WriteLn(' maximum norm of matrix b: ', omvnmm(b[0, 0], m, n, n):0:3);
767	767
768	768		ReadLn;
769	769		end.
770	770		```
771	771
772	772		## Determinant of a matrix (unit det)
773	773
774	774		Unit det provides a series of routines for calculation of the [determinant](https://en.wikipedia.org/wiki/Determinant) of a square matrix. Besides a procedure for general matrices there are also routines optimized for special matrix types (tridiagonal, symmetric).
775	775
776	776		In order to avoid possible overflows and underflows a factor may be split off during calculation, and the determinant is returned as a floating point number f times a multiplier 8<sup>k</sup>. In many cases, however, k = 0, and the determinant simply is equal to f.
777	777
778	778		### Determinant of a standard matrix
779	779
780	780		The term "standard matrix" means here a matrix with is stored in the standard NumLib storage as described here, i.e. the matrix is stored as a 2D or 1D array of real numbers. The 2D array has at least n rows and n columns, the 1D array is dimensioned to contain at least n<sup>2</sup> elements.
781	781
782	782		Without making any assumptions on the properties of the matrix the determinant can be calculated by the procedure `detgen`. More efficient algorithms can be used, however, if the matrix has these particular properties:
783	783
784	784		- symmetric matrix --\> use procedure `detgsy`.
785	785		- symmetric positive definite matrix --\> use procedure `getgpd`.
786	786
787	787		``` pascal
788	788		procedure detgen(n, rwidth: ArbInt; var a, f: ArbFloat; var k, term: ArbInt); // --> generic matrix
789	789		procedure detgsy(n, rwidth: ArbInt; var a, f: ArbFloat; var k, term: ArbInt); // --> symmetric matrix
790	790		procedure detgpd(n, rwidth: ArbInt; var a, f: ArbFloat; var k, term: ArbInt); // --> symmetric positive definite matrix
791	791		```
792	792
793	793		- `n` determines the number of rows and columns of the input matrix (it must be square matrix).
794	794		- `rwidth` specifies the allocated row length of the input matrix. This takes care of the fact that matrices can be allocated larger than actually needed.
795	795		- `a` is the first (i.e. top/left) element of the input matrix. The entire matrix must be contained within the same allocated memory block.
796	796		- `f` and `k` return the determinant of the matrix. The determinant is given by f 8<sup>k</sup>. Usually k = 0, and the determinant simply is equal to f*.
797	797
798	798		Example
799	799
800	800		Calculate the determinant of the 4x4 matrix
801	801
802	802
803	803		$`A=
804	804		\left[
805	805		\begin{array}{rrrr}
806	806		4 & 2 & 4 & 1 \\
807	807		30 & 20 & 45 & 12 \\
808	808		20 & 15 & 36 & 10 \\
809	809		35 & 28 & 70 & 20
810	810		\end{array}
811	811		\right]`$
812	812
813	813		The matrix is positive definite because all elements are positive and the product b<sup>T</sup> A b cannot become zero or negative for any non-zero vector b. While the routine `detgpd` is the best for this matrix type the code below contains also commented-out calls to the other routines.
814	814
815	815		``` pascal
816	816		program determinant_spd_matrix; // "spd" = symmetric positive definite
817	817
818	818		uses
819	819		math, typ, det;
820	820
821	821		const
822	822		n = 4;
823	823
824	824		const
825	825		a: array[1..n, 1..n] of ArbFloat = (
826	826		( 4, 2, 4, 1)
827	827		(30, 20, 45, 12)
828	828		(20, 15, 36, 10)
829	829		(35, 28, 70, 20)
830	830		);
831	831
832	832		var
833	833		f: ArbFloat;
834	834		k: Integer;
835	835		term: Integer = 0;
836	836		i, j: Integer;
837	837
838	838		begin
839	839		// write input matrix
840	840		WriteLn('A = ');
841	841		for i := 1 to n do begin
842	842		for j := 1 to n do
843	843		Write(a[i, j]:10:0);
844	844		WriteLn;
845	845		end;
846	846		WriteLn;
847	847
848	848		// Calculate determinant of symmetric positive definite matrix
849	849		detgpd(n, n, a[1, 1], f, k, term);
850	850
851	851		// or uncomment to use one of the other methods
852	852		// detgen(n, n, a[1, 1], f, k, term);
853	853		// detgsy(n, n, a[1, 1], f, k, term);
854	854
855	855		if term = 1 then
856	856		WriteLn('Determinant of A = ', f * intpower(8, k):0:9, ' (k = ', k, ')')
857	857		else
858	858		WriteLn('Incorrect input');
859	859		end.
860	860		```
861	861
862	862		### Determinant of a band matrix
863	863
864	864		If a matrix has the shape of a band matrix then its determinant is best calculated by using the procedure `detgba`
865	865
866	866		``` pascal
867	867		procedure detgba(n, l, r: ArbInt; var a, f: ArbFloat; var k, term:ArbInt);
868	868		```
869	869
870	870		- `n` denotes the number of columns and rows of the matrix (it must be square).
871	871		- `l` identifies the left bandwidth, i.e. how many diagonals the band extends to the left, or below, the main diagonal.
872	872		- `r` identifies the right bandwidth, i.e. how many diagons the band extends to the right, or above, the main diagonsl.
873	873		- `a` specifies the first element of a 1D array which contains the elements of the diagonal band, see this. This array contains only the band elements and is obtained by running across the rows of the band matrix from left to right and top to bottom, starting at element A<sub>11</sub>. It must be dimensioned to contain at least `n(l+1+r) - (r(r+1)+l*(l+1)) div 2` elements.
874	874		- `f` and `k` return the determinant of the matrix. The determinant is given by f * 8<sup>k</sup>. Usually k = 0, and the determinant simply is equal to f.
875	875
876	876		Example
877	877
878	878		Calculate the determinant of the 5x5 band matrix with left bandwidth l = 2 and right bandwidth r = 1
879	879
880	880		``` math
881	881		A = \begin{bmatrix}
882	882		1 & 2 & 0 & 0 & 0 \\
883	883		2 & 1 & 1 & 0 & 0 \\
884	884		1 & 3 & 0 & 1 & 0 \\
885	885		0 & 1 & 1 & 1 & 1 \\
886	886		0 & 0 & 1 & 1 & 1
887	887		\end{bmatrix}
888	888		```
889	889
890	890		``` pascal
891	891		program determinant_band_matrix;
892	892
893	893		uses
894	894		math, typ, det;
895	895
896	896		const
897	897		n = 5;
898	898		L = 2;
899	899		R = 1;
900	900		nelements = n(L+1+R) - (R(R+1) + L*(L+1)) div 2;
901	901
902	902		const
903	903		// Diagonal elements arranged in rows, elements outside of the band omitted!
904	904		a: array[1..nelements] of Arbfloat = (
905	905		1, 2,
906	906		2, 1, 1,
907	907		1, 3, 0, 1,
908	908		1, 1, 1, 1,
909	909		1, 1, 1
910	910		);
911	911		var
912	912		f: ArbFloat;
913	913		term: Integer = 0;
914	914		i, i1, i2, j, k: Integer;
915	915
916	916		begin
917	917		// Write diagonal elements
918	918		WriteLn('n = ', n);
919	919		Writeln('Left band width L = ', L);
920	920		WriteLn('Right band width R = ', R);
921	921		Write('Diagonal elements of A = ', a[1]:0:0);
922	922		for k := 2 to nelements do
923	923		Write(a[k]:5:0);
924	924		WriteLn;
925	925		WriteLn;
926	926
927	927		// write reconstructed band input matrix (not needed for calculation)
928	928		WriteLn('Reconstructed A = ');
929	929		k := 1;
930	930		for j := 1 to n do begin
931	931		i1 := j - L;
932	932		i2 := i1 + L + R;
933	933		for i:=1 to i1-1 do
934	934		Write(0.0:5:0);
935	935		for i := i1 to i2 do
936	936		if (i >= 1) and (i <= n) then begin
937	937		Write(a[k]:5:0);
938	938		inc(k);
939	939		end;
940	940		for i:= i2 + 1 to n do
941	941		if i <= n then
942	942		Write(0.0:5:0);
943	943		WriteLn;
944	944		end;
945	945		WriteLn;
946	946
947	947		// Calculate determinant
948	948		detgba(n, l, r, a[1], f, k, term);
949	949		if term = 1 then
950	950		WriteLn('Determinant of A = ', f * intpower(8, k):0:9, ' (k = ', k, ')')
951	951		else
952	952		WriteLn('Incorrect input');
953	953		end.
954	954		```
955	955
956	956		### Determinant of a symmetric positive definite band matrix
957	957
958	958		If a matrix has the shape of a band matrix and is symmetric positive definite then `detgpb` is the dedicated routine for calculation of its determinant.
959	959
960	960		``` pascal
961	961		procedure detgpb(n, w: ArbInt; var a, f: ArbFloat; var k, term:ArbInt);
962	962		```
963	963
964	964		- `n` denotes the number of columns and rows of the matrix (it must be square).
965	965		- `w` identifies the bandwidth, i.e. how many diagonals the band extends to the left or to the right of the main diagonal. Left and right bandwidth must be equal.
966	966		- `a` specifies the first element of a 1D array which contains the elements of the diagonal band, see this. This array contains only the left band elements and the elements of the main diagonal. It is obtained by running across the rows of the band matrix from left to right and top to bottom, but stopping at the main diagonal - the right band is ignored since it is equal to the left band. The array must be dimensioned to contain at least `n(w+1) - (w(w+1)) div 2` elements.
967	967		- `f` and `k` return the determinant of the matrix. The determinant is given by f * 8<sup>k</sup>. Usually k = 0, and the determinant simply is equal to f.
968	968
969	969		Example
970	970
971	971		Calculate the determinant of the 7x7 symmetric band matrix with bandwidth w = 2 $`A= \begin{bmatrix}
972	972		5 & -4 & 1 & & & & \\
973	973		-4 & 6 & -4 & 1 & & & \\
974	974		1 & -4 & 6 & -4 & 1 & & \\
975	975		& 1 & -4 & 6 & -4 & 1 & \\
976	976		& & 1 & -4 & 6 & -4 & 1 \\
977	977		& & & 1 & -4 & 6 & -4 \\
978	978		& & & & 1 & -4 & 5
979	979		\end{bmatrix}`$
980	980
981	981		``` pascal
982	982		program determinant_spd_band_matrix; // "spd" = symmetric positive definite
983	983
984	984		uses
985	985		math, typ, det;
986	986
987	987		const
988	988		n = 7;
989	989		w = 2;
990	990		nelements = n * (w + 1) - (w * (w + 1)) div 2;
991	991
992	992		const
993	993		// Band elements arragned in rows, elements outside of the band omitted!
994	994		a: array[1..nelements] of Arbfloat = (
995	995		5,
996	996		-4, 6,
997	997		1, -4, 6,
998	998		1, -4, 6,
999	999		1, -4, 6,
1000	1000		1, -4, 6,
1001	1001		1, -4, 5
1002	1002		);
1003	1003
1004	1004		function MatrixIndexToArrayIndex(i, j, n, w: Integer): Integer;
1005	1005
1006	1006		function DiagElementIndex(i: Integer): Integer;
1007	1007		var
1008	1008		k: Integer;
1009	1009		begin
1010	1010		Result := 1;
1011	1011		if i = 1 then
1012	1012		exit;
1013	1013		// Rows truncated at left
1014	1014		for k := 2 to w do begin
1015	1015		Result := Result + k;
1016	1016		if k = i then
1017	1017		exit;
1018	1018		end;
1019	1019		// full rows and those truncated at right
1020	1020		for k := w+1 to n do begin
1021	1021		Result := Result + w + 1;
1022	1022		if k = i then
1023	1023		exit;
1024	1024		end;
1025	1025		Result := n;
1026	1026		end;
1027	1027
1028	1028		var
1029	1029		d: Integer;
1030	1030		begin
1031	1031		if j > i then begin
1032	1032		Result := MatrixIndexToArrayIndex(j, i, n, w);
1033	1033		exit;
1034	1034		end;
1035	1035
1036	1036		Result := DiagElementIndex(i);
1037	1037		if (i = j) then
1038	1038		exit;
1039	1039
1040	1040		d := i - j;
1041	1041		if d > w then
1042	1042		Result := -1
1043	1043		else begin
1044	1044		dec(Result, d);
1045	1045		if (Result < 1) then
1046	1046		Result := -1
1047	1047		end;
1048	1048		end;
1049	1049
1050	1050		var
1051	1051		i, j, k: Integer;
1052	1052		f: ArbFloat;
1053	1053		term: Integer = 0;
1054	1054		begin
1055	1055		// Write diagonal elements
1056	1056		WriteLn('n = ', n);
1057	1057		Writeln('(One-sided) band width W = ', w);
1058	1058		Write('Diagonal elements of A = ', a[1]:0:0);
1059	1059		for k := 2 to nelements do
1060	1060		Write(a[k]:3:0);
1061	1061		WriteLn;
1062	1062		WriteLn;
1063	1063
1064	1064		// write reconstructed band input matrix (not needed for calculation)
1065	1065		WriteLn('Reconstructed A = ');
1066	1066		k := 1;
1067	1067		for i := 1 to n do begin
1068	1068		for j := 1 to n do begin
1069	1069		k := MatrixIndexToArrayIndex(i, j, n, w);
1070	1070		if k = -1 then
1071	1071		Write(0.0:3:0)
1072	1072		else
1073	1073		Write(a[k]:3:0);
1074	1074		end;
1075	1075		WriteLn;
1076	1076		end;
1077	1077		WriteLn;
1078	1078
1079	1079		// Calculate determinant
1080	1080		detgpb(n, w, a[1], f, k, term);
1081	1081		if term = 1 then
1082	1082		WriteLn('Determinant of A = ', f * intpower(8, k):0:9, ' (k = ', k, ')')
1083	1083		else
1084	1084		WriteLn('Incorrect input');
1085	1085		end.
1086	1086		```
1087	1087
1088	1088		### Determinant of a tridiagonal matrix
1089	1089
1090	1090		Calculation of the determinant of a tridiagonal matrix is done by the procedure `detgtr`:
1091	1091
1092	1092		``` pascal
1093	1093		procedure detgtr(n: ArbInt; var l, d, u, f: ArbFloat; var k, term:ArbInt);
1094	1094		```
1095	1095
1096	1096		- `n` denotes the number of columns and rows of the matrix (it must be square).
1097	1097		- `l` specifies the first element in the subdiagonal of the matrix. This 1D array must be dimensioned to at least `n-1` elements.
1098	1098		- `d` specifies the first element along the main diagonal of the matrix. This 1D array must be dimensioned to contain at least `n` elements.
1099	1099		- `u` specifies the first element along the superdiagonal of the matrix. The array must be dimensioned to have at least `n-1` elements.
1100	1100		- `f` and `k` return the determinant of the matrix. The determinant is given by f 8<sup>k</sup>. Usually k = 0, and the determinant simply is equal to f*.
1101	1101
1102	1102		Example Calculate the determinant of the matrix
1103	1103
1104	1104
1105	1105		$`A=
1106	1106		\left(
1107	1107		\begin{array}{rrrr}
1108	1108		1 & 1 & & \\
1109	1109		2 & 1 & 1 & \\
1110	1110		& 1 & 1 & 0 \\
1111	1111		& & 1 & 2
1112	1112		\end{array}
1113	1113		\right)`$
1114	1114
1115	1115		``` pascal
1116	1116		program determinant_tridiagonal_matrix;
1117	1117
1118	1118		{$mode objfpc}{$H+}
1119	1119
1120	1120		uses
1121	1121		math, typ, det;
1122	1122
1123	1123		const
1124	1124		n = 4;
1125	1125
1126	1126		const
1127	1127		u: array[0..n-2] of ArbFloat = (-1, -1, -1 ); // superdiagonal
1128	1128		d: array[0..n-1] of ArbFloat = ( 2, 2, 2, 2); // main diagonal
1129	1129		l: array[1..n-1] of ArbFloat = ( -1, -1, -1); // subdiagonal
1130	1130
1131	1131		var
1132	1132		f: ArbFloat;
1133	1133		k: Integer;
1134	1134		term: Integer = 0;
1135	1135		i, j: Integer;
1136	1136
1137	1137		begin
1138	1138		// write input matrix
1139	1139		WriteLn('A = ');
1140	1140		for i := 0 to n-1 do begin
1141	1141		for j := 0 to i-2 do
1142	1142		Write(0.0:5:0);
1143	1143		if i > 0 then
1144	1144		Write(l[i]:5:0);
1145	1145		Write(d[i]:5:0);
1146	1146		if i < n-1 then
1147	1147		Write(u[i]:5:0);
1148	1148		for j:=i+2 to n-1 do
1149	1149		Write(0.0:5:0);
1150	1150		WriteLn;
1151	1151		end;
1152	1152		WriteLn;
1153	1153
1154	1154		// Calculate determinant
1155	1155		detgtr(n, l[1], d[0], u[0], f, k, term);
1156	1156		if term = 1 then
1157	1157		WriteLn('Determinant of A = ', f * intpower(8, k):0:9, ' (k = ', k, ')')
1158	1158		else
1159	1159		WriteLn('Incorrect input');
1160	1160
1161	1161		ReadLn;
1162	1162
1163	1163		end.
1164	1164		```
1165	1165
1166	1166		## Inverse of a matrix (unit inv)
1167	1167
1168	1168		Suppose a square matrix A. Then the matrix A<sup>-1</sup> is the inverse of A if the product A<sup>-1</sup> A is the identity matrix I (all elements zero, except for the diagonal elements which are 1).
1169	1169
1170	1170		The inverse matrix exists if the determinant of A is not zero.
1171	1171
1172	1172
1173	1173		$`A^{-1} A = I =
1174	1174		\begin{bmatrix}
1175	1175		1 & & 0 \\
1176	1176		& \ddots & \\
1177	1177		0 & & 1
1178	1178		\end{bmatrix}`$
1179	1179
1180	1180		NumLib provides various routines for calculation of the inverse matrix. Which algorithm is used depends on the properties of the matrix A. All routines discussed here have in common that the matrix is stored in the standard NumLib way, i.e. as a 2D array of real numbers with n rows and m columns, or as a 1D array with n m elements.
1181	1181
1182	1182		No special requirements are made for the matrix passed to `invgen`, except that it is square. The calculation in this routine is based on [LU decomposition with partial pivoting](https://en.wikipedia.org/wiki/LU_decomposition). If the matrix is symmetric, it can be processed by `invgsy` which uses a reduction to tridiagonal form internally. A symmetric positive definite matrix can be passed to `invgpd` where the inverse is calculated by Cholesky decomposition.
1183	1183
1184	1184		``` pascal
1185	1185		procedure invgen(n, rwidth: ArbInt; var ai: ArbFloat; var term: ArbInt); // --> generic matrix
1186	1186		procedure invgsy(n, rwidth: ArbInt; var ai: ArbFloat; var term: ArbInt); // --> symmetric matrix
1187	1187		procedure invgpd(n, rwidth: ArbInt; var ai: ArbFloat; var term: ArbInt); // --> symmetric positive definite matrix
1188	1188		```
1189	1189
1190	1190		- `n` indicates the size of the matrix.
1191	1191		- `rwidth` is length of the rows of the two-dimensional array which holds the matrix. Normally `n = rwidth`, but the matrix could have been declared larger than actually needed; in this case `rwidth` refers to the declared size of the data array.
1192	1192		- `ai` is the first (top-left) element of the matrix in the data array. After the calculation the input data are overwritten by the elements of the inverse matrix in case of successful completion (`term = 1`); in other cases, the matrix elements are undefined.
1193	1193		- `term` provides an error code:
1194	1194		- 1 - successful completion, the elements of the inverse matrix can be found at the same place as the input matrix.
1195	1195		- 2 - the inverse could not be calculated because the input matrix is (almost) singular.
1196	1196		- 3 - incorrect input data, `n < 1`.
1197	1197
1198	1198		Example
1199	1199
1200	1200		Calculate the inverse of the symmetric matrix $`A =
1201	1201		\left[
1202	1202		\begin{array}{rrrr}
1203	1203		5 & 7 & 6 & 5 \\
1204	1204		7 & 10 & 8 & 7 \\
1205	1205		6 & 8 & 10 & 9 \\
1206	1206		5 & 7 & 9 & 10
1207	1207		\end{array}
1208	1208		\right]`$ This matrix is symmetric and positive definite (the product b<sup>T</sup> A b with any vector b = \[b<sub>1</sub>, b<sub>2</sub>, b<sub>3</sub>\]<sup>T</sup> cannot become zero or negative since all elements of this matrix are positive). Therefore, `invgpd` is best-suited for this task although the other routines can be used as well (uncomment their calls below).
1209	1209
1210	1210		``` pascal
1211	1211		program inverse_matrix;
1212	1212
1213	1213		uses
1214	1214		typ, inv, omv;
1215	1215
1216	1216		const
1217	1217		n = 4;
1218	1218		D = 5;
1219	1219
1220	1220		var
1221	1221		// a is the input matrix to be inverted
1222	1222		// Note that this is matrix must be symmetric positive definite.
1223	1223		a: array[1..n, 1..n] of ArbFloat = (
1224	1224		(5, 7, 6, 5),
1225	1225		(7, 10, 8, 7),
1226	1226		(6, 8, 10, 9),
1227	1227		(5, 7, 9, 10)
1228	1228		);
1229	1229		b: array[1..n, 1..n] of Arbfloat;
1230	1230		c: array[1..n, 1..n] of ArbFloat;
1231	1231		term: Integer = 0;
1232	1232		i, j: Integer;
1233	1233
1234	1234		begin
1235	1235		// write input matrix
1236	1236		WriteLn('a = ');
1237	1237		for i := 1 to n do begin
1238	1238		for j := 1 to n do
1239	1239		Write(a[i, j]:10:D);
1240	1240		WriteLn;
1241	1241		end;
1242	1242		WriteLn;
1243	1243
1244	1244		// Store input matrix because it will be overwritten by the inverse of a
1245	1245		for i := 1 to n do
1246	1246		for j := 1 to n do
1247	1247		b[i, j] := a[i, j];
1248	1248
1249	1249		// Calculate inverse -- uncomment the procedure to be used.
1250	1250		//invgen(n, n, a[1, 1], term);
1251	1251		//invgsy(n, n, a[1, 1], term);
1252	1252		invgpd(n, n, a[1, 1], term);
1253	1253
1254	1254		// write inverse
1255	1255		WriteLn('a^(-1) = ');
1256	1256		for i := 1 to n do begin
1257	1257		for j := 1 to n do
1258	1258		Write(a[i, j]:10:D);
1259	1259		WriteLn;
1260	1260		end;
1261	1261		WriteLn;
1262	1262
1263	1263		// Check validity of result by multiplying inverse with saved input matrix
1264	1264		omvmmm(a[1, 1], n, n, n,
1265	1265		b[1, 1], n, n,
1266	1266		c[1, 1], n);
1267	1267
1268	1268		// write inverse
1269	1269		WriteLn('a^(-1) x a = ');
1270	1270		for i := 1 to n do begin
1271	1271		for j := 1 to n do
1272	1272		Write(c[i, j]:10:D);
1273	1273		WriteLn;
1274	1274		end;
1275	1275		end.
1276	1276		```
1277	1277
1278	1278		## Solving systems of linear equations (unit sle)
1279	1279
1280	1280		A [system of linear equations](https://en.wikipedia.org/wiki/System_of_linear_equations) (or linear system) is a collection of two or more linear equations involving the same set of variables, x<sub>1</sub> ... x<sub>n</sub>.
1281	1281
1282	1282
1283	1283		$`\begin{cases}
1284	1284		a_{11} x_1 + a_{12} x_2 + \dots + a_{1n} x_n = b_1 \\
1285	1285		a_{21} x_1 + a_{22} x_2 + \dots + a_{2n} x_n = b_n \\
1286	1286		\ \ \ \vdots \\
1287	1287		a_{m1} x_1 + a_{m2} x_2 + \dots + a_{mn} x_n = b_m
1288	1288		\end{cases}`$
1289	1289
1290	1290		In order to find a solution which satisfies all equations simultaneously the system of equations is usually converted to a matrix equation A x = b, where
1291	1291
1292	1292
1293	1293		$`A = \left [ \begin{array}{ccc}
1294	1294		a_{11} & a_{12} & \dots & a_{1n} \\
1295	1295		a_{21} & a_{22} & \dots & a_{2n} \\
1296	1296		\vdots & \vdots & \ddots & \vdots \\
1297	1297		a_{m1} & a_{m2} & \dots & a_{mn}
1298	1298		\end{array}
1299	1299		\right ],
1300	1300		\
1301	1301		\mathbf{x} = \left [ \begin{array}{c}
1302	1302		x_1 \\
1303	1303		x_2 \\
1304	1304		\vdots \\
1305	1305		x_n
1306	1306		\end{array}
1307	1307		\right ],
1308	1308		\
1309	1309		\mathbf{b} = \left [ \begin{array}{c}
1310	1310		b_1 \\
1311	1311		b_2 \\
1312	1312		\vdots \\
1313	1313		b_m
1314	1314		\end{array}
1315	1315		\right ]`$
1316	1316
1317	1317		NumLib offers several procedures to solve the matrix equation depending on the properties of the matrix.
1318	1318
1319	1319		### Square matrices with standard storage
1320	1320
1321	1321		Here we describe several routines in which the matrix of the equation system is square, i.e. there are as many equations as unknowns. This matrix is assumed to be passed to the procedures in the standard NumLib way as the first element of a 2D or 1D array. The 2D array must be dimensioned to contain at least n rows and n columns, the 1D array must contain at least n<sup>2</sup> elements.
1322	1322
1323	1323		The most general procedure is `slegen` and applied to any square matrix A. The system is solved using the Gauss elimination procedure with partial pivoting.
1324	1324
1325	1325		If the matrix A is symmetric then the reduction to tridiagonal form generally is more stable. This is done in routine `slegsy`.
1326	1326
1327	1327		If A, furthermore, is a symmetric positive definite matrix, then the Cholesky method is best-suited. This is applied in routine `slegpd`.
1328	1328
1329	1329		``` pascal
1330	1330		procedure slegen(n, rwidth: ArbInt; var a, b, x, ca: ArbFloat; var term: ArbInt); // generic matrix
1331	1331		procedure slegsy(n, rwidth: ArbInt; var a, b, x, ca: ArbFloat; var term: ArbInt); // symmetric matrix
1332	1332		procedure slegpd(n, rwidth: ArbInt; var a, b, x, ca: ArbFloat; var term: ArbInt); // symmetric positive definite matrix
1333	1333		```
1334	1334
1335	1335		- `n` is the number of unknown variables. It must be the same as the number of columns and rows of the coefficent matrix.
1336	1336		- `rwidth` identifies the allocated number of columns of the coefficient matrix. This is useful if the matrix is allocated larger than needed. It is required that `n <= rwidth`.
1337	1337		- `a` is the first (i.e. top left) element of the coefficient matrix A; the other elements must follow within the same allocated memory block. The matrix will not be changed during the calculation.
1338	1338		- `b` is the first element of the array containing the constant vector b. The array length at least must be equal to `n`. The vector will not be changed during the calculation.
1339	1339		- `x` is the first element of the array to receive the solution vector x. It must be allocated to contain at least `n` values.
1340	1340		- `ca` is a parameter to describe the accuracy of the solution.
1341	1341		- `term` returns an error code:
1342	1342		- 1 - successful completion, the solution vector x is valid
1343	1343		- 2 - the solution could not have been determined because the matrix is (almost) singular.
1344	1344		- 3 - error in input values: `n < 1`.
1345	1345
1346	1346		Example
1347	1347
1348	1348		Solve this system of linear equations:
1349	1349
1350	1350
1351	1351		$`\begin{cases}
1352	1352		5 x_1 + 7 x_2 + 6 x_3 + 5 x_4 = 57 \\
1353	1353		7 x_1 + 10 x_2 + 8 x_3 + 7 x_4 = 79 \\
1354	1354		6 x_1 + 8 x_2 + 10 x_3 + 9 x_4 = 88 \\
1355	1355		5 x_1 + 7 x_2 + 9 x_3 + 10 x_4 = 86
1356	1356		\end{cases}
1357	1357		\qquad \Rightarrow \qquad
1358	1358		A=
1359	1359		\left[
1360	1360		\begin{array}{rrrr}
1361	1361		5 & 7 & 6 & 5 \\
1362	1362		7 & 10 & 8 & 7 \\
1363	1363		6 & 8 & 10 & 9 \\
1364	1364		5 & 7 & 9 & 10
1365	1365		\end{array}
1366	1366		\right]
1367	1367		, \;
1368	1368		\mathbf{b} =
1369	1369		\left[
1370	1370		\begin{array}{r}
1371	1371		57 \\
1372	1372		79 \\
1373	1373		88 \\
1374	1374		86
1375	1375		\end{array}
1376	1376		\right]`$
1377	1377
1378	1378		This matrix is symmetric positive definite. Therefore, all routines presented here can be applied. Uncomment the requested line in the code below. As a test the result vector x is multiplied to matrix A, and the result must equal the constant vector b.
1379	1379
1380	1380		``` pascal
1381	1381		program solve_sys_lin_eq;
1382	1382
1383	1383		uses
1384	1384		typ, sle, omv;
1385	1385
1386	1386		const
1387	1387		n = 4;
1388	1388		A: array[1..n, 1..n] of ArbFloat = (
1389	1389		( 5, 7, 6, 5),
1390	1390		( 7, 10, 8, 7),
1391	1391		( 6, 8, 10, 9),
1392	1392		( 5, 7, 9, 10)
1393	1393		);
1394	1394		b: array[1..n] of ArbFloat = (
1395	1395		57, 79, 88, 86
1396	1396		);
1397	1397
1398	1398		var
1399	1399		x: array[1..n] of ArbFloat;
1400	1400		b_test: array[1..n] of ArbFloat;
1401	1401		ca: ArbFloat;
1402	1402		i, j: Integer;
1403	1403		term: Integer;
1404	1404
1405	1405		begin
1406	1406		WriteLn('Solve matrix system A x = b');
1407	1407		WriteLn;
1408	1408
1409	1409		// Print A
1410	1410		WriteLn('Matrix A = ');
1411	1411		for i:= 1 to n do begin
1412	1412		for j := 1 to n do
1413	1413		Write(A[i, j]:10:0);
1414	1414		WriteLn;
1415	1415		end;
1416	1416		WriteLn;
1417	1417
1418	1418		// Print b
1419	1419		WriteLn('Vector b = ');
1420	1420		for i:= 1 to n do
1421	1421		Write(b[i]:10:0);
1422	1422		WriteLn;
1423	1423		WriteLn;
1424	1424
1425	1425		// Solve - select one of these methods
1426	1426		slegen(n, n, A[1, 1], b[1], x[1], ca, term);
1427	1427		// slegsy(n, n, A[1, 1], b[1], x[1], ca, term);
1428	1428		// slegpd(n, n, A[1, 1], b[1], x[1], ca, term);
1429	1429
1430	1430		if term = 1 then begin
1431	1431		WriteLn('Solution vector x = ');
1432	1432		for i:= 1 to n do
1433	1433		Write(x[i]:10:0);
1434	1434		WriteLn;
1435	1435		WriteLn;
1436	1436
1437	1437		omvmmv(A[1,1], n, n, n, x[1], b_test[1]);
1438	1438
1439	1439		WriteLn('Check result: A x = (must be equal to b)');
1440	1440		for i:= 1 to n do
1441	1441		Write(b_test[i]:10:0);
1442	1442		WriteLn;
1443	1443		end
1444	1444		else
1445	1445		WriteLn('Error');
1446	1446
1447	1447		ReadLn;
1448	1448
1449	1449		end.
1450	1450		```
1451	1451
1452	1452		### Band matrix
1453	1453
1454	1454		A special solution, `slegba`, is implemented for a band matrix, i.e. a matrix in which all elements are zero outside a band of width l below and of width r above the main diagonal.
1455	1455
1456	1456		``` pascal
1457	1457		procedure slegba(n, l, r: ArbInt; var a, b, x, ca: ArbFloat; var term:ArbInt);
1458	1458		```
1459	1459
1460	1460		- `n` denotes the number of columns and rows of the matrix (it must be square).
1461	1461		- `l` identifies the left bandwidth, i.e. the number of diagonals the band extends below (or to the left of) the main diagonal.
1462	1462		- `r` identifies the right bandwidth, i.e. the number of diagonals the band extends abofe (or to the right of) the main diagonsl.
1463	1463		- `a` specifies the first element of a 1D array which contains the elements of the diagonal band, see this. This array contains only the band elements and is obtained by running across the rows of the band matrix from left to right and top to bottom, starting at element A<sub>11</sub>. It must be dimensioned to contain at least `n(l+1+r) - (r(r+1)+l*(l+1)) div 2` elements. Note that a 2D array cannot be used for this routine.
1464	1464		- `b` is the first element of the array containing the constant vector b. The array length at least must be equal to `n`. The vector will not be changed during the calculation.
1465	1465		- `x` is the first element of the array to receive the solution vector x. It must be allocated to contain at least `n` values.
1466	1466		- `ca` is a parameter to describe the accuracy of the solution.
1467	1467		- `term` returns an error code:
1468	1468		- 1 - successful completion, the solution vector x is valid
1469	1469		- 2 - the solution could not have been determined because the matrix is (almost) singular.
1470	1470		- 3 - error in input values: `n < 1`, `l < 0`, `l >= n`, `r < 0`, or `r >= n`
1471	1471
1472	1472		Example
1473	1473
1474	1474		Solve this system of linear equations:
1475	1475
1476	1476
1477	1477		$`\begin{array}{ccccccc}
1478	1478		5 x_1 & - & 4 x_2 & + & x_3 & & & & & & & & & = 0 \\
1479	1479		-4 x_1 & + & 6 x_2 & - & 4 x_3 & + & x_4 & & & & & & & = 0 \\
1480	1480		x_1 & - & 4 x_2 & + & 6 x_3 & - & 4 x_4 & & & & & & & = 0 \\
1481	1481		& & x_2 & - & 4 x_3 & + & 6 x_4 & - & 4 x_5 & & & & & = 1 \\
1482	1482		& & & & x_3 & - & 4 x_4 & + & 6 x_5 & - & 4 x_6 & + & x_7 & = 0 \\
1483	1483		& & & & & & x_4 & - & 4 x_5 & + & 6 x_6 & - & 4 x_7 & = 0 \\
1484	1484		& & & & & & & & x_5 & - & 4 x_6 & + & 5 x_7 & = 0 \\
1485	1485		\end{array}
1486	1486		\qquad \Rightarrow \qquad
1487	1487		A=
1488	1488		\left[
1489	1489		\begin{array}{rrrr}
1490	1490		5 & -4 & 1 & 0 & 0 & 0 & 0 \\
1491	1491		-4 & 6 & -4 & 1 & 0 & 0 & 0 \\
1492	1492		1 & -4 & 6 & -4 & 1 & 0 & 0 \\
1493	1493		0 & 1 & -4 & 6 & -4 & 1 & 0 \\
1494	1494		0 & 0 & 1 & -4 & 6 & -4 & 1 \\
1495	1495		0 & 0 & 0 & 1 & -4 & 6 & -4 \\
1496	1496		0 & 0 & 0 & 0 & 1 & -4 & 5
1497	1497		\end{array}
1498	1498		\right]
1499	1499		, \;
1500	1500		\mathbf{b} =
1501	1501		\left[
1502	1502		\begin{array}{r}
1503	1503		0 \\
1504	1504		0 \\
1505	1505		0 \\
1506	1506		1 \\
1507	1507		0 \\
1508	1508		0 \\
1509	1509		0
1510	1510		\end{array}
1511	1511		\right]`$
1512	1512
1513	1513		``` pascal
1514	1514		program solve_band;
1515	1515
1516	1516		uses
1517	1517		typ, sle, omv;
1518	1518
1519	1519		const
1520	1520		n = 7;
1521	1521		L = 2;
1522	1522		R = 2;
1523	1523		nelements = n * (L + 1 + R) - (L * (L + 1) + R * (R + 1)) div 2;
1524	1524
1525	1525		// Band elements arranged in rows, elements outside of the band omitted!
1526	1526		a: array[1..nelements] of Arbfloat = (
1527	1527		5, -4, 1,
1528	1528		-4, 6, -4, 1,
1529	1529		1, -4, 6, -4, 1,
1530	1530		1, -4, 6, -4, 1,
1531	1531		1, -4, 6, -4, 1,
1532	1532		1, -4, 6, -4,
1533	1533		1, -4, 5
1534	1534		);
1535	1535		b: array[1..n] of ArbFloat = (
1536	1536		0, 0, 0, 1, 0, 0, 0
1537	1537		);
1538	1538
1539	1539		var
1540	1540		x: array[1..n] of ArbFloat;
1541	1541		a_test: array[1..n, 1..n] of ArbFloat;
1542	1542		b_test: array[1..n] of ArbFloat;
1543	1543		ca: ArbFloat;
1544	1544		i, j, k: Integer;
1545	1545		term: Integer;
1546	1546
1547	1547		begin
1548	1548		WriteLn('Solve matrix system A x = b where A is a band matrix');
1549	1549		WriteLn;
1550	1550
1551	1551		// Write diagonal elements
1552	1552		WriteLn('Matrix A:');
1553	1553		WriteLn(' n = ', n);
1554	1554		Writeln(' Left band width L = ', L);
1555	1555		WriteLn(' Right band width R = ', R);
1556	1556		WriteLn(' Diagonal elements of A = ');
1557	1557		for k := 1 to nelements do
1558	1558		Write(a[k]:5:0);
1559	1559		WriteLn;
1560	1560		WriteLn;
1561	1561
1562	1562		// Solve
1563	1563		slegba(n, l, r, A[1], b[1], x[1], ca, term);
1564	1564
1565	1565		if term = 1 then begin
1566	1566		// To test the result we multiply A with the solution vector x and check
1567	1567		// whether the result is equal to b.
1568	1568		// Since mutliplication of a matrix stored like a band matrix is quite
1569	1569		// cumbersome we copy the band matrix into a standard matrix.
1570	1570		FillChar(a_test, SizeOf(a_test), 0);
1571	1571		k := 1;
1572	1572		for i:=1 to L do
1573	1573		for j := i-L to i+R do
1574	1574		if (j >= 1) and (j <= n) then begin
1575	1575		a_test[i,j] := A[k];
1576	1576		inc(k);
1577	1577		end;
1578	1578
1579	1579		for i:= L+1 to n-R do
1580	1580		for j := i-L to i+R do begin
1581	1581		a_test[i,j] := A[k];
1582	1582		inc(k);
1583	1583		end;
1584	1584		for i := n-R+1 to n do
1585	1585		for j := i-L to i+R do
1586	1586		if j <= n then begin
1587	1587		a_test[i,j] := A[k];
1588	1588		inc(k);
1589	1589		end;
1590	1590
1591	1591		// Print A
1592	1592		WriteLn('Matrix A =');
1593	1593		for i:=1 to n do begin
1594	1594		for j:=1 to n do
1595	1595		Write(a_test[i,j]:5:0);
1596	1596		WriteLn;
1597	1597		end;
1598	1598		WriteLn;
1599	1599
1600	1600		// Print b
1601	1601		WriteLn('Vector b = ');
1602	1602		for i:= 1 to n do
1603	1603		Write(b[i]:10:0);
1604	1604		WriteLn;
1605	1605		WriteLn;
1606	1606
1607	1607		// Print solution
1608	1608		WriteLn('Solution vector x = ');
1609	1609		for i:= 1 to n do
1610	1610		Write(x[i]:10:3);
1611	1611		WriteLn;
1612	1612		WriteLn;
1613	1613
1614	1614		omvmmv(a_test[1,1], n, n, n, x[1], b_test[1]);
1615	1615
1616	1616		WriteLn('Check result: A x = (must be equal to b)');
1617	1617		for i:= 1 to n do
1618	1618		Write(b_test[i]:10:0);
1619	1619		WriteLn;
1620	1620		end
1621	1621		else
1622	1622		WriteLn('Error');
1623	1623		end.
1624	1624		```
1625	1625
1626	1626		### Symmetric positive definite band matrix
1627	1627
1628	1628		The optimized solution method for a symmetric positive band matrix is called `slegpb`.
1629	1629
1630	1630		``` pascal
1631	1631		procedure slegpb(n, w: ArbInt; var a, b, x, ca: ArbFloat; var term: ArbInt);
1632	1632		```
1633	1633
1634	1634		- `n` denotes the number of columns and rows of the matrix (it must be square).
1635	1635		- `w` identifies the (one-sided) bandwidth, i.e. the number of diagonals the band extends below (or to the left of) the main diagonal. This is the same value as the number of diagonals the band extends above (or to the right of) the main diagonsl.
1636	1636		- `a` specifies the first element of a 1D array which contains the elements of the diagonal band, see this. This array contains only the band elements and is obtained by running across the rows of the band matrix from left to right and top to bottom, starting at element A<sub>11</sub>, but stopping at the main diagonal. This means that the array contains only the left band and the main diagonal elements. It must be dimensioned to contain at least `n(w+1) - (w(w+1)) div 2` elements. Note that a 2D array cannot be used for this routine.
1637	1637		- `b` is the first element of the array containing the constant vector b. The array length at least must be equal to `n`. The vector will not be changed during the calculation.
1638	1638		- `x` is the first element of the array to receive the solution vector x. It must be allocated to contain at least `n` values. If `term` is not 1 then the solution array contains undefined data.
1639	1639		- `ca` is a parameter to describe the accuracy of the solution.
1640	1640		- `term` returns an error code:
1641	1641		- 1 - successful completion, the solution vector x is valid
1642	1642		- 2 - the solution could not have been determined because the matrix is (almost) singular.
1643	1643		- 3 - error in input values: `n<1`, `w<0`, or `w>=n`
1644	1644
1645	1645		Example
1646	1646
1647	1647		The matrix of the band matrix example is symmetrical positive definite. Therefore, procedure `slegpb` can be applied as well:
1648	1648
1649	1649		``` pascal
1650	1650		program solve_spd; // "spd" = symmetric positive definite
1651	1651
1652	1652		uses
1653	1653		typ, sle, omv;
1654	1654
1655	1655		const
1656	1656		n = 7;
1657	1657		w = 2;
1658	1658		nelements = n * (w + 1) - (w * (w + 1)) div 2;
1659	1659
1660	1660		// Band elements arranged in rows, elements outside of the band as well as right diagonals omitted!
1661	1661		a: array[1..nelements] of Arbfloat = (
1662	1662		5,
1663	1663		-4, 6,
1664	1664		1, -4, 6,
1665	1665		1, -4, 6,
1666	1666		1, -4, 6,
1667	1667		1, -4, 6,
1668	1668		1, -4, 5
1669	1669		);
1670	1670		b: array[1..n] of ArbFloat = (
1671	1671		0, 0, 0, 1, 0, 0, 0
1672	1672		);
1673	1673
1674	1674		var
1675	1675		x: array[1..n] of ArbFloat;
1676	1676		a_test: array[1..n, 1..n] of ArbFloat;
1677	1677		b_test: array[1..n] of ArbFloat;
1678	1678		ca: ArbFloat;
1679	1679		i, j, k: Integer;
1680	1680		term: Integer;
1681	1681
1682	1682		begin
1683	1683		WriteLn('Solve matrix system A x = b where A is a symmetric positive definite band matrix');
1684	1684		WriteLn;
1685	1685
1686	1686		// Write diagonal elements
1687	1687		WriteLn('Matrix A:');
1688	1688		WriteLn(' n = ', n);
1689	1689		Writeln(' (One-sided) band width w = ', w);
1690	1690		WriteLn(' Diagonal elements of A = ');
1691	1691		for k := 1 to nelements do
1692	1692		Write(a[k]:5:0);
1693	1693		WriteLn;
1694	1694		WriteLn;
1695	1695
1696	1696		// Solve
1697	1697		slegpb(n, w, A[1], b[1], x[1], ca, term);
1698	1698
1699	1699		if term = 1 then begin
1700	1700		// To test the result we multiply A with the solution vector x and check
1701	1701		// whether the result is equal to b.
1702	1702		// Since mutliplication of a matrix stored like a band matrix is quite
1703	1703		// cumbersome we copy the band matrix into a standard matrix.
1704	1704		FillChar(a_test, SizeOf(a_test), 0);
1705	1705		i := 1;
1706	1706		k := 1;
1707	1707		while (k <= nelements) and (i <= n) do begin
1708	1708		for j := i - w to i do
1709	1709		if (j >= 1) and (j <= n) then begin
1710	1710		a_test[i, j] := a[k];
1711	1711		a_test[j, i] := a[k];
1712	1712		inc(k);
1713	1713		end;
1714	1714		inc(i);
1715	1715		end;
1716	1716
1717	1717		// Print A
1718	1718		WriteLn('Matrix A =');
1719	1719		for i:=1 to n do begin
1720	1720		for j:=1 to n do
1721	1721		Write(a_test[i,j]:5:0);
1722	1722		WriteLn;
1723	1723		end;
1724	1724		WriteLn;
1725	1725
1726	1726		// Print b
1727	1727		WriteLn('Vector b = ');
1728	1728		for i:= 1 to n do
1729	1729		Write(b[i]:10:0);
1730	1730		WriteLn;
1731	1731		WriteLn;
1732	1732
1733	1733		// Print solution
1734	1734		WriteLn('Solution vector x = ');
1735	1735		for i:= 1 to n do
1736	1736		Write(x[i]:10:3);
1737	1737		WriteLn;
1738	1738		WriteLn;
1739	1739
1740	1740		omvmmv(a_test[1,1], n, n, n, x[1], b_test[1]);
1741	1741
1742	1742		WriteLn('Check result: A x = (must be equal to b)');
1743	1743		for i:= 1 to n do
1744	1744		Write(b_test[i]:10:3);
1745	1745		WriteLn;
1746	1746		end
1747	1747		else
1748	1748		WriteLn('Error');
1749	1749		end.
1750	1750		```
1751	1751
1752	1752		### Tridiagonal matrix
1753	1753
1754	1754		NumLib offers two specialized procedures for solving a system of linear equations based on a tridiagonal matrix A. Both rely on the memory-optimized storage of tridiagonal matrices as described here. They differ in the internal calculation algorithm and in the behavior with critical matrices.
1755	1755
1756	1756		``` pascal
1757	1757		procedure sledtr(n: ArbInt; var l, d, u, b, x: ArbFloat; var term: ArbInt);
1758	1758		procedure slegtr(n: ArbInt; var l, d, u, b, x, ca: ArbFloat; var term: ArbInt);
1759	1759		```
1760	1760
1761	1761		- `n` is the number of unknown variables. It must be the same as the number of columns and rows of the coefficent matrix.
1762	1762		- `l` specifies the first element in the subdiagonal of the matrix A. This 1D array must be dimensioned to at least `n-1` elements.
1763	1763		- `d` specifies the first element along the main diagonal of the matrix A. This 1D array must be dimensioned to contain at least `n` elements.
1764	1764		- `u` specifies the first element along the superdiagonal of the matrix A. The array must be dimensioned to have at least `n-1` elements.
1765	1765		- `b` is the first element of the array containing the constant vector b. The array length at least must be equal to `n`. The vector will not be changed during the calculation.
1766	1766		- `x` is the first element of the array to receive the solution vector x. It must be allocated to contain at least `n` values.
1767	1767		- `ca` is a parameter to describe the accuracy of the solution.
1768	1768		- `term` returns an error code:
1769	1769		- 1 - successful completion, the solution vector x is valid
1770	1770		- 2 - the solution could not have been determined because the matrix is (almost) singular.
1771	1771		- 3 - error in input values: `n < 1`.
1772	1772
1773	1773		`sledtr` is numerically stable if matrix A fulfills one of these conditions:
1774	1774
1775	1775		- A is regular (i.e. its inverse matrix exists), and A is columnar-diagonally dominant; this means:
1776	1776		- \\|d<sub>1</sub>\\| ≥ \\|l<sub>2</sub>\\|,
1777	1777		- \\|d<sub>i</sub>\\| ≥ \\|u<sub>i-1</sub>\\| + \\|l<sub>i</sub>+1\\|,i= 2, ...,n''-1,
1778	1778		- \\|d<sub>n</sub>\\| ≥ \\|u<sub>n-1</sub>\\|
1779	1779		- A is regular, and A is diagonally dominant; this means:
1780	1780		- \\|d<sub>1</sub>\\| ≥ \\|l<sub>2</sub>\\|,
1781	1781		- \\|d<sub>i</sub>\\| ≥ \\|u<sub>i</sub>\\| + \\|l<sub>i</sub>\\|,i= 2, ...,n''-1,
1782	1782		- \\|d<sub>n</sub>\\| ≥ \\|u<sub>n</sub>\\|
1783	1783		- A is symmetric and positive-definite.
1784	1784
1785	1785		However, this procedure does not provide the parameter `ca` from which the accuracy of the determined solution can be evaluated. If this is needed the (less stable) procedure `slegtr` must be used.
1786	1786
1787	1787		Example Solve this tridiagonal system of linear equations for n = 8:
1788	1788
1789	1789
1790	1790		$`\begin{cases}
1791	1791		188 x_1 - 100 x_2 = 0 \\
1792	1792		-100 x_1 + 188 x_2 -100 x_3 = 0 \\
1793	1793		\vdots \\
1794	1794		-100 x_{n-2} + 188 x_{n-1} - 100 x_n = 0 \\
1795	1795		-100 x_{n-1} + 188 x_n = 0
1796	1796		\end{cases}
1797	1797		\qquad \Rightarrow \qquad
1798	1798		A=
1799	1799		\left[
1800	1800		\begin{array}{rrrrr}
1801	1801		188 & -100 & & & 0 \\
1802	1802		-100 & 188 & -100 & & \\
1803	1803		& \ddots & \ddots & \ddots & \\
1804	1804		& & -100 & 188 & -100 \\
1805	1805		& & & -100 & 188 \\
1806	1806		\end{array}
1807	1807		\right]
1808	1808		, \;
1809	1809		\mathbf{b} =
1810	1810		\left[
1811	1811		\begin{array}{r}
1812	1812		88 \\
1813	1813		-12 \\
1814	1814		\vdots \\
1815	1815		-12 \\
1816	1816		88
1817	1817		\end{array}
1818	1818		\right]`$
1819	1819
1820	1820		``` pascal
1821	1821		program solve_tridiag;
1822	1822
1823	1823		uses
1824	1824		typ, sle;
1825	1825
1826	1826		const
1827	1827		n = 8;
1828	1828		u: array[1..n-1] of ArbFloat = (-100, -100, -100, -100, -100, -100, -100 );
1829	1829		d: array[1..n] of ArbFloat = ( 188, 188, 188, 188, 188, 188, 188, 188);
1830	1830		l: array[2..n] of ArbFloat = ( -100, -100, -100, -100, -100, -100, -100);
1831	1831		b: array[1..n] of ArbFloat = ( 88, -12, -12, -12, -12, -12, -12, 88);
1832	1832
1833	1833		var
1834	1834		A: array[1..n, 1..n] of ArbFloat;
1835	1835		x: array[1..n] of ArbFloat;
1836	1836		b_test: array[1..n] of ArbFloat;
1837	1837		ca: ArbFloat;
1838	1838		i, j: Integer;
1839	1839		term: Integer;
1840	1840
1841	1841		begin
1842	1842		WriteLn('Solve tridiagonal matrix system A x = b');
1843	1843		WriteLn;
1844	1844
1845	1845		Write('Superdiagonal of A:':25);
1846	1846		for i := 1 to n-1 do
1847	1847		Write(u[i]:10:0);
1848	1848		WriteLn;
1849	1849
1850	1850		Write('Main diagonal of A:':25);
1851	1851		for i:= 1 to n do
1852	1852		Write(d[i]:10:0);
1853	1853		WriteLn;
1854	1854
1855	1855		Write('Subdiagonal of A:':25);
1856	1856		Write('':10);
1857	1857		for i:=2 to n do
1858	1858		Write(l[i]:10:0);
1859	1859		WriteLn;
1860	1860
1861	1861		Write('Vector b:':25);
1862	1862		for i:=1 to n do
1863	1863		Write(b[i]:10:0);
1864	1864		WriteLn;
1865	1865
1866	1866		// Solve
1867	1867		slegtr(n, l[2], d[1], u[1], b[1], x[1], ca, term);
1868	1868		{ or:
1869	1869		sledtr(n, l[2], d[1], u[1], b[1], x[1], term);
1870	1870		}
1871	1871
1872	1872		if term = 1 then begin
1873	1873		Write('Solution vector x:':25);
1874	1874		for i:= 1 to n do
1875	1875		Write(x[i]:10:0);
1876	1876		WriteLn;
1877	1877
1878	1878		// Multiply A with x to test the result
1879	1879		// NumLib does not have a routine to multiply a tridiagonal matrix with a
1880	1880		// vector... Let's do it manually.
1881	1881		b_test[1] := d[1]x[1] + u[1]x[2];
1882	1882		for i := 2 to n-1 do
1883	1883		b_test[i] := l[i]x[i-1] + d[i]x[i] + u[i]*x[i+1];
1884	1884		b_test[n] := l[n]x[n-1] + d[n]x[n];
1885	1885
1886	1886		Write('Check b = A x:':25);
1887	1887		for i:= 1 to n do
1888	1888		Write(b_test[i]:10:0);
1889	1889		WriteLn;
1890	1890		end
1891	1891		else
1892	1892		WriteLn('Error');
1893	1893		end.
1894	1894		```
1895	1895
1896	1896		### Overdetermined systems (least squares)
1897	1897
1898	1898		Unlike the other routines in unit sle which require a square matrix A, `slegls` can solve linear systems described by a rectangular matrix which has more rows than colums, or speaking in terms of equations, has more equations than unknowns. Such a system generally is not solvable exactly. But an approximate solution can be found such that the sum of squared residuals of each equation, i.e. the norm $`\\|\mathbf{b} - A \mathbf{x}\\|_2`$ is as small as possible. The most prominent application of this technique is fitting of equations to data (regression analysis).
1899	1899
1900	1900		``` pascal
1901	1901		procedure slegls(var a: ArbFloat; m, n, rwidtha: ArbInt; var b, x: ArbFloat; var term: ArbInt);
1902	1902		```
1903	1903
1904	1904		- `a` is the first (i.e. top left) element of an array with the coefficient matrix A; the other elements must follow within the same allocated memory block. The array will not be changed during the calculation.
1905	1905		- `m` is the number of rows of the matrix A (or: the number of equations).
1906	1906		- `n` is the number of columns of the matrix A (or: the number of unknown variables). n must not be larger than m.
1907	1907		- `rwidth` identifies the allocated number of columns of the coefficient matrix A. This is useful if the matrix is allocated larger than needed. It is required that `n <= rwidth`.
1908	1908		- `b` is the first element of the array containing the constant vector b. The array length must correspond to the number of matrix rows `m`, but the array can be allocated to be larger.
1909	1909		- `x` is the first element of the array to receive the solution vector x. The array length must correspond to the number of matrix columns `n`, but the array can be allocated to be larger.
1910	1910		- `term` returns an error code:
1911	1911		- 1 - successful completion, the solution vector x is valid
1912	1912		- 2 - there is no unambiguous solution because the columns of the matrix are linearly dependant..
1913	1913		- 3 - error in input values: `n < 1`, or `n` \> `m`.
1914	1914
1915	1915		The method is based on reduction of the array A to upper triangle shape through Householder transformations.
1916	1916
1917	1917		Example
1918	1918
1919	1919		Find the least-squares solution for the system A x = b of 4 equations and 3 unknowns with
1920	1920
1921	1921		$`A=
1922	1922		\left(
1923	1923		\begin{array}{rrr}
1924	1924		1 & 0 & 1 \\
1925	1925		1 & 1 & 1 \\
1926	1926		0 & 1 & 0 \\
1927	1927		1 & 1 & 0
1928	1928		\end{array}
1929	1929		\right), \ \
1930	1930		b=
1931	1931		\left(
1932	1932		\begin{array}{r}
1933	1933		21 \\ 39 \\ 21 \\ 30
1934	1934		\end{array}
1935	1935		\right).`$
1936	1936
1937	1937		``` pascal
1938	1938		program solve_leastsquares;
1939	1939
1940	1940		uses
1941	1941		typ, sle, omv;
1942	1942
1943	1943		const
1944	1944		m = 4;
1945	1945		n = 3;
1946	1946
1947	1947		A: array[1..m, 1..n] of ArbFloat = (
1948	1948		(1, 0, 1),
1949	1949		(1, 1, 1),
1950	1950		(0, 1, 0),
1951	1951		(1, 1, 0)
1952	1952		);
1953	1953		b: array[1..m] of ArbFloat = (
1954	1954		21, 39, 21, 30);
1955	1955
1956	1956		var
1957	1957		x: array[1..n] of ArbFloat;
1958	1958		term: ArbInt;
1959	1959		i, j: Integer;
1960	1960		b_test: array[1..m] of ArbFloat;
1961	1961		sum: ArbFloat;
1962	1962
1963	1963		begin
1964	1964		WriteLn('Solve A x = b with the least-squares method');
1965	1965		WriteLn;
1966	1966
1967	1967		// Display input data
1968	1968		WriteLn('A = ');
1969	1969		for i := 1 to m do
1970	1970		begin
1971	1971		for j := 1 to n do
1972	1972		Write(A[i, j]:10:0);
1973	1973		WriteLn;
1974	1974		end;
1975	1975		WriteLn;
1976	1976		WriteLn('b = ');
1977	1977		for i := 1 to m do
1978	1978		Write(b[i]:10:0);
1979	1979		WriteLn;
1980	1980
1981	1981		// Calculate and show solution
1982	1982		slegls(a[1,1], m, n, n, b[1], x[1], term);
1983	1983
1984	1984		WriteLn;
1985	1985		WriteLn('Solution x = ');
1986	1986		for j:= 1 to n do
1987	1987		Write(x[j]:10:0);
1988	1988		WriteLn;
1989	1989
1990	1990		// Calculate and display residuals
1991	1991		WriteLn;
1992	1992		WriteLn('Residuals A x - b = ');
1993	1993		sum := 0;
1994	1994		omvmmv(a[1,1], m, n, n, x[1], b_test[1]);
1995	1995		for i:=1 to m do begin
1996	1996		Write((b_test[i] - b[i]):10:0);
1997	1997		sum := sum + sqr(b_test[i] - b[i]);
1998	1998		end;
1999	1999		WriteLn;
2000	2000
2001	2001		// Sum of squared residuals
2002	2002		WriteLn;
2003	2003		WriteLn('Sum of squared residuals');
2004	2004		WriteLn(sum:10:0);
2005	2005
2006	2006		WriteLn;
2007	2007		WriteLn('----------------------------------------------------------------------------');
2008	2008		WriteLn;
2009	2009
2010	2010		// Modify solution to show that the sum of squared residuals increases';
2011	2011		WriteLn('Modified solution x'' (to show that it has a larger sum of squared residuals)');
2012	2012		x[1] := x[1] + 1;
2013	2013		x[2] := x[2] - 1;
2014	2014		WriteLn;
2015	2015		for j:=1 to n do
2016	2016		Write(x[j]:10:0);
2017	2017		omvmmv(a[1,1], m, n, n, x[1], b_test[1]);
2018	2018		sum := 0;
2019	2019		for i:=1 to m do
2020	2020		sum := sum + sqr(b_test[i] - b[i]);
2021	2021		WriteLn;
2022	2022		WriteLn;
2023	2023		WriteLn('Sum of squared residuals');
2024	2024		WriteLn(sum:10:0);
2025	2025		end.
2026	2026		```
2027	2027
2028	2028		## Eigenvalues and eigenvectors (unit eig)
2029	2029
2030	2030		A vector is an [eigenvector](https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors) if it does not change its direction after a linear transformation has been applied to it. In matrix terms: if a non-zero vector x of dimension n is multiplied to a square nxn matrix A and the product is the same vector multiplied by a scalar λ then x is an eigenvector, and λ is called an eigenvalue:
2031	2031
2032	2032
2033	2033		$`A \mathbf{x} = \lambda \mathbf{x} \qquad \Rightarrow \qquad A \mathbf{x} - \lambda \mathbf{x} = 0`$
2034	2034
2035	2035		NumLib offers a variety of routines for calculation of eigenvectors and eigenvalues. They are optimized for the matrix type. Every matrix type will be supported by two or four procedures of various degree of complexity; they are identified by appended numbers 1 and 2, or 1 to 4:
2036	2036
2037	2037		- The procedures with appended 1 calculate all eigenvalues.
2038	2038		- The procedures with the appended 2 do the same but only for some eigenvalues which have an index in a specified range (λ<sub>k1</sub>...λ<sub>k2</sub> out of λ<sub>1</sub>...λ<sub>n</sub>).
2039	2039		- The procedures with the appended 3 calculate all eigenvalues and all eigenvectors.
2040	2040		- The procedures with the appended 4 do the same, but again return only some eigenvalues and eigenvectors having an index in the specified interval.
2041	2041
2042	2042		### Matrices with general storage
2043	2043
2044	2044		The routines discussed here assume that the n x n matrix for which the eigenvalues are to be calculated is stored in the conventional way as a 2D (or 1D) array of `ArbFloat` values.
2045	2045
2046	2046		The NumLib method for determining the eigenvalues and eigenvectors of a generic matrix are called `eigge1` and `eigge3` ("ge" = "generic"). Optimized calculation schemes exists for generic symmetric matrices and symmetric positive definite matrices; the corresponding routines are called `eiggs1`, `eiggs2`, `eiggs3`, and `eiggs4` for generic symmetric matrices ("gs" = "generic symmetric")., and `eiggg1`, `eiggg2`, `eiggg3`, `eiggg4` for symmetric positive definite matrices ("gg" = ??)
2047	2047
2048	2048		``` pascal
2049	2049		// Generic matrix (without any special symmetries)
2050	2050		procedure eigge1(var a: ArbFloat; n, rwidth: ArbInt; var lam: complex; var term: ArbInt); // all eigenvalues
2051	2051		procedure eigge3(var a: ArbFloat; n, rwidtha: ArbInt; var lam, x: complex; rwidthx: ArbInt; var term: ArbInt); // all eigenvalues and eigenvectors
2052	2052
2053	2053		// Generic symmetric matrix
2054	2054		procedure eiggs1(var a: ArbFloat; n, rwidth: ArbInt; var lam: ArbFloat; var term: ArbInt); // all eigenvalues
2055	2055		procedure eiggs2(var a: ArbFloat; n, rwidth, k1, k2: ArbInt; var lam: ArbFloat; var term: ArbInt); // some eigenvalues (index k1..k2)
2056	2056		procedure eiggs3(var a: ArbFloat; n, rwidtha: ArbInt; var lam, x: ArbFloat; var term: ArbInt); // all eigenvalues and eigenvectors
2057	2057		procedure eiggs4(var a: ArbFloat; n, rwidtha, k1, k2: ArbInt; var lam, x: ArbFloat; var term: ArbInt); // some eigenvalues and eigenvectors (index k1..k2)
2058	2058
2059	2059		// Symmetric positive definite matrix
2060	2060		procedure eiggg1(var a: ArbFloat; n, rwidth: ArbInt; var lam: ArbFloat; var term: ArbInt); // all eigenvalues
2061	2061		procedure eiggg2(var a: ArbFloat; n, rwidth, k1, k2: ArbInt; var lam: ArbFloat; var term: ArbInt); // some eigenvalues (index k1..k2)
2062	2062		procedure eiggg3(var a: ArbFloat; n, rwidtha: ArbInt; var lam, x: ArbFloat; var term: ArbInt); // all eigenvalues and eigenvectors
2063	2063		procedure eiggg4(var a: ArbFloat; n, rwidtha, k1, k2: ArbInt; var lam, x: ArbFloat; var term: ArbInt); // some eigenvalues and eigenvectors (index k1..k2)
2064	2064		```
2065	2065
2066	2066		- `a` is the first element of an array containing the matrix A for which the eigenvalue/eigenvector has to be calculated. The array must be dimensioned to provide space for at least n<sup>2</sup> floating point values.
2067	2067		- `n` specifies the size of the matrix A, i.e. the number of rows or columns. Note that the input matrix must be square, i.e. the number of rows and columns is equal.
2068	2068		- `rwidth` is the allocated row length of the array `a`. It can be larger than n if the array is allocated larger than necessary, but normally `rwidth = n`.
2069	2069		- `lam` is the first element of an array receiving the calculated eigenvalues. In case of a generic matrix (`eigge1`, `eigge2`) the eigenvalues can be complex; therefore, the array must be dimensioned for values of type `complex` as declared in unit typ. In the other cases (`eiggs1..4` or `eiggg1..4`) the eigenvalues are real, and the array must be dimensioned for datatype `ArbFloat`. Since a nxn matrix has n eigenvalues the array must be allocated for at least `n` values, in case of the procedures with appended 2 or 4 only `k2-k1+1` values are sufficient (see below).
2070	2070		- `term` returns an error code:
2071	2071		- 1 -- successful calculation
2072	2072		- 2 -- calculation failed
2073	2073		- 3 -- error in input data: `n<1`, `k1<1`, `k1>k2`, or `k2>n`.
2074	2074
2075	2075		Additionally, in case of procedure `eigge3`, `eigge4` , `eiggg3` or `eiggg4`:
2076	2076
2077	2077		- `x` is the first element of a matrix to receive the calculated eigenvectors. Again, the eigenvectors of a generic matrix can have complex components. Therefore, the matrix must be declared for the datatype `complex`, and it must be large enough to hold at least n<sup>2</sup> values, in case of `eigge4` or `eiggg4` `n * (k2-k1+1)`. If the matrix is symmetric or positive definite, the eigenvectors are real, and the array must be declared for datatype `ArbFloat`. In any case, the eigenvectors are normalized to unit length and arranged in the columns of this matrix.
2078	2078		- `rwidthx` denotes the allocated row length of the matrix `x`. Thus it is possible to dimension the result matrix larger than actually needed.
2079	2079
2080	2080		Additionally, in case of procedures `eiggs2`, `eiggs4`, `eiggg2` and `eiggg4`:
2081	2081
2082	2082		- `k1` and `k2` define the interval of indexes k for which the eigenvalues (λ<sub>k</sub>) and eigenvalues (c<sub>k</sub>) are to be calculated. They are integers and must be ordered such that `1<=k1<=k2<=n`.
2083	2083
2084	2084		Example
2085	2085
2086	2086		Calculate the eigenvalues and eigenvectors of the matrix
2087	2087
2088	2088
2089	2089		$`A=
2090	2090		\begin{bmatrix}
2091	2091		8 & -1 & -5 \\
2092	2092		-4 & 4 & -2 \\
2093	2093		18 & -5 & -7
2094	2094		\end{bmatrix}`$
2095	2095
2096	2096		``` pascal
2097	2097		program eig_general_matrix;
2098	2098
2099	2099		uses
2100	2100		SysUtils, math,
2101	2101		typ, eig;
2102	2102
2103	2103		const
2104	2104		n = 3;
2105	2105		D = 3;
2106	2106
2107	2107		var
2108	2108		// a is the input matrix
2109	2109		a: array[1..n, 1..n] of ArbFloat = (
2110	2110		( 8, -1, -5),
2111	2111		(-4, 4, -2),
2112	2112		(18, -5, -7)
2113	2113		);
2114	2114		lambda: array[1..n] of complex;
2115	2115		x: array[1..n, 1..n] of complex;
2116	2116		term: Integer = 0;
2117	2117		i, j: Integer;
2118	2118
2119	2119		function ComplexToStr(z: complex; Decimals: Integer): String;
2120	2120		var
2121	2121		sgn: array[boolean] of string = ('+', '-');
2122	2122		begin
2123	2123		Result := Format('%.f %s %.fi', [Decimals, z.Re, sgn[z.Im < 0], Decimals, abs(z.Im)]);
2124	2124		end;
2125	2125
2126	2126		begin
2127	2127		// write input matrix
2128	2128		WriteLn('a = ');
2129	2129		for i := 1 to n do begin
2130	2130		for j := 1 to n do
2131	2131		Write(a[i, j]:10:D);
2132	2132		WriteLn;
2133	2133		end;
2134	2134		WriteLn;
2135	2135
2136	2136		// Calculate eigenvalues/vectors
2137	2137		eigge3(a[1,1], n, n, lambda[1], x[1,1], n, term);
2138	2138
2139	2139		// write eigenvalues
2140	2140		WriteLn('Eigenvalues: lambda = ');
2141	2141		for i := 1 to n do
2142	2142		Write(ComplexToStr(lambda[i], D):25);
2143	2143		WriteLn;
2144	2144		WriteLn;
2145	2145
2146	2146		// Write eigenvectors
2147	2147		WriteLn('Eigenvectors (as columns): x = ');
2148	2148		for i := 1 to n do begin
2149	2149		for j := 1 to n do
2150	2150		Write(ComplexToStr(x[i, j], D):25);
2151	2151		WriteLn;
2152	2152		end;
2153	2153		end.
2154	2154		```
2155	2155
2156	2156		### Symmetric band matrices
2157	2157
2158	2158		NumLib provides four routines to calculate eigenvalues and eigenvectors of symmetric band matrices: `eigbs1`...`eigbs4` ("bs" = band and symmetric):
2159	2159
2160	2160		``` pascal
2161	2161		// All eigenvalues
2162	2162		procedure eigbs1(var a: ArbFloat; n, w: ArbInt; var lam: ArbFloat; var term: ArbInt);
2163	2163
2164	2164		// Some eigenvalues (index k1..k2)
2165	2165		procedure eigbs2(var a: ArbFloat; n, w, k1, k2: ArbInt; var lam: ArbFloat; var term: ArbInt);
2166	2166
2167	2167		// All eigenvalues and all eigenvectors
2168	2168		procedure eigbs3(var a: ArbFloat; n, w: ArbInt; var lam, x: ArbFloat; rwidthx: ArbInt; var term: ArbInt);
2169	2169
2170	2170		// Some eigenvalues, some eigenvectors (index k1..k2)
2171	2171		procedure eigbs4(var a: ArbFloat; n, w, k1, k2: ArbInt; var lam, x: ArbFloat; rwidthx: ArbInt; var m2, term: ArbInt);
2172	2172		```
2173	2173
2174	2174		- `a` specifies the first element of a 1D array which contains the diagonal elements of the left band and the main diagonal. It is obtained by running across the rows of the band matrix from left to right and top to bottom, but stopping at the main diagonal - the right band is ignored since it is equal to the left band. The array must be dimensioned to contain at least `n`(`w`+1) - (`w`(`w`+1)) div 2 elements.
2175	2175		- `n` denotes the number of columns and rows of the matrix (it must be square).
2176	2176		- `w` identifies the bandwidth, i.e. how many diagonals the band extends to the left or to the right of the main diagonal. Left and right bandwidth must be equal because the matrix is symmetric.
2177	2177		- `lam` is the first element of a real-valued array receiving the calculated eigenvalues. The eigenvalues of a symmetric matric are real values. Since an n x n matrix has n eigenvalues the array must be allocated for at least n `ArbFloat` values.
2178	2178		- `x`, in case of routines `eigbs3` and `eigbs4`, is the first element of a matrix to receive the calculated eigenvectors. The matrix must be declared for the datatype `Arbfloat`, and it must be large enough to hold at least n x n values. The eigenvectors are normalized to unit length and are arranged in the columns of this matrix.
2179	2179		- `rwidthx`, in case of routines `eigbs3` and `eigbs4`, denotes the allocated row length of the matrix x. Thus it is possible to dimension the result matrix larger than actually needed.
2180	2180		- `k1` and `k2`, in case of procedures `eigbs2` and `eigbs4`, define the interval of indexes k for which the eigenvalues λ<sub>k</sub> and eigenvalues x<sub>k</sub> are to be calculated. They are integers and must be ordered such that `1<=k1<=k2<=n`.
2181	2181		- `m2`, in case of routine `eigbs4`, indicates the index of the largest eigenvalue for which an eigenvector has been calculated.
2182	2182		- `term` returns an error code:
2183	2183		- 1 -- successful calculation
2184	2184		- 2 -- calculation failed
2185	2185		- 3 -- error in input parameters: `n<1`, `w<0`, `w>=n`, `k1<1`, `k1>k2`, or `k2>n`.
2186	2186
2187	2187		If the bandwidth `w > n/3` then it is better to calculate all eigenvalues and eigenvectors even if not all of them are needed.
2188	2188
2189	2189		Example
2190	2190
2191	2191		Calculate the eigenvalues and eigenvectors of the symmetric 7 x 7 matrix $`A=
2192	2192		\begin{bmatrix}
2193	2193		5 & -4 & 1 & & & & 0 \\
2194	2194		-4 & 6 & -4 & 1 & & & \\
2195	2195		1 & -4 & 6 & -4 & 1 & & \\
2196	2196		& 1 & -4 & 6 & -4 & 1 & \\
2197	2197		& & 1 & -4 & 6 & -4 & 1 \\
2198	2198		& & & 1 & -4 & 6 & -4 \\
2199	2199		0 & & & & 1 & -4 & 5
2200	2200		\end{bmatrix}`$ The eigenvalues are $`\lambda_k = 16\sin^{4}\frac{k\pi}{16}, \qquad k=1,2,\ldots,7`$
2201	2201
2202	2202		``` pascal
2203	2203		program eig_symmband_matrix;
2204	2204
2205	2205		uses
2206	2206		SysUtils, math,
2207	2207		typ, eig;
2208	2208
2209	2209		const
2210	2210		n = 7;
2211	2211		w = 2;
2212	2212		nelements = n * (w + 1) - (w * (w + 1)) div 2;
2213	2213		D = 3;
2214	2214
2215	2215		// a contains the elements of the left and main diagonals of the input matrix
2216	2216		a: array[1..nelements] of ArbFloat = (
2217	2217		5,
2218	2218		-4, 6,
2219	2219		1, -4, 6,
2220	2220		1, -4, 6,
2221	2221		1, -4, 6,
2222	2222		1, -4, 6,
2223	2223		1, -4, 5 );
2224	2224
2225	2225		function MatrixIndexToArrayIndex(i, j, n, w: Integer): Integer;
2226	2226
2227	2227		function DiagElementIndex(i: Integer): Integer;
2228	2228		var
2229	2229		k: Integer;
2230	2230		begin
2231	2231		Result := 1;
2232	2232		if i = 1 then
2233	2233		exit;
2234	2234		// Rows truncated at left
2235	2235		for k := 2 to w do begin
2236	2236		Result := Result + k;
2237	2237		if k = i then
2238	2238		exit;
2239	2239		end;
2240	2240		// full rows and those truncated at right
2241	2241		for k := w+1 to n do begin
2242	2242		Result := Result + w + 1;
2243	2243		if k = i then
2244	2244		exit;
2245	2245		end;
2246	2246		Result := n;
2247	2247		end;
2248	2248
2249	2249		var
2250	2250		d: Integer;
2251	2251		begin
2252	2252		if j > i then begin
2253	2253		Result := MatrixIndexToArrayIndex(j, i, n, w);
2254	2254		exit;
2255	2255		end;
2256	2256
2257	2257		Result := DiagElementIndex(i);
2258	2258		if (i = j) then
2259	2259		exit;
2260	2260
2261	2261		d := i - j;
2262	2262		if d > w then
2263	2263		Result := -1
2264	2264		else begin
2265	2265		dec(Result, d);
2266	2266		if (Result < 1) then
2267	2267		Result := -1
2268	2268		end;
2269	2269		end;
2270	2270
2271	2271		var
2272	2272		lambda: array[1..n] of ArbFloat;
2273	2273		x: array[1..n, 1..n] of ArbFloat;
2274	2274		term: Integer = 0;
2275	2275		i, j, k: Integer;
2276	2276
2277	2277		begin
2278	2278		// Write diagonal elements
2279	2279		WriteLn('n = ', n);
2280	2280		Writeln('(One-sided) band width w = ', w);
2281	2281		Write('Diagonal elements of A = ', a[1]:0:0);
2282	2282		for k := 2 to nelements do
2283	2283		Write(a[k]:3:0);
2284	2284		WriteLn;
2285	2285		WriteLn;
2286	2286
2287	2287		// write reconstructed band input matrix (not needed for calculation)
2288	2288		WriteLn('Reconstructed A = ');
2289	2289		k := 1;
2290	2290		for i := 1 to n do begin
2291	2291		for j := 1 to n do begin
2292	2292		k := MatrixIndexToArrayIndex(i, j, n, w);
2293	2293		if k = -1 then
2294	2294		Write(0.0:3:0)
2295	2295		else
2296	2296		Write(a[k]:3:0);
2297	2297		end;
2298	2298		WriteLn;
2299	2299		end;
2300	2300		WriteLn;
2301	2301
2302	2302		// Calculate eigenvalues/vectors
2303	2303		eigbs3(a[1], n, w, lambda[1], x[1,1], n, term);
2304	2304		if term <> 1 then begin
2305	2305		WriteLn('term = ', term, ' --> ERROR');
2306	2306		halt;
2307	2307		end;
2308	2308
2309	2309		// Write expected results of eigenvalues
2310	2310		WriteLn('Expected eigenvalues:');
2311	2311		for i := 1 to n do
2312	2312		Write(16 * intpower(sin(i*pi/16), 4):15:D);
2313	2313		WriteLn;
2314	2314		WriteLn;
2315	2315
2316	2316		// write eigenvalues
2317	2317		WriteLn('Calculated eigenvalues: lambda = ');
2318	2318		for i := 1 to n do
2319	2319		Write(lambda[i]:15:D);
2320	2320		WriteLn;
2321	2321		WriteLn;
2322	2322
2323	2323		// Write eigenvectors
2324	2324		WriteLn('Eigenvectors (as columns): x = ');
2325	2325		for i := 1 to n do begin
2326	2326		for j := 1 to n do
2327	2327		Write(x[i, j]:15:D);
2328	2328		WriteLn;
2329	2329		end;
2330	2330		end.
2331	2331		```
2332	2332
2333	2333		### Symmetric tridiagonal matrices
2334	2334
2335	2335		This type of matrices is treated best by the routines `eigts1`..`eigts4`:
2336	2336
2337	2337		``` pascal
2338	2338		// All eigenvalues
2339	2339		procedure eigts1(var d, cd: ArbFloat; n: ArbInt; var lam: ArbFloat; var term: ArbInt);
2340	2340
2341	2341		// Some eigenvalues (with indices k1..k2)
2342	2342		procedure eigts2(var d, cd: ArbFloat; n, k1, k2: ArbInt; var lam: ArbFloat; var term: ArbInt);
2343	2343
2344	2344		// All eigenvalues and all eigenvectors
2345	2345		procedure eigts3(var d, cd: ArbFloat; n: ArbInt; var lam, x: ArbFloat; rwidth: ArbInt; var term: ArbInt);
2346	2346
2347	2347		// Some eigenvalues and eigenvectors (with indices k1..k2)
2348	2348		procedure eigts4(var d, cd: ArbFloat; n, k1, k2: ArbInt; var lam, x: ArbFloat; rwidth: ArbInt; var m2, term: ArbInt);
2349	2349		```
2350	2350
2351	2351		- `d` specifies the first element along the main diagonal of the matrix. This 1D array must be dimensioned to contain at least `n` elements.
2352	2352		- `c` specifies the first element in the subdiagonal of the matrix. This 1D array must be dimensioned to at least `n-1` elements.
2353	2353		- `n` denotes the number of columns or rows of the matrix (it must be square).
2354	2354		- `k1` and `k2` define the interval of indexes k for which the eigenvalues (λ<sub>k</sub>) and eigenvalues (c<sub>k</sub>) are to be calculated. They are integers and must be ordered such that `1<=k1<=k2<=n`.
2355	2355		- `lam` is the first element of an array receiving the calculated eigenvalues. Since there are n real eigenvalues the array must be dimensioned to contain at least n values of type `ArbFloat`. In case of the routines with appended 2 and 4, it is sufficient to provide space only for `k2 - k1 + 1` values.
2356	2356		- `x` is the first element of an array to receive the calculated eigenvectors. Since the eigenvectors are real the array must be prepared for data type `ArbFloat`. In case of routine `eigts3`, the array must be dimensioned for at least n<sup>2</sup> elements, in case of routine `eigts4` it must provide space for at least `n * (k2-k1+1)` values. The eigenvectors are normalized to unit length and arranged in the columns of this matrix.
2357	2357		- `rwidthx` denotes the allocated row length of the matrix `x`. Thus it is possible to dimension the result matrix larger than actually needed.
2358	2358		- `term` returns an error code:
2359	2359		- 1 -- successful calculation
2360	2360		- 2 -- calculation failed
2361	2361		- 3 -- error in input data: `n<1`, `k1<1`, `k1>k2`, or `k2>n`.
2362	2362
2363	2363		Example
2364	2364
2365	2365		Calculate eigenvalues and eigenvectors of the matrix $`A=
2366	2366		\begin{bmatrix}
2367	2367		1 & 1 & 0 & 0 \\
2368	2368		1 & 1 & 2 & 0 \\
2369	2369		0 & 2 & 1 & 1 \\
2370	2370		0 & 0 & 1 & 1
2371	2371		\end{bmatrix}`$ The expected eigenvalues are $`-\sqrt{2},\ 2-\sqrt{2},\ \sqrt{2},\ 2+\sqrt{2}`$
2372	2372
2373	2373		``` pascal
2374	2374		program eig_symmtridiag_matrix;
2375	2375
2376	2376		uses
2377	2377		typ, eig;
2378	2378
2379	2379		const
2380	2380		n = 4;
2381	2381
2382	2382		// a contains the elements of the main diagonal of the input matrix
2383	2383		d: array[1..n] of ArbFloat = (1, 1, 1, 1);
2384	2384		// c contains the elements of the subdiagonal
2385	2385		c: array[2..n] of ArbFloat = (1, 2, 1);
2386	2386
2387	2387		var
2388	2388		lambda: array[1..n] of ArbFloat;
2389	2389		x: array[1..n, 1..n] of ArbFloat;
2390	2390		term: Integer = 0;
2391	2391		i, j, k: Integer;
2392	2392
2393	2393		begin
2394	2394		// Write diagonals elements
2395	2395		WriteLn('n = ', n);
2396	2396		Write('Elements of main diagonal = ', d[1]:0:0);
2397	2397		for k := 2 to n do
2398	2398		Write(d[k]:3:0);
2399	2399		WriteLn;
2400	2400		Write('Elements of subdiagonal = ', ' ':3, c[2]:0:0);
2401	2401		for k := 3 to n do
2402	2402		Write(c[k]:3:0);
2403	2403		WriteLn;
2404	2404		WriteLn;
2405	2405
2406	2406		// write reconstructed band input matrix (not needed for calculation)
2407	2407		WriteLn('Reconstructed A = ');
2408	2408		for i := 1 to n do begin
2409	2409		for j := 1 to n do begin
2410	2410		if j = i then
2411	2411		Write(d[i]:3:0)
2412	2412		else if (j = i-1) then
2413	2413		Write(c[i]:3:0)
2414	2414		else if (j = i+1) then
2415	2415		Write(c[i+1]:3:0)
2416	2416		else
2417	2417		Write(0.0:3:0);
2418	2418		end;
2419	2419		WriteLn;
2420	2420		end;
2421	2421		WriteLn;
2422	2422
2423	2423		// Calculate eigenvalues/vectors
2424	2424		eigts3(d[1], c[2], n, lambda[1], x[1,1], n, term);
2425	2425		if term <> 1 then begin
2426	2426		WriteLn('term = ', term, ' --> ERROR');
2427	2427		halt;
2428	2428		end;
2429	2429
2430	2430		// Write expected results of eigenvalues
2431	2431		WriteLn('Expected eigenvalues:');
2432	2432		Write(-sqrt(2):15:3, 2-sqrt(2):15:3, sqrt(2):15:3, 2+sqrt(2):15:3);
2433	2433		WriteLn;
2434	2434		WriteLn;
2435	2435
2436	2436		// write eigenvalues
2437	2437		WriteLn('Calculated eigenvalues: lambda = ');
2438	2438		for i := 1 to n do
2439	2439		Write(lambda[i]:15:3);
2440	2440		WriteLn;
2441	2441		WriteLn;
2442	2442
2443	2443		// Write eigenvectors
2444	2444		WriteLn('Eigenvectors (as columns): x = ');
2445	2445		for i := 1 to n do begin
2446	2446		for j := 1 to n do
2447	2447		Write(x[i, j]:15:3);
2448	2448		WriteLn;
2449	2449		end;
2450	2450		end.
2451	2451		```
2452	2452
2453	2453		## Finding the roots of a function (unit roo)
2454	2454
2455	2455		The roots are the x values at which a function f(x) is zero.
2456	2456
2457	2457		### Roots of a polynomial
2458	2458
2459	2459		A polynomial of degree n
2460	2460
2461	2461		``` math
2462	2462		z^n + a_1 z^{n-1} + a_2 z^{n-2} + ... + a_{n-1} z + a_n = 0
2463	2463		```
2464	2464
2465	2465		always has n, not necessarily distinct, complex solutions. The datatype `complex` has been described in section complex numbers. These roots can be calculated by the function `roopol`:
2466	2466
2467	2467		``` pascal
2468	2468		procedure roopol(var a: ArbFloat; n: ArbInt; var z: complex; var k, term: ArbInt);
2469	2469		```
2470	2470
2471	2471		- `a` is the first element of an array containing the polynomial coefficients. Note that it is assumed that the polynomial has been normalized such that the coefficient of the highest-degree term is 1; this coefficient is not contained in the array. Therefore, the array must be dimensioned for at least n values. Since only real polynomials are handled here the array elements must be of datatype `ArbFloat`. The polynomial coefficients must be ordered from highest to lowest degree. Note that the data organization of the array is different from other polynomial routines in this library.
2472	2472		- `n` is the degree of the polynomial. Must be a positive integer.
2473	2473		- `z` is the first element in an array of `complex` values which returns the roots of the polynomial. The array must be dimensioned to contain at least n values. The returned values are undefined if an error had occurred (`term` \<\> 1).
2474	2474		- `k` returns how many roots were found. Always must be equal to n, otherwise an error has occured.
2475	2475		- `term` returns an error code:
2476	2476		- 1 -- sucessful completion, the result array contains valid data.
2477	2477		- 2 -- Not all roots could be detected, i.e. `k < n`.
2478	2478		- 3 -- Error in input data: `n < 1`.
2479	2479
2480	2480		Example
2481	2481
2482	2482		Calculate the roots of the polynomial z<sup>5</sup> + 3 z<sup>4</sup> + 4 z<sup>3</sup> - 8 z<sup>2</sup>. The expected zero points are:
2483	2483
2484	2484		``` math
2485	2485		z_1=0,\ z_2=0,\ z_3=1,\ z_4=-2+2i,\ z_5=-2-2i
2486	2486		```
2487	2487
2488	2488		``` pascal
2489	2489		program polynomial;
2490	2490
2491	2491		uses
2492	2492		SysUtils, typ, roo;
2493	2493
2494	2494		const
2495	2495		n = 5;
2496	2496		a: array[1..n] of ArbFloat = (3, 4, -8, 0, 0);
2497	2497
2498	2498		var
2499	2499		z: array[1..n] of complex;
2500	2500		k: ArbInt;
2501	2501		term: ArbInt;
2502	2502		i: Integer;
2503	2503		c: complex;
2504	2504
2505	2505		function ComplexToStr(z: complex; Decimals: Integer): string;
2506	2506		const
2507	2507		SIGN: array[boolean] of string = ('+', '-');
2508	2508		begin
2509	2509		Result := Format('%.f %s %.f i', [Decimals, z.re, SIGN[z.im <0], Decimals, abs(z.im)]);
2510	2510		end;
2511	2511
2512	2512		begin
2513	2513		// Solve equation
2514	2514		roopol(a[1], n, z[1], k, term);
2515	2515
2516	2516		if term = 1 then begin
2517	2517		// Display results
2518	2518		WriteLn('Results of procedure roopol:');
2519	2519		for i:=1 to n do
2520	2520		WriteLn(' Solution #', i, ': ', ComplexToStr(z[i], 6):20);
2521	2521		WriteLn;
2522	2522
2523	2523		// Display expected results
2524	2524		Writeln('Expected results:');
2525	2525		c.Init(0, 0); // z1 = 0
2526	2526		WriteLn(' Solution #1: ', complexToStr(c, 6):20);
2527	2527		c.Init(0, 0); // z2 = 0
2528	2528		WriteLn(' Solution #2: ', complexToStr(c, 6):20);
2529	2529		c.Init(1, 0); // z3 = 1
2530	2530		WriteLn(' Solution #3: ', complexToStr(c, 6):20);
2531	2531		c.Init(-2, +2); // z4 = -2 + 2 i
2532	2532		WriteLn(' Solution #4: ', complexToStr(c, 6):20);
2533	2533		c.Init(-2, -2); // z4 = -2 - 2 i
2534	2534		WriteLn(' Solution #5: ', complexToStr(c, 6):20);
2535	2535		end else
2536	2536		WriteLn('ERROR');
2537	2537		end.
2538	2538		```
2539	2539
2540	2540		### Roots of the quadratic equation
2541	2541
2542	2542		The quadratic equation
2543	2543
2544	2544
2545	2545		$`{z}^2 + {p} {z} + {q} = 0`$
2546	2546
2547	2547		is a special polynomal of degree 2. It always has two, not necessarily different, complex roots. These solutions can be determined by the procedure `rooqua`.
2548	2548
2549	2549		``` pascal
2550	2550		procedure rooqua(p, q: ArbFloat; var z1, z2: complex);
2551	2551		```
2552	2552
2553	2553		- `p` and `q` are the (real) coefficients of the quadratic equation
2554	2554		- `z1` and `z2` return the two complex roots. See unit typ for the declaration of the type `complex`
2555	2555
2556	2556		Note that `rooqua` assumes that the coefficient of the quadratic term has been normalized to be 1.
2557	2557
2558	2558		Example
2559	2559
2560	2560		Determine the roots of the equation z<sup>2</sup> + 2 z + 5 = 0.
2561	2561
2562	2562		``` pascal
2563	2563		program quadratic_equ;
2564	2564
2565	2565		uses
2566	2566		SysUtils, typ, roo;
2567	2567
2568	2568		var
2569	2569		z1, z2: complex;
2570	2570
2571	2571		const
2572	2572		SIGN: array[boolean] of string = ('+', '-');
2573	2573
2574	2574		begin
2575	2575		rooqua(2, 5, z1, z2);
2576	2576
2577	2577		WriteLn(Format('1st solution: %g %s %g i', [z1.re, SIGN[z1.im < 0], abs(z1.im)]));
2578	2578		WriteLn(Format('2nd solution: %g %s %g i', [z2.re, SIGN[z2.im < 0], abs(z2.im)]));
2579	2579		end.
2580	2580		```
2581	2581
2582	2582		### Roots of the binomial equation
2583	2583
2584	2584		The equation
2585	2585
2586	2586		``` math
2587	2587		z^n = a
2588	2588		```
2589	2589
2590	2590		is another special polynomial in which all terms except for the highest- and lowest-order terms have zero coefficients. It has n complex solutions which can be calculated by the procedure `roobin`:
2591	2591
2592	2592		``` pascal
2593	2593		procedure roobin(n: ArbInt; a: complex; var z: complex; var term: ArbInt);
2594	2594		```
2595	2595
2596	2596		- `n` is the exponent in the binomial equation used. n must be a positive integer.
2597	2597		- `a` is the constant term a at the right side of the binomial equation. It is expected to be a complex value (see section on complex numbers).
2598	2598		- `z` is the first element of an array of data type `complex` to receive the results of the procedure. It must be dimensioned for at least n `complex` values.
2599	2599		- `term` returns an error code:
2600	2600		- 1 -- successful termination
2601	2601		- 2 -- error in input data: `n < 1`
2602	2602
2603	2603		Example
2604	2604
2605	2605		Calculate the roots of the equation
2606	2606
2607	2607		``` math
2608	2608		z^4 = -1
2609	2609		```
2610	2610
2611	2611		The exact solutions are
2612	2612
2613	2613		``` math
2614	2614		z_1 = \frac{1}{2} \sqrt{2} (1+i) ,\ z_2 = \frac{1}{2} \sqrt{2} (1-i) ,\ z_3= \frac{1}{2} \sqrt{2} (-1+i) ,\ z_4= \frac{1}{2} \sqrt{2} (-1-i)
2615	2615		```
2616	2616
2617	2617		``` pascal
2618	2618		program binomial;
2619	2619
2620	2620		uses
2621	2621		SysUtils, typ, roo;
2622	2622
2623	2623		const
2624	2624		n = 4;
2625	2625
2626	2626		var
2627	2627		z: array[1..n] of complex;
2628	2628		term: ArbInt;
2629	2629		i: Integer;
2630	2630		a: complex;
2631	2631		c: complex;
2632	2632
2633	2633		function ComplexToStr(z: complex; Decimals: Integer): string;
2634	2634		const
2635	2635		SIGN: array[boolean] of string = ('+', '-');
2636	2636		begin
2637	2637		Result := Format('%.g %s %.g i', [Decimals, z.re, SIGN[z.im < 0], Decimals, abs(z.im)]);
2638	2638		end;
2639	2639
2640	2640		begin
2641	2641		// Prepare constant term as a complex value
2642	2642		a.Init(-1, 0);
2643	2643
2644	2644		// Solve equation
2645	2645		roobin(n, a, z[1], term);
2646	2646
2647	2647		if term = 1 then begin
2648	2648		// Display results
2649	2649		WriteLn('Results of procedure roobin:');
2650	2650		for i:=1 to n do
2651	2651		WriteLn(' Solution #', i, ': ', ComplexToStr(z[i], 6):20);
2652	2652		WriteLn;
2653	2653
2654	2654		// Display expected results
2655	2655		Writeln('Expected results:');
2656	2656		c.Init(1, 1);
2657	2657		c.Scale(0.5*sqrt(2));
2658	2658		WriteLn(' Solution #1: ', complexToStr(c, 6):20);
2659	2659		c.Init(1, -1);
2660	2660		c.Scale(0.5*sqrt(2));
2661	2661		WriteLn(' Solution #2: ', complexToStr(c, 6):20);
2662	2662		c.Init(-1, 1);
2663	2663		c.Scale(0.5*sqrt(2));
2664	2664		WriteLn(' Solution #3: ', complexToStr(c, 6):20);
2665	2665		c.Init(-1, -1);
2666	2666		c.Scale(0.5*sqrt(2));
2667	2667		WriteLn(' Solution #4: ', complexToStr(c, 6):20);
2668	2668		end else
2669	2669		WriteLn('ERROR');
2670	2670		end.
2671	2671		```
2672	2672
2673	2673		### Bisection
2674	2674
2675	2675		In the [bisection method](https://en.wikipedia.org/wiki/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.
2676	2676
2677	2677		In NumLib, this approach is supported by the procedure `roof1r`:
2678	2678
2679	2679		``` pascal
2680	2680		procedure roof1r(f: rfunc1r; a, b, ae, re: ArbFloat; var x: ArbFloat; var term: ArbInt);
2681	2681		```
2682	2682
2683	2683		- `f` is the function for which the root is to be determined. It must be a function of one floating point argument (type `ArbFloat`). The type of the function, `rfunc1r`, is declared in unit `typ`.
2684	2684		- `a` and `b` 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.
2685	2685		- `ae` and `re` 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 if `ae` is given as `MachEps` (see unit `typ`). Both parameters must not be negative.
2686	2686		- `x` returns the value of the found root.
2687	2687		- `term` returns whether the process has been successful:
2688	2688		- 1 - successful termination, a zero point has been found with absolute accuracy `ae` or relative accuracy `re`
2689	2689		- 2 - the required accuracy of the root could not be reached; However, the value of x is called the "best achievable" approach
2690	2690		- 3 - error in input parameters: `ae` \< 0 or `re` \< 0, or f(a)*f(b) \> 0
2691	2691
2692	2692		Example
2693	2693
2694	2694		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.
2695	2695
2696	2696		``` pascal
2697	2697		program bisection_demo;
2698	2698
2699	2699		uses
2700	2700		typ, roo;
2701	2701
2702	2702		function f(x: ArbFloat): ArbFloat;
2703	2703		begin
2704	2704		Result := x*x - 2;
2705	2705		end;
2706	2706
2707	2707		var
2708	2708		x: ArbFloat = 0.0;
2709	2709		term : ArbInt;
2710	2710
2711	2711		begin
2712	2712		roof1r(@f, 1.0, 2.0, 1e-9, 0, x, term);
2713	2713		WriteLn('Bisection result ', x);
2714	2714		WriteLn('sqrt(2) ', sqrt(2.0));
2715	2715		end.
2716	2716		```
2717	2717
2718	2718		### Roots of a system of nonlinear equations
2719	2719
2720	2720		Procedure `roofnr` finds the roots of a system of (nonlinear) equations:
2721	2721
2722	2722
2723	2723		$`f_{i}(x_1,x_2,\ldots,x_n)=0, \; i=1,2,\ldots,n`$
2724	2724
2725	2725		``` pascal
2726	2726		procedure roofnr(f: roofnrfunc; n: ArbInt; var x, residu: ArbFloat; ra: ArbFloat; var term: ArbInt);
2727	2727		```
2728	2728
2729	2729		- `f` contains the address of the procedure that calculates the function values f<sub>i</sub>(x<sub>1</sub>, x<sub>2</sub>, ..., x<sub>n</sub>). The type of the function, `roofnrfunc`, is declared in unit typ to be
2730	2730
2731	2731		``` pascal
2732	2732		type
2733	2733		roofnrfunc = procedure(var x, fx: ArbFloat; var deff: boolean);
2734	2734
2735	2735		```
2736	2736
2737	2737		:* Here, `x` is the first element of an array of at least n `ArbFloat` values which provides the x<sub>j</sub> values at which the functions are to be calculated.
2738	2738
2739	2739		:* `fx` is the first element of another array of n `ArbFloat` values to receive the result of each function f<sub>i</sub>. Suppose, we want to solve the following system of equations
2740	2740
2741	2741
2742	2742		$`\begin{array}{l}
2743	2743		f_1(x_1, x_2) = x_1^2 + x_2^2 - 2 = 0 \\
2744	2744		f_2(x_1, x_2) = -(x_1 - 1)^2 + x_2 = 0
2745	2745		\end{array}`$
2746	2746
2747	2747
2748	2748		Then the function to be used as `f` can be written like this:
2749	2749
2750	2750		``` pascal
2751	2751		procedure func(var x1, f1x: real; var deff: boolean);
2752	2752		var
2753	2753		x: array[1..2] of real absolute x1;
2754	2754		f: array[1..2] of real absolute f1x;
2755	2755		begin
2756	2756		f[1] := sqr(x[1]) + sqr(x[2]) - 2;
2757	2757		f[2] := -sqr(x[1] - 1) + x[2];
2758	2758		end;
2759	2759		```
2760	2760
2761	2761		:* The boolean parameter `deff` provides a way to stop the root finding process if some condition is met. The above example has two solutions, namely the intersection of a circle with a parabola. If only the intersection with `x[1] >= 1` is requested then the process can be stopped as soon as a value `x[1] < 1` is reached by setting `deff` to `false`:
2762	2762
2763	2763		``` pascal
2764	2764		procedure func(var x1, f1x: real; var deff: boolean); far;
2765	2765		var
2766	2766		x: array[1..2] of real absolute x1;
2767	2767		f: array[1..2] of real absolute f1x;
2768	2768		begin
2769	2769		deff := x[1] >= 1;
2770	2770		if deff then begin
2771	2771		f[1] := sqr(x[1]) + sqr(x[2]) - 2;
2772	2772		f[2] := -sqr(x[1] - 1) + x[2];
2773	2773		end;
2774	2774		end;
2775	2775		```
2776	2776
2777	2777		- `n` is the number of equations or independent x values (x<sub>i</sub>) in the system. Both must be the same.
2778	2778		- `x` is the first element of an array of at least n `ArbFloat` values which serves for both input and output. On input the array must contain estimates of the roots to be used as starting points of the iterative calculation. On output, the array is overwritten with the zero points found. Inspect the returned value of `term` to decide whether the output array contains a valid solution.
2779	2779		- `residu` is the 2-norm of the vector of residuals in the calculated solution
2780	2780
2781	2781
2782	2782		$`\\|f\\|_2 = \sqrt{\Sigma_{i=1}^{n} f_i^2}`$
2783	2783
2784	2784		- `ra` contains the relative accuracy with which the desired solution is to be calculated. Recommended values: 10<sup>-3</sup>, 10<sup>-5</sup>, 10<sup>-8</sup> if `ArbFloat` is `single`, `real`, or `double`.
2785	2785		- `term` returns an error code:
2786	2786		- 1 -- successfull completion, a solution has been found with the desired accuracy
2787	2787		- 2 -- the desired accuracy could not be reached, the returned values of x are the best ones achievable.
2788	2788		- 3 -- Incorrect input data: `n < 0`, or `re < 0`
2789	2789		- 4 -- The calculation process has been stopped because an internal limit for the count of function calls has been reached. Retry with other starting values for x.
2790	2790		- 5 -- Not enough progess, there may be no solution at all. If the solution is close to 0 choose another starting value.
2791	2791		- 6 -- The procedure wants to calculate a value outside the range defined by `deff`.
2792	2792
2793	2793		Example
2794	2794
2795	2795		Calculate a solutions of the equation system
2796	2796
2797	2797
2798	2798		$`\begin{array}{l}
2799	2799		f_1(x_1, x_2) = x_1^2 + x_2^2 - 2 = 0 \\
2800	2800		f_2(x_1, x_2) = -(x_1 - 1)^2 + x_2 = 0
2801	2801		\end{array}`$
2802	2802
2803	2803		``` pascal
2804	2804		program nonlinsystem;
2805	2805
2806	2806		uses
2807	2807		SysUtils, typ, roo;
2808	2808
2809	2809		const
2810	2810		n = 2;
2811	2811		ra = 1e-10;
2812	2812
2813	2813		var
2814	2814		x: array[1..n] of ArbFloat;
2815	2815		f_check: array[1..n] of ArbFloat;
2816	2816		residu: ArbFloat;
2817	2817		i: Integer;
2818	2818		term: Integer;
2819	2819		deff: Boolean;
2820	2820
2821	2821		procedure funcs(var x0, fx: ArbFloat; var deff: boolean);
2822	2822		var
2823	2823		xloc: array[1..n] of ArbFloat absolute x0;
2824	2824		f: array[1..n] of ArbFloat absolute fx;
2825	2825		x, y, z: ArbFloat;
2826	2826		begin
2827	2827		x := xloc[1];
2828	2828		y := xloc[2];
2829	2829		f[1] := sqr(x) + sqr(y) - 2;
2830	2830		f[2] := -sqr(x - 1) + y;
2831	2831		end;
2832	2832
2833	2833		begin
2834	2834		// Initial guessed values
2835	2835		x[1] := 0;
2836	2836		x[2] := 0;
2837	2837
2838	2838		// Solve equation
2839	2839		roofnr(@funcs, n, x[1], residu, ra, term);
2840	2840
2841	2841		WriteLn('term = ', term);
2842	2842		WriteLn;
2843	2843
2844	2844		if term in [1, 2] then begin
2845	2845		// Display results
2846	2846		WriteLn('Results found by procedure roofnr:');
2847	2847		for i:=1 to n do
2848	2848		WriteLn('Solution #' + IntToStr(i)+': ':20, x[i]:0:6);
2849	2849		WriteLn;
2850	2850
2851	2851		WriteLn('Norm of residuals: ', residu:0:15);
2852	2852		WriteLn;
2853	2853
2854	2854		// Test functions:
2855	2855		f(x[1], f_check[1], deff);
2856	2856		WriteLn('Check if the function values really are 0 at the found roots:');
2857	2857		for i := 1 to n do
2858	2858		WriteLn('Function #' + IntToStr(i) + ': ':20, f_check[i]:0:15);
2859	2859		end else
2860	2860		WriteLn('ERROR');
2861	2861
2862	2862		ReadLn;
2863	2863		end.
2864	2864		```
2865	2865
2866	2866		## Numerical integration of a function (unit int)
2867	2867
2868	2868		The NumLib function `int1fr` calculates the [integral](https://en.wikipedia.org/wiki/Integral) of a given function between limits `a` and `b` with a specified absolute accuracy `ae`:
2869	2869
2870	2870		``` pascal
2871	2871		procedure int1fr(f: rfunc1r; a, b, ae: ArbFloat; var integral, err: ArbFloat; var term: ArbInt);
2872	2872		```
2873	2873
2874	2874		- `f` points to the function to be integrated. It must be a function of a single real variable (type `ArbFloat`) and return an `Arbfloat`. See the type `rfunc1r` declared in unit typ.
2875	2875		- `a, b` are the limits of integration. `a` and/or `b` can attain the values `+/-Infinity` for integrating over an infinite interval. The order of `a` and `b` is handled in the mathematically correct sense.
2876	2876		- `ae` determines the absolute accuracy requested.
2877	2877		- `integreal` returns the value of the integral. It is only valid if `term = 1`.
2878	2878		- `err` returns the achieved accuracy if the specified accuracy could not be reached. `term` has the value 2 in this case.
2879	2879		- `term` returns an error code:
2880	2880		- 1 - successful termination, the integral could be calculated with absolute accuracy `ae`.
2881	2881		- 2 - the requested accuracy could not be reached. But the integral is approximated within the accuracy `err`.
2882	2882		- 3 - incorrect input data: `ae < 0`, or `a = b = infinity`, or `a = b = -infinity`
2883	2883		- 4 - the integral could not be calculated: divergence, or too-slow convergence.
2884	2884
2885	2885		Example
2886	2886
2887	2887		Calculate the integral
2888	2888
2889	2889
2890	2890		$`\int_a^b \frac 1 {x^2} \mathrm{d}x`$
2891	2891
2892	2892		for several integration limits `a` and `b`. Since the function diverges at x = 0 the interval from a to b must not contain this point. The analytic result is $`-1/{b} + 1/{a}`$
2893	2893
2894	2894		``` pascal
2895	2895		program integrate;
2896	2896
2897	2897		uses
2898	2898		SysUtils, typ, int;
2899	2899
2900	2900		function recipx2(x: ArbFloat): ArbFloat;
2901	2901		begin
2902	2902		Result := 1.0 / sqr(x);
2903	2903		end;
2904	2904
2905	2905		function integral_recipx2(a, b: ArbFloat): Arbfloat;
2906	2906		begin
2907	2907		if a = 0 then
2908	2908		a := Infinity
2909	2909		else if a = Infinity then
2910	2910		a := 0.0
2911	2911		else
2912	2912		a := -1/a;
2913	2913		if b = 0 then
2914	2914		b := Infinity
2915	2915		else if b = Infinity then
2916	2916		b := 0.0
2917	2917		else
2918	2918		b := -1/b;
2919	2919		Result := b - a;
2920	2920		end;
2921	2921
2922	2922		procedure Execute(a, b: ArbFloat);
2923	2923		var
2924	2924		err: ArbFloat = 0.0;
2925	2925		term: ArbInt = 0;
2926	2926		integral: ArbFloat = 0.0;
2927	2927		begin
2928	2928		try
2929	2929		int1fr(@recipx2, a, b, 1e-9, integral, err, term);
2930	2930		except
2931	2931		term := 4;
2932	2932		end;
2933	2933		Write(' The integral from ' + FloatToStr(a) + ' to ' + FloatToStr(b));
2934	2934		case term of
2935	2935		1: WriteLn(' is ', integral:0:9, ', exected: ', integral_recipx2(a, b):0:9);
2936	2936		2: WriteLn(' is ', integral:0:9, ', error: ', err:0:9, ', exected: ', integral_recipx2(a, b):0:9);
2937	2937		3: WriteLn(' cannot be calculated: Error in input data');
2938	2938		4: WriteLn(' cannot be calculated: Divergent, or calculation converges too slowly.');
2939	2939		end;
2940	2940		end;
2941	2941
2942	2942		begin
2943	2943		WriteLn('Integral of f(x) = 1/x^2');
2944	2944		Execute(1.0, 2.0);
2945	2945		Execute(1.0, 1.0);
2946	2946		Execute(2.0, 1.0);
2947	2947		Execute(1.0, Infinity);
2948	2948		Execute(-Infinity, -1.0);
2949	2949		Execute(0.0, Infinity);
2950	2950		// NOTE: The next line will raise an exception if run in the IDE. This will not happen outside the IDE.
2951	2951		// Integrate(-1.0, Infinity);
2952	2952		end.
2953	2953		```
2954	2954
2955	2955		## Ordinary differential equations (unit ode)
2956	2956
2957	2957		### Solving a single first-order differential equation
2958	2958
2959	2959		The routine `odeiv1` solves an initial value problem described by a first-order differential equation of the form $`\begin{cases}
2960	2960		y' = f(x, y), \qquad x \in [a, b]\\
2961	2961		y(a) = \alpha
2962	2962		\end{cases}`$ with given a, b, f and initial condition y(a) = α.
2963	2963
2964	2964		``` pascal
2965	2965		procedure odeiv1(f: rfunc2r; a, ya: ArbFloat; var b, yb: ArbFloat; ae: ArbFloat; var term: ArbInt);
2966	2966		```
2967	2967
2968	2968		- `f` is the funtion to be processed. It depends on two real variables and returns a real value, as declared in unit typ as `type rfunc2r = function(x, y: ArbFloat): ArbFloat`.
2969	2969		- `a` is the x value of the starting point of the interval in which the solution is determined.
2970	2970		- `ya` is the value α of the function at the starting point a (initial value).
2971	2971		- `b` is the x value of the end point of the interval in which the solution is determined. The case `a > b` is allowed. If, after the calculation, the error code `term` is 2 then the interval has been changed, and `b` contains the new endpoint to which the solution could be calculated with the desired accuracy, `ae`. In all other cases, `b` is unchanged.
2972	2972		- `yb` returns the result of the calculation if `term` \< 3. It is the value of the function at the end point b of the interval. If `term = 3` then the result is undefined.
2973	2973		- `ae` defines the absolute accuracy with which the value y(b) must be determined.
2974	2974		- `term` returns an error code:
2975	2975		- 1 - successful completion
2976	2976		- 2 - the solution could not reach the point b with the required accuracy. However, `yb` is a sufficiently accurate approximation of y at the delivered `b`.
2977	2977		- 3 - input error: `ae <= 0`
2978	2978
2979	2979		The algorithm is adaptive and is based on an [explicit one-step Runge-Kutta method](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) of order five.
2980	2980
2981	2981		The procedure is not suitable for a rigid differential comparison.
2982	2982
2983	2983		The accuracy can not be guaranteed in all cases. In so-called unstable problems, e.g. small variations in y(a) give large variations in y(b), the error may be significantly larger than what the result of the procedure suggests. Call the procedure with different values of `ae` to investigate this case.
2984	2984
2985	2985		If one wants to solve the initial value problem for a number of points, e.g. from x = 0 to x = 1 with step size 0.1, then it is advisable to "integrate" and thus not to restart at x = 0 with every new step.
2986	2986
2987	2987		Example
2988	2988
2989	2989		Solve the differential equation y" = -10 (y - x<sup>2</sup>) with the initial condition y(0) = 0 for x = 0, 0.5, ..., 5. Compare with the exact solution y(x) = -0.02 exp(-10 x) + x<sup>2</sup> - 0.2 x + 0.02
2990	2990
2991	2991		``` pascal
2992	2992		program solve_ode;
2993	2993
2994	2994		uses
2995	2995		typ, ode;
2996	2996
2997	2997		function f(x, y: ArbFloat): Arbfloat;
2998	2998		begin
2999	2999		Result := -10 * (y - sqr(x));
3000	3000		end;
3001	3001
3002	3002		function exact(x: real): real; far;
3003	3003		begin
3004	3004		Result := -0.02exp(-10x) + sqr(x) - 0.2*x + 0.02;
3005	3005		end;
3006	3006
3007	3007		const
3008	3008		d = 0.5; // interval length
3009	3009		ae = 1e-5;
3010	3010		n = 10;
3011	3011
3012	3012		var
3013	3013		a, b: ArbFloat;
3014	3014		ya, yb: ArbFloat;
3015	3015		term: ArbInt;
3016	3016		i: Arbint;
3017	3017
3018	3018		begin
3019	3019		// Set initial conditions
3020	3020		a := 0.0;
3021	3021		b := a + d;
3022	3022		ya := 0.0;
3023	3023
3024	3024		WriteLn('x':12, 'y':12, 'exact':12, 'error code':17);
3025	3025		WriteLn;
3026	3026		WriteLn(a:12:5, ya:12:5, exact(a):12:5, '-':17);
3027	3027
3028	3028		for i := 1 to n do begin
3029	3029		odeiv1(@f, a, ya, b, yb, ae, term);
3030	3030		WriteLn(b:12:5, yb:12:5, exact(b):12:5, term:17);
3031	3031		a := b;
3032	3032		ya := yb;
3033	3033		b := b + d;
3034	3034		end;
3035	3035		end.
3036	3036		```
3037	3037
3038	3038		### Solving a system of first-order differential equations
3039	3039
3040	3040		An initial value problem described by system of first-order differential equations can be solved by the procedure `odeiv2`: $`\begin{cases}
3041	3041		\mathbf{y}' = \mathbf{f}(x, \mathbf{y}), \qquad x \in [a, b]\\
3042	3042		\mathbf{y}(a) = \mathbf{\alpha}
3043	3043		\end{cases}`$
3044	3044
3045	3045		where
3046	3046
3047	3047		- *y* is an unknown n vector
3048	3048
3049	3049
3050	3050
3051	3051		y: \[a, b\] --\> R<sup>n</sup>
3052	3052
3053	3053		or: y = \[y<sub>1</sub>(x), y<sub>2</sub>(x), ..., y<sub>n</sub>(x)\]
3054	3054
3055	3055		- *f* is a vector function
3056	3056
3057	3057
3058	3058
3059	3059		*f: \[a, b\] --\> R*<sup>n</sup>
3060	3060
3061	3061		or: *f(x, y) = \[f<sub>1</sub>(x,y), f<sub>2</sub>(x,y), ..., f<sub>n</sub>(x,y*)\]
3062	3062
3063	3063		- a, b, and the initial conditions *α* = *y(a) = \[y<sub>1</sub>(a), y<sub>2</sub>(a), ... y<sub>n</sub>(a*)\] are given.
3064	3064
3065	3065		The algorithm is based on an [explicit one-step Runge-Kutta method](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) of order five with variable step size.
3066	3066
3067	3067		``` pascal
3068	3068		procedure odeiv2(f: oderk1n; a: ArbFloat; var ya, b, yb: ArbFloat; n: ArbInt; ae: ArbFloat; var term: ArbInt);
3069	3069		```
3070	3070
3071	3071		- The parameter `f` contains the name of the procedure that calculates the function values *f<sub>i</sub>(x, y<sub>i</sub>). The declaration of the procedure must match the `type oderk1n = procedure(x: ArbFloat; var y, fxy: ArbFloat)` declared in unit typ. Here `x` is the coordinate in the interval \[a, b\], and `y` is the first value of an array of y* values to be used in each equation of the ODE system. The function result is returned in another array which is specified by its first value `fxy`. Both arrays must be dimensioned to provide at least n elements.
3072	3072		- `a` is the begin of the calculation interval.
3073	3073		- `ya` is the first element of an array containing the initial conditions for each ODE of the system, i.e. the function values y<sub>i</sub>(a) = α<sub>i</sub>. The array must be allocated to containt at least n elements.
3074	3074		- `b` is the end of the calculation interval. The case `a > b` is allowed. If after the calculation `term` is 2 then the interval has been changed, and `b` contains the new endpoint to which the solution could have been calculated with the desired accuracy, `ae`. In all other cases, `b` is unchanged.
3075	3075		- `yb` is the first element of an array which receives the results of the calculation. It must be dimensioned to contain at least n elements.
3076	3076		- `n` denotes the number of equations in the system.
3077	3077		- `ae` specifies the absolute accuracy with which the solution must be calculated.
3078	3078		- `term` is an error code:
3079	3079		- 1 - successful completion.
3080	3080		- 2 - the solution could not reach the point b with the required accuracy. However, the values in the `yb` array are sufficiently accurate approximations of the y values at the delivered `b`.
3081	3081		- 3 - input error: `n < 1`, or `ae <= 0`.
3082	3082
3083	3083		Example
3084	3084
3085	3085		Integrate the following initial value problem between x = 0 and x = 1 with step size 0.1: $`\begin{cases}
3086	3086		y'_{1}=2xy_{1}+y_{2} \\
3087	3087		y'_{2}=-y_{1}+2xy_{2}\\
3088	3088		y_1(0)=0,\ y_2(0)=1
3089	3089		\end{cases}`$
3090	3090
3091	3091		The exact solutions are $`\begin{cases}
3092	3092		y_1(x)=\exp({x^2}) \sin(x) \\
3093	3093		y_2(x)=\exp({x^2}) \cos(x)
3094	3094		\end{cases}`$
3095	3095
3096	3096		``` pascal
3097	3097		program solve_ode_sys;
3098	3098
3099	3099		uses
3100	3100		typ, ode;
3101	3101
3102	3102		const
3103	3103		ae = 1e-5;
3104	3104
3105	3105		procedure f(x: ArbFloat; var y, fxy: ArbFloat);
3106	3106		var
3107	3107		_y: array[1..2] of ArbFloat absolute y;
3108	3108		_fxy: array[1..2] of ArbFloat absolute fxy;
3109	3109		begin
3110	3110		_fxy[1] := 2x_y[1] + _y[2];
3111	3111		_fxy[2] := -_y[1] + 2x_y[2]
3112	3112		end;
3113	3113
3114	3114		function exact1(x: ArbFloat): ArbFloat;
3115	3115		begin
3116	3116		Result := exp(xx) sin(x);
3117	3117		end;
3118	3118
3119	3119		function exact2(x: ArbFloat): ArbFloat;
3120	3120		begin
3121	3121		Result := exp(xx) cos(x);
3122	3122		end;
3123	3123
3124	3124		var
3125	3125		a, b, d: ArbFloat;
3126	3126		ya: array[1..2] of ArbFloat;
3127	3127		yb: array[1..2] of ArbFloat;
3128	3128		term: ArbInt;
3129	3129		i: Arbint;
3130	3130		n: ArbInt;
3131	3131
3132	3132		begin
3133	3133		// Set initial conditions
3134	3134		a := 0.0;
3135	3135		b := 0.1;
3136	3136		d := b - a;
3137	3137		ya[1] := 0.0;
3138	3138		ya[2] := 1.0;
3139	3139		n := 10;
3140	3140
3141	3141		WriteLn('x':12, 'y[1]':12, 'y[2]':12, 'exact[1]':12, 'exact[2]':12, 'error code':17);
3142	3142		WriteLn;
3143	3143		WriteLn(a:12:5, ya[1]:12:5, ya[2]:12:5, exact1(a):12:5, exact2(a):12:5, '-':17);
3144	3144
3145	3145		for i := 1 to n do begin
3146	3146		odeiv2(@f, a, ya[1], b, yb[1], 2, ae, term);
3147	3147		WriteLn(b:12:5, yb[1]:12:5, yb[2]:12:5, exact1(b):12:5, exact2(b):12:5, term:17);
3148	3148		a := b;
3149	3149		ya[1] := yb[1];
3150	3150		ya[2] := yb[2];
3151	3151		b := b + d;
3152	3152		end;
3153	3153		end.
3154	3154		```
3155	3155
3156	3156		## Interpolation and fitting (unit ipf)
3157	3157
3158	3158		Unit ipf contains routines for
3159	3159
3160	3160		- fitting a set of data points with a polynomial,
3161	3161		- interpolating or fitting a set of data points with a natural cubic spline. A [spline](https://en.wikipedia.org/wiki/Spline_(mathematics)) is a piecewise defined function of polynomials which has a high degree of smoothness at the connection points of the polynomials ("knots"). In case of a natural cubic spline, the polynomials have degree 3 and their second derivative is zero at the first and last knot.
3162	3162
3163	3163		Fitting means: determine an approximating smooth function such that the deviations from the data points are at minimum.\
3164	3164		Interpolation means: determine a function that runs exactly through the data points.
3165	3165
3166	3166		Use of polynomials or splines is recommended unless the data are known to belong to a known function in which there are still some unknown parameters, for example, the data are measurements of the function ''f(x) = a + b e<sup>-c x</sup>.
3167	3167
3168	3168		In the most common case, the data contain measurement errors and appear to follow the shape of a "simple" function. Therefore, it is recommended to first try to fit with a low-degree polynomial, e.g. not higher than degree 5. If this is not successful or if the shape of the function is more complicated try to fit with a spline.
3169	3169
3170	3170		If the data contain no (or very little) measurement noise, it is almost always wise to interpolate the data with a spline.
3171	3171
3172	3172		### Fitting with a polynomial
3173	3173
3174	3174		For this purpose, the routine `ipfpol` can be applied. It uses given data points (x<sub>i</sub>, y<sub>i</sub>), i = 1, ..., m, and n given coefficients b<sub>0</sub>, ..., b<sub>n</sub> to calculate a polynomial $`\operatorname{p}(x) = b_0 + b_1 x + \dots + b_n x^n`$
3175	3175
3176	3176		In the sense of the least squares fitting, the coefficients b<sub>0</sub> ... b<sub>n</sub> are adjusted to minimize the expression
3177	3177
3178	3178		$`\sum_{i=1}^m {(\operatorname{p}(x) - y_i)^2}`$
3179	3179
3180	3180		``` pascal
3181	3181		procedure ipfpol(m, n: ArbInt; var x, y, b: ArbFloat; var term: ArbInt);
3182	3182		```
3183	3183
3184	3184		- `m` denotes the number of data points. It is required that `m ≥ 1`
3185	3185		- `n` is the degree of the polynomial to be fitted. It is required that `n > 0`.
3186	3186		- `x` is the first element of a one-dimensional array containing the x values of the data points. The array must be allocated to contain at least `m` values.
3187	3187		- `y` is the first element of a one-dimensional array containing the y values of the data points. The array must be dimensioned to contain at least `m` values.
3188	3188		- `b` is the first element of a one-dimenisonal array in which the determined polynomial coefficients are returned if `term` is 1. The array must have been dimensioned to provide space for at least `n + 1` values. The returned coefficients are ordered from lowest to highest degree.
3189	3189		- `term` returns an error code:
3190	3190		- 1 - successful completion, the array `b` contains valid coefficients
3191	3191		- 3 - input error: `n < 0`, or `m < 1`.
3192	3192
3193	3193		After the coefficients b<sub>0</sub>, ..., b<sub>n</sub> once have been determined the value of p(x) can be calculated for any x using the procedure `spepol` from the unit spe.
3194	3194
3195	3195		> Note: With this routine, of course, an interpolating polynomial can also be determined by equaling n to m - 1.
3196	3196
3197	3197		Example
3198	3198
3199	3199		Fit a polynomial of degree 2 to these data points: $`\begin{array}{ccc}
3200	3200		\hline
3201	3201		i & x_i & y_i \\
3202	3202		\hline
3203	3203
3204	3204		1 & 0.00 & 1.26 \\
3205	3205		2 & 0.08 & 1.37 \\
3206	3206		3 & 0.22 & 1.72 \\
3207	3207		4 & 0.33 & 2.08 \\
3208	3208		5 & 0.46 & 2.31 \\
3209	3209		6 & 0.52 & 2.64 \\
3210	3210		7 & 0.67 & 3.12 \\
3211	3211		8 & 0.81 & 3.48
3212	3212		\end{array}`$ ![](/wiki-images/a/ae/Numlib_ipfpol.png)
3213	3213
3214	3214		``` pascal
3215	3215		program polynomial_fit;
3216	3216
3217	3217		uses
3218	3218		SysUtils, StrUtils, typ, ipf, spe;
3219	3219
3220	3220		const
3221	3221		m = 8;
3222	3222		n = 2;
3223	3223		x: array[1..m] of Arbfloat = (
3224	3224		0.00, 0.08, 0.22, 0.33, 0.46, 0.52, 0.67, 0.81
3225	3225		);
3226	3226		y: array[1..m] of Arbfloat = (
3227	3227		1.26, 1.37, 1.72, 2.08, 2.31, 2.64, 3.12, 3.48
3228	3228		);
3229	3229
3230	3230		var
3231	3231		b: array[0..n] of ArbFloat;
3232	3232		i: Integer;
3233	3233		term: ArbInt;
3234	3234		xint, yint: ArbFloat;
3235	3235		begin
3236	3236		WriteLn('Fitting a polynomial of degree 2: p(x) = b[0] + b[1] x + b[2] x^2');
3237	3237		WriteLn;
3238	3238
3239	3239		WriteLn('Data points');
3240	3240		WriteLn('i':10, 'xi':10, 'yi':10);
3241	3241		for i:=1 to m do
3242	3242		WriteLn(i:10, x[i]:10:2, y[i]:10:2);
3243	3243		WriteLn;
3244	3244
3245	3245		// Execute fit
3246	3246		ipfpol(m, n, x[1], y[1], b[0], term);
3247	3247
3248	3248		// Display fitted polynomial coefficients
3249	3249		WriteLn('Fitted coefficients b = ');
3250	3250		WriteLn(b[0]:10:6, b[1]:10:6, b[2]:10:6);
3251	3251		WriteLn;
3252	3252
3253	3253		// Interpolate and display fitted polynomial for some x values
3254	3254		WriteLn('Some interpolated data');
3255	3255		WriteLn('':10, 'x':10, 'y':10);
3256	3256		for i:= 1 to 5 do begin
3257	3257		xint := 0.2*i;
3258	3258		yint := spepol(xint, b[0], n);
3259	3259		WriteLn('':10, xint:10:2, yint:10:2);
3260	3260		end;
3261	3261		end.
3262	3262		```
3263	3263
3264	3264		### Interpolation with a natural cubic spline
3265	3265
3266	3266		The routine `ipfisn` can be used for this purpose. It calculates the parameters of the spline. Once these parameters are known, the value of the spline can be calculated at any point using the `ipfspn` procedure.
3267	3267
3268	3268		Assume that (x<sub>i</sub>, y<sub>i</sub>),i= 0, ...,n' are the given data points with x<sub>0</sub> \< x<sub>1</sub> \< ... \< x<sub>n</sub>. Then assume a function g(x) which is twice continuously differentiable (i.e., the first and second derivatives exist and are continuous) and which interpolates the data (i.e., g(x<sub>i</sub> = y<sub>i</sub>). Imagine g(x) as the shape of a thin flexible bar. Then the bending energy (in first approximation) is proportional to
3269	3269
3270	3270		$`\int_{-\infty}^{\infty} g''(x)^2 dx`$
3271	3271
3272	3272		The interpolating natural spline s can be understood as the function g which minimizes this integral and therefore the bending energy.
3273	3273
3274	3274		The spline s has the following properties:
3275	3275
3276	3276		- s is a cubic polynomial on the interval \[x<sub>i</sub>, x<sub>i+1</sub>\]
3277	3277		- it is twice continuously differentiable on (-∞, ∞)
3278	3278		- s"(x<sub>0</sub>) = s"(x<sub>n</sub>) = 0
3279	3279		- on (-∞, x<sub>0</sub>) and (x<sub>n</sub>, ∞) s is linear.
3280	3280
3281	3281		``` pascal
3282	3282		procedure ipfisn(n: ArbInt; var x, y, d2s: ArbFloat; var term: ArbInt);
3283	3283		```
3284	3284
3285	3285		- Calculates the second derivatives of the splines at the knots x<sub>i</sub>.
3286	3286		- `n` is the index of the last data point, i.e. there are `n+1` data points.
3287	3287		- `x` is the first element of the array with the x coordinates of the data points. It must be dimensioned to contain at last `n+1` values. The standard declaration would be `var x: array[0..n] of ArbFloat`. The data points must be in strictly ascending order, i.e. there must not be any data points for which x<sub>i</sub> \>= x<sub>i+1</sub>.
3288	3288		- `y` is the first element of the array with the y coordinates of the data points. It must be dimensioned to contain at last `n+1` values. The standard declaration would be `var y: array[0..n] of ArbFloat`.
3289	3289		- `d2s` is the first element of the array which receives the second derivatives of the splines at the knots if the returned error code `term` is 1. Note that the array does not contain values for the first and last data point because they are automatically set to 0 (natural splines). Therefore the array must be dimensioned to contain at least `n-1` elements. The standard declaration would be: `var d2s: array[1..n-1] of ArbFloat`
3290	3290		- `term` returns an error code:
3291	3291		- 1 - successful completion, the array beginning at `d2s` contains valid data
3292	3292		- 2 - calculation of the second derivatives was not successful
3293	3293		- 3 - incorrect input parameters: n ≤ 1, or x<sub>i</sub> ≥ x<sub>i+1</sub>
3294	3294
3295	3295		``` pascal
3296	3296		function ipfspn(n: ArbInt; var x, y, d2s: ArbFloat; t: ArbFloat; var term: ArbInt): ArbFloat;
3297	3297		```
3298	3298
3299	3299		- Calculates the value of the spline function at any coordinate t. The calculation uses the second derivatives determined by `ipnisn`.
3300	3300		- `n`, `x`, `y`, `d2s` have the same meaning as above
3301	3301		- `term` returns an error code:
3302	3302		- 1 - successful completion
3303	3303		- 3 - The calculation could not be executed because n ≤ 1.
3304	3304
3305	3305		Example
3306	3306
3307	3307		Interpolate a natural cubic spline through the following data points: $`\begin{array}{ccc}
3308	3308		\hline
3309	3309		i & x_i & y_i \\
3310	3310		\hline
3311	3311
3312	3312		0 & 0.00 & 0.980 \\
3313	3313		1 & 0.09 & 0.927 \\
3314	3314		2 & 0.22 & 0.732 \\
3315	3315		3 & 0.34 & 0.542 \\
3316	3316		4 & 0.47 & 0.385 \\
3317	3317		5 & 0.58 & 0.292 \\
3318	3318		6 & 0.65 & 0.248 \\
3319	3319		7 & 0.80 & 0.179 \\
3320	3320		8 & 0.93 & 0.139
3321	3321		\end{array}`$ ![](/wiki-images/e/e6/Numlib_ipfisn.png)
3322	3322
3323	3323		``` pascal
3324	3324		program spline_interpolation;
3325	3325
3326	3326		uses
3327	3327		typ, ipf;
3328	3328
3329	3329		const
3330	3330		n = 8;
3331	3331		x: array[0..n] of Arbfloat = (
3332	3332		0.00, 0.09, 0.22, 0.34, 0.47, 0.58, 0.65, 0.80, 0.93
3333	3333		);
3334	3334		y: array[0..n] of Arbfloat = (
3335	3335		0.990, 0.927, 0.734, 0.542, 0.388, 0.292, 0.248, 0.179, 0.139
3336	3336		);
3337	3337
3338	3338		var
3339	3339		s: array[0..n] of ArbFloat;
3340	3340		d2s: array[1..n-1] of ArbFloat;
3341	3341		i: Integer;
3342	3342		term: ArbInt;
3343	3343		xint, yint: ArbFloat;
3344	3344		begin
3345	3345		WriteLn('Interpolation with a natural cubic spline');
3346	3346		WriteLn;
3347	3347
3348	3348		WriteLn('Data points');
3349	3349		WriteLn('i':10, 'xi':10, 'yi':10);
3350	3350		for i:=0 to n do
3351	3351		WriteLn(i:10, x[i]:10:2, y[i]:10:3);
3352	3352		WriteLn;
3353	3353
3354	3354		// Interpolation
3355	3355		ipfisn(n, x[0], y[0], d2s[1], term);
3356	3356
3357	3357		// Display 2nd derivatives of splines at xi
3358	3358		WriteLn('s"(xi)');
3359	3359		Write(0.0:10:6); // 2nd derivative of 1st point is 0
3360	3360		for i:=1 to n-1 do
3361	3361		Write(d2s[i]:10:6);
3362	3362		WriteLn(0.0:10:6); // 2nd derivative of last point is 0
3363	3363		WriteLn;
3364	3364
3365	3365		// Interpolate for some x values
3366	3366		WriteLn('Some interpolated data');
3367	3367		WriteLn('':10, 'x':10, 'y':10);
3368	3368		for i:= 1 to 5 do begin
3369	3369		xint := 0.2 * i;
3370	3370		yint := ipfspn(n, x[0], y[0], d2s[1], xint, term);
3371	3371		WriteLn('':10, xint:10:2, yint:10:2);
3372	3372		end;
3373	3373
3374	3374		end.
3375	3375		```
3376	3376
3377	3377		### Fitting with a natural cubic spline
3378	3378
3379	3379		The routine `ipffsn` can be used for this purpose. This routine calculates the spline-fit parameters. Once these parameters are known the value of the spline can be calculated at any point using the procedure `ipfspn`.
3380	3380
3381	3381		Following the model with the bending energy in the previous chapter the calculated fit minimizes
3382	3382
3383	3383		$`\int_{-\infty}^{\infty} g''(x)^2 dx + \lambda \sum_{i=0}^n (g(x_i) - y_i)^2`$
3384	3384
3385	3385		where the parameter λ must be determined. It is a measure of the compromise between the "smoothness" of the fit - which is measured by the integral (the bending energy) - and the deviation of the fit to the given data - which is measured by the sum of squares. The solution is again a natural cubic spline.
3386	3386
3387	3387		> Warning: The procedure `ipffsn` is buggy and crashes. Therefore, it is not documented at the moment.
3388	3388
3389	3389		## Special functions (unit spe)
3390	3390
3391	3391		### Evaluation of a polynomial
3392	3392
3393	3393		[Horner's scheme](https://en.wikipedia.org/wiki/Horner%27s_method) provides an efficient method for evaluating a [polynomial](https://en.wikipedia.org/wiki/Polynomial) at a specific x value:
3394	3394
3395	3395
3396	3396		$`\operatorname{p}(x) = a_0 + a_1 x + a_2 x^2 + ... + a_n x^n \ \Rightarrow \ \operatorname{p}(x) = a_0 + x (a_1 + x (a_2 + \dots + x (a_{n-1} + x a_n)))`$
3397	3397
3398	3398		This method is applied in NumLib's procedure `spepol`:
3399	3399
3400	3400		``` pascal
3401	3401		function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;
3402	3402		```
3403	3403
3404	3404		- `x` is the argument at which the polynomial is evaluated.
3405	3405		- `a` is the first element of an array containing the polynomial coefficients a<sub>0</sub>, a<sub>1</sub>, a<sub>2</sub>, ... a<sub>n</sub>. The coefficient must be ordered from lowest to highest degree.
3406	3406		- `n` is the degree of the polynomial, i.e. there are `n+1` elements in the array.
3407	3407
3408	3408		Example
3409	3409
3410	3410		Evaluate the polynomial
3411	3411
3412	3412
3413	3413		$`\operatorname{p}(x) = 2 x^3 - 6 x^2 + 2 x - 1`$
3414	3414
3415	3415		at x = 3. Note that the terms are not in the correct order for `spepol` and rearrange them to
3416	3416
3417	3417
3418	3418		$`\operatorname{p}(x) = -1 + 2 x - 6 x^2 + 2 x^3`$
3419	3419
3420	3420		``` pascal
3421	3421		program polynomial;
3422	3422
3423	3423		uses
3424	3424		typ, spe;
3425	3425
3426	3426		const
3427	3427		n = 3;
3428	3428		a: array[0..n] of ArbFloat = (-1, 2, -6, 2);
3429	3429
3430	3430		var
3431	3431		x, y: ArbFloat;
3432	3432		begin
3433	3433		x := 3.0;
3434	3434		y := spepol(x, a[0], n);
3435	3435		WriteLn(y:0:6);
3436	3436		end.
3437	3437		```
3438	3438
3439	3439		### Error function
3440	3440
3441	3441		![](/wiki-images/2/20/Numlib_erf_erfc.png) The [error function erf(x)](https://en.wikipedia.org/wiki/Error_function) and its complement, the [complementary error function erfc(x)](https://en.wikipedia.org/wiki/Error_function#Complementary_error_function), are the lower and upper integrals of the Gaussian function, normalized to unity:
3442	3442
3443	3443
3444	3444		$`\operatorname{erf}(x) = \frac {2} {\sqrt{\pi}} \int_0^x e^{-{t}^2} dt`$
3445	3445
3446	3446
3447	3447		$`\operatorname{erfc}(x) = \frac {2} {\sqrt{\pi}} \int_x^\infty e^{-{t}^2} dt = {1} - \operatorname{erf}(x)`$
3448	3448
3449	3449		Both functions can be calculated by the NumLib functions `speerf` and `speefc`, respectively
3450	3450
3451	3451		``` pascal
3452	3452		function speerf(x: ArbFloat): ArbFloat; // --> erf(x)
3453	3453		function speefc(x: ArbFloat): ArbFloat; // --> erfc(x)
3454	3454		```
3455	3455
3456	3456		Example
3457	3457
3458	3458		A table of some values of the error function and the complementary error function is calculated in the following demo project:
3459	3459
3460	3460		``` pascal
3461	3461		program erf_Table;
3462	3462
3463	3463		uses
3464	3464		SysUtils, StrUtils, typ, spe;
3465	3465
3466	3466		const
3467	3467		Wx = 7;
3468	3468		Wy = 20;
3469	3469		D = 6;
3470	3470
3471	3471		var
3472	3472		i: Integer;
3473	3473		x: ArbFloat;
3474	3474		fs: TFormatSettings;
3475	3475		begin
3476	3476		fs := DefaultFormatSettings;
3477	3477		fs.DecimalSeparator := '.';
3478	3478		WriteLn('x':Wx, 'erf(x)':Wy+1, 'erfc(x)':Wy+1);
3479	3479		WriteLn;
3480	3480		for i:= 0 to 20 do begin
3481	3481		x := 0.1*i;
3482	3482		WriteLn(Format('%.1f %.f %.*f', [
3483	3483		Wx, x,
3484	3484		Wy, D, speerf(x),
3485	3485		Wy, D, speefc(x)
3486	3486		], fs));
3487	3487		end;
3488	3488		end.
3489	3489		```
3490	3490
3491	3491		### Normal distribution
3492	3492
3493	3493		![](/wiki-images/7/7b/Numlib_normaldist.png) The (cumulative) [normal distribution](https://en.wikipedia.org/wiki/Normal_distribution) is the integral of the Gaussian function with zero mean and unit standard deviation, i.e.
3494	3494
3495	3495		$`\operatorname{N}(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{x} \operatorname{exp}(\frac{t^2}{2}) dt
3496	3496		= \frac{1}{2}[1 + \operatorname{erf} (\frac{x}{\sqrt{2}})]`$
3497	3497
3498	3498		Suppose a set of random floating point numbers which are distributed according to this formula. Then N(x) measures the probability that a number from this set is smaller than x.
3499	3499
3500	3500		In NumLib, the normal distribution is calculated by the function `normaldist`:
3501	3501
3502	3502		``` pascal
3503	3503		function normaldist(x: ArbFloat): ArbFloat;
3504	3504		```
3505	3505
3506	3506		- `x` can take any value between -∞ and +∞.
3507	3507		- The function result is between 0 and 1.
3508	3508
3509	3509		The inverse normal distribution of y returns the x value for which the normal distribution N(x) has the value y.
3510	3510
3511	3511		``` pascal
3512	3512		function invnormaldist(y: ArbFloat): ArbFloat;
3513	3513		```
3514	3514
3515	3515		- `y` is allowed only between 0 and 1.
3516	3516		- The function can take any value between -∞ and +∞.
3517	3517
3518	3518		> Note: The functions `normaldist` and `invnormaldist` are not available in fpc 3.0.x and older.
3519	3519
3520	3520		Example
3521	3521
3522	3522		Calculate the probabilities for x = -10, -8, ..., 12 for a normal distribution with mean μ = 2 and standard deviation σ = 5. In order to reduce this distribution to the standard normal distribution implemented by NumLib, we introduce a transformation x<sub>r</sub> = (x - μ)/σ.
3523	3523
3524	3524		``` pascal
3525	3525		program normdist;
3526	3526
3527	3527		uses
3528	3528		typ, spe;
3529	3529
3530	3530		const
3531	3531		mu = 2.0;
3532	3532		sigma = 5.0;
3533	3533		x: array[0..12] of ArbFloat = (-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12, 14);
3534	3534
3535	3535		var
3536	3536		xr, y: ArbFloat;
3537	3537		i: Integer;
3538	3538		begin
3539	3539		WriteLn('Normal distribution with mean ', mu:0:1, ' and standard deviation ', sigma:0:1);
3540	3540		WriteLn;
3541	3541		WriteLn('x':20, 'y':20);
3542	3542
3543	3543		for i:= 0 to High(x) do begin
3544	3544		xr := (x[i] - mu) / sigma;
3545	3545		y := normaldist(xr);
3546	3546		WriteLn(x[i]:20:1, y:20:6);
3547	3547		end;
3548	3548		end.
3549	3549		```
3550	3550
3551	3551		### Gamma function
3552	3552
3553	3553		![](/wiki-images/8/8d/Numlib_gamma.png) The [gamma function](https://en.wikipedia.org/wiki/Gamma_function) is needed by many probability functions. It is defined by the integral
3554	3554
3555	3555
3556	3556		$`\Gamma({x}) = \int_0^{\infty}t^{x-1} e^{-t} dt`$
3557	3557
3558	3558		NumLib provides two functions for its calculation:
3559	3559
3560	3560		``` pascal
3561	3561		function spegam(x: ArbFloat): ArbFloat;
3562	3562		function spelga(x: ArbFloat): ArbFloat;
3563	3563		```
3564	3564
3565	3565		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.
3566	3566
3567	3567		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.
3568	3568
3569	3569		Example
3570	3570
3571	3571		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`.
3572	3572
3573	3573		``` pascal
3574	3574		program Gamma_Table;
3575	3575
3576	3576		uses
3577	3577		SysUtils, StrUtils, typ, spe;
3578	3578
3579	3579		const
3580	3580		VALUES: array[0..2] of ArbFloat = (1, 2, 5);
3581	3581		Wx = 7;
3582	3582		Wy = 30;
3583	3583		Wln = 20;
3584	3584
3585	3585		var
3586	3586		i: Integer;
3587	3587		x: ArbFloat;
3588	3588		magnitude: ArbFloat;
3589	3589		begin
3590	3590		WriteLn('x':Wx, 'Gamma(x)':Wy, 'ln(Gamma(x))':Wln);
3591	3591		WriteLn;
3592	3592		magnitude := 1E-3;
3593	3593		while magnitude <= 1000 do begin
3594	3594		for i := 0 to High(VALUES) do begin
3595	3595		x := VALUES[i] * magnitude;
3596	3596		Write(PadLeft(FloatToStr(x), Wx));
3597	3597		if x <= 170 then // Extended overflow above 170
3598	3598		Write(FormatFloat('0.000', spegam(x)):Wy)
3599	3599		else
3600	3600		Write('overflow':Wy);
3601	3601		WriteLn(spelga(x):Wln:3);
3602	3602		if abs(x-1000) < 1E-6 then
3603	3603		break;
3604	3604		end;
3605	3605		magnitude := magnitude * 10;
3606	3606		end;
3607	3607		end.
3608	3608		```
3609	3609
3610	3610		### Incomplete gamma function
3611	3611
3612	3612		![](/wiki-images/8/81/Numlib_gammapq.png) The lower and upper [incomplete gamma functions](https://en.wikipedia.org/wiki/Incomplete_Gamma_function) are defined by
3613	3613
3614	3614
3615	3615		$`\operatorname{P}({s},{x}) = \frac{1}{\Gamma({s})} \int_0^{x}t^{s-1} e^{-t} dt`$
3616	3616
3617	3617
3618	3618		$`\operatorname{Q}({s},{x}) = \frac{1}{\Gamma({s})} \int_{x}^{\infty}t^{s-1} e^{-t} dt = 1 - \operatorname{P}({s}, {x})`$
3619	3619
3620	3620		``` pascal
3621	3621		function gammap(s, x: ArbFloat): ArbFloat;
3622	3622		function gammaq(s, x: ArbFloat): ArbFloat;
3623	3623		```
3624	3624
3625	3625		- The parameter `s` must be positive: s \> 0
3626	3626		- The argument `x` must not be negative: x ≥ 0
3627	3627		- The function results are in the interval \[0, 1\].
3628	3628
3629	3629		> Note: The functions `gammap` and `gammaq` are not available in fpc v3.0.x and older.
3630	3630
3631	3631		Example
3632	3632
3633	3633		The following project prints out a table of `gammap` and `gammaq` values:
3634	3634
3635	3635		``` pascal
3636	3636		program IncompleteGamma_Table;
3637	3637
3638	3638		uses
3639	3639		SysUtils, StrUtils, typ, spe;
3640	3640
3641	3641		const
3642	3642		s = 1.25;
3643	3643
3644	3644		Wx = 7;
3645	3645		Wy = 20;
3646	3646		D = 6;
3647	3647
3648	3648		var
3649	3649		i: Integer;
3650	3650		x: ArbFloat;
3651	3651		fs: TFormatSettings;
3652	3652		begin
3653	3653		fs := DefaultFormatSettings;
3654	3654		fs.DecimalSeparator := '.';
3655	3655		WriteLn('x':Wx, 'GammaP('+FloatToStr(s,fs)+', x)':Wy+1, 'GammaQ('+FloatToStr(s,fs)+', x)':Wy+1);
3656	3656		WriteLn;
3657	3657		for i := 0 to 20 do begin
3658	3658		x := 0.5*i;
3659	3659		WriteLn(Format('%.1f %.f %.*f', [
3660	3660		Wx, x,
3661	3661		Wy, D, gammap(s, x),
3662	3662		Wy, D, gammaq(s, x)
3663	3663		], fs));
3664	3664		end;
3665	3665		end.
3666	3666		```
3667	3667
3668	3668		### Beta function
3669	3669
3670	3670		![](/wiki-images/3/32/Numlib_beta.png) The [Beta function](https://en.wikipedia.org/wiki/Beta_function) is another basis of important functions used in statistics. It is defined as
3671	3671
3672	3672
3673	3673		$`\operatorname{B}(a, b) = \frac{\Gamma(a+b)} = \int_0^1{t^{x-1} (1-t)^{y-1} dt}`$
3674	3674
3675	3675		The NumLib routine for its calculation is
3676	3676
3677	3677		``` pascal
3678	3678		function beta(a, b: ArbFloat): ArbFloat;
3679	3679		```
3680	3680
3681	3681		- The arguments `a` and `b` must be positive: `a`` > 0` and `b`` > 1`.
3682	3682		- The calculation is based on the logarithms of the gamma function and therefore does not suffer from the overflow issues of the gamma function.
3683	3683
3684	3684		> Note: The function `beta` is not available in fpc v3.0.x or older.
3685	3685
3686	3686		### Incomplete beta function
3687	3687
3688	3688		![](/wiki-images/3/3a/Numlib_betai.png) The (regularized) [incomplete beta function](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) is defined by
3689	3689
3690	3690
3691	3691		$`\operatorname{I}_x(a,b) = \frac {1}{\operatorname{B}(ab)} \int_0^x{t^{x-1} (1-t)^{y-1} dt}`$
3692	3692
3693	3693		``` pascal
3694	3694		function betai(a, b, x: ArbFloat): ArbFloat;
3695	3695		```
3696	3696
3697	3697		- The arguments `a` and `b` must be positive: `a`` > 0` and `b`` > 1`.
3698	3698		- `x` is only allowed to be in the range 0 ≤ x ≤ 1.
3699	3699		- The function result is in the interval \[0, 1\].
3700	3700
3701	3701		The inverse incomplete beta function determines the argument x for which the passed parameter is equal to `I`<sub>`x`</sub>`(a, b)`:
3702	3702
3703	3703		``` pascal
3704	3704		function invbetai(a, b, y: ArbFloat): ArbFloat;
3705	3705		```
3706	3706
3707	3707		- The arguments `a` and `b` must be positive: `a`` > 0` and `b`` > 1`.
3708	3708		- `y` is only allowed to be in the range 0 ≤ x ≤ 1.
3709	3709		- The function result is in the interval \[0, 1\].
3710	3710
3711	3711		> Note: The functions `betai` and `invbetai` are not available in fpc v3.0.x and older.
3712	3712
3713	3713		Example
3714	3714
3715	3715		Calculate the incomplete beta function for some values of x, a, b, and apply the result to the inverse incomplete beta function to demonstrate the validity of the calculations.
3716	3716
3717	3717		``` pascal
3718	3718		program incomplete_beta_func;
3719	3719
3720	3720		uses
3721	3721		SysUtils, typ, spe;
3722	3722
3723	3723		procedure TestBetai(a, b, x: ArbFloat);
3724	3724		var
3725	3725		xc, y: ArbFloat;
3726	3726		begin
3727	3727		y := betai(a, b, x);
3728	3728		xc := invbetai(a, b, y, macheps);
3729	3729		WriteLn(' x = ', x:0:6, ' ---> y = ', y:0:9, ' ---> x = ', xc:0:6);
3730	3730		end;
3731	3731
3732	3732		begin
3733	3733		WriteLn('Incomplete beta function and its inverse');
3734	3734		WriteLn;
3735	3735		WriteLn(' y = betai(1, 3, x) x = invbetai(1, 3, y)');
3736	3736		TestBetaI(1, 3, 0.001);
3737	3737		TestBetaI(1, 3, 0.01);
3738	3738		TestBetaI(1, 3, 0.1);
3739	3739		TestBetaI(1, 3, 0.2);
3740	3740		TestBetaI(1, 3, 0.9);
3741	3741		TestBetaI(1, 3, 0.99);
3742	3742		TestBetaI(1, 3, 0.999);
3743	3743		WriteLn;
3744	3744
3745	3745		WriteLn(' y = betai(1, 1, x) x = invbetai(1, 1, y)');
3746	3746		TestBetaI(1, 1, 0.001);
3747	3747		TestBetaI(1, 1, 0.01);
3748	3748		TestBetaI(1, 1, 0.1);
3749	3749		TestBetaI(1, 1, 0.2);
3750	3750		TestBetaI(1, 1, 0.9);
3751	3751		TestBetaI(1, 1, 0.99);
3752	3752		TestBetaI(1, 1, 0.999);
3753	3753		WriteLn;
3754	3754
3755	3755		WriteLn(' y = betai(3, 1, x) x = invbetai(3, 1, y)');
3756	3756		TestBetaI(3, 1, 0.001);
3757	3757		TestBetaI(3, 1, 0.01);
3758	3758		TestBetaI(3, 1, 0.1);
3759	3759		TestBetaI(3, 1, 0.2);
3760	3760		TestBetaI(3, 1, 0.9);
3761	3761		TestBetaI(3, 1, 0.99);
3762	3762		TestBetaI(3, 1, 0.999);
3763	3763		WriteLn;
3764	3764		end.
3765	3765		```
3766	3766
3767	3767		### Bessel functions
3768	3768
3769	3769		The [Bessel functions](https://en.wikipedia.org/wiki/Bessel_function) are solutions of the Bessel differential equation:
3770	3770
3771	3771
3772	3772		$`{x}^2 y'' + x y' + ({x}^2 - \alpha^2) {y} = 0`$
3773	3773
3774	3774		The Bessel functions of the first kind, or cylindrical harmonics, J<sub>α</sub>(x), are the solutions which are not singular at the origin, in contrast to the Bessel functions of the second kind, Y<sub>α</sub>(x). Y<sub>α</sub>(x) is only defined for x \> 0.
3775	3775
3776	3776		Example
3777	3777
3778	3778		- Calculate the probability that random normal-distributed numbers are within 1, 2, and 3 standard deviations around the mean.
3779	3779		- Calculate the maximum deviation of normal-distributed deviates (in units of standard deviations) from the mean at probabilities 0.90, 0.95, 0.98, 0.99 and 0.999.
3780	3780
3781	3781		![](/wiki-images/0/0c/Numlib_bessel_1st.png) ![](/wiki-images/d/de/Numlib_bessel_2nd.png)
3782	3782
3783	3783		The solutions for a purely imaginary argument are called modified Bessel functions of the first kind, I<sub>α</sub>(x), and modified Bessel function of the second kind, K<sub>α</sub>(x). K<sub>α</sub>(x) is only defined for x \> 0.
3784	3784
3785	3785		![](/wiki-images/b/b6/Numlib_mod_bessel_1st.png) ![](/wiki-images/a/ae/Numlib_mod_bessel_2nd.png)
3786	3786
3787	3787		NumLib implements only the solutions for the parameters α = 0 and α = 1:
3788	3788
3789	3789		``` pascal
3790	3790		function spebj0(x: ArbFloat): ArbFloat; // Bessel function of the first kind, J0 (alpha = 0)
3791	3791		function spebj1(x: ArbFloat): ArbFloat; // Bessel function of the first kind, J1 (alpha = 1)
3792	3792
3793	3793		function speby0(x: ArbFloat): ArbFloat; // Bessel function of the second kind, Y0 (alpha = 0)
3794	3794		function speby1(x: ArbFloat): ArbFloat; // Bessel function of the second kind, Y1 (alpha = 1)
3795	3795
3796	3796		function spebi0(x: ArbFloat): ArbFloat; // modified Bessel function of the first kind, I0 (alpha = 0)
3797	3797		function spebi1(x: ArbFloat): ArbFloat; // modified Bessel function of the first kind, I1 (alpha = 1)
3798	3798
3799	3799		function spebk0(x: ArbFloat): ArbFloat; // modified Bessel function of the second kind, K0 (alpha = 0)
3800	3800		function spebk1(x: ArbFloat): ArbFloat; // modified Bessel function of the second kind, K1 (alpha = 1)
3801	3801		```
3802	3802
3803	3803		### Other functions
3804	3804
3805	3805		The unit spe contains some other functions which will not be documented here because they are already available in th FPC's standard unit math.
3806	3806
3807		-	~~<thead>~~
3808		-
3809		-	~~in <em>math</em></th>~~
3810		-
3811		-	~~</thead>~~
3812		-	~~<tbody>~~
3813		-
3814		-	~~</tbody>~~
	3807	+	\| Function \| Equivalent function · in math \| Description \|
	3808	+	\|---\|---\|---\|
	3809	+	\| `function`` speent(x: ArbFloat): LongInt;` \| `floor(x)` \| Entier function, calculates first integer smaller than or equal to `x` \|
	3810	+	\| `function`` spemax(a, b: Arbfloat): ArbFloat;` \| `max(a, b)` \| Maximum of two floating point values \|
	3811	+	\| `function`` spepow(a, b: ArbFloat): ArbFloat;` \| `power(a, b)` \| Calculates ab \|
	3812	+	\| `function`` spesgn(x: ArbFloat): ArbInt;` \| `sign(x)` \| Returns the sign of `x` (`-1` for `x < 0`, `0` for `x = 0`, `+1` for `x > 0`) \|
	3813	+	\| `function`` spears(x: ArbFloat): ArbFloat;` \| `arcsin(x)` \| Inverse function of `sin(x)` \|
	3814	+	\| `function`` spearc(x: ArbFloat): ArbFloat;` \| `arccos(x)` \| Inverse function of `cos(x)` \|
	3815	+	\| `function`` spesih(x: ArbFloat): ArbFloat;` \| `sinh(x)` \| Hyperbolic sine \|
	3816	+	\| `function`` specoh(x: ArbFloat): ArbFloat;` \| `cosh(x)` \| Hyperbolic cosine \|
	3817	+	\| `function`` spetah(x: ArbFloat): ArbFloat;` \| `tanh(x)` \| Hyperbolic tangent \|
	3818	+	\| `function`` speash(x: ArbFloat): ArbFloat;` \| `arcsinh(x)` \| Inverse of the hyperbolic sine \|
	3819	+	\| `function`` speach(x: ArbFloat): ArbFloat;` \| `arccosh(x)` \| Inverse of the hyperbolic cosine \|
	3820	+	\| `function`` speath(x: ArbFloat): ArbFloat;` \| `arctanh(x)` \| Inverse of the hyperbolic tangent \|
3815	3821
3816	3822		---
3817	3823		From the Free Pascal wiki, page "NumLib Documentation", as it stood on 2024-02-14 after a change by Wp, out of the wiki's dump of 2024-05-26. Copyright stays with the author who wrote it, according to the old wiki's [[Licensing]] page. If you have concerns about this material, [message the admin](/wiki-contact).
3818	3824		Category: Mathematics
3819	3825		Category: Statistical algorithms

Back to history