1 {
2 This file is part of the Numlib package.
3 Copyright (c) 1986-2000 by
4 Kees van Ginneken, Wil Kortsmit and Loek van Reij of the
5 Computational centre of the Eindhoven University of Technology
6
7 FPC port Code by Marco van de Voort (marco@freepascal.org)
8 documentation by Michael van Canneyt (Michael@freepascal.org)
9
10 This is the most basic unit from NumLib.
11 The most important items this unit defines are matrix types and machine
12 constants
13
14 See the file COPYING.FPC, included in this distribution,
15 for details about the copyright.
16
17 This program is distributed in the hope that it will be useful,
18 but WITHOUT ANY WARRANTY; without even the implied warranty of
19 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
20
21 **********************************************************************}
22
23 {
24 In the FPC revision, instead of picking a certain floating point type,
25 a new type "ArbFloat" is defined, which is used as floating point type
26 throughout the entire library. If you change the floating point type, you
27 should only have to change ArbFloat, and the machineconstants belonging to
28 the type you want.
29 However for IEEE Double (64bit) and Extended(80bit) these constants are
30 already defined, and autoselected by the library. (the library tests the
31 size of the float type in bytes for 8 and 10 and picks the appropiate
32 constants
33
34 Also some stuff had to be added to get ipf running (vector object and
35 complex.inp and scale methods)
36 }
37
38 {$mode objfpc}{$H+}
39 {$modeswitch nestedprocvars}
40
41 unit typ;
42
43 {$I DIRECT.INC} {Contains "global" compilerswitches which
44 are imported into every unit of the library }
45
46
47 interface
48
49 uses
50 Math;
51
52 {$if sizeof(extended)=10}
53 {$DEFINE ArbExtended}
54 {$endif}
55
56
57 CONST numlib_version=2; {used to detect version conflicts between
58 header unit and dll}
59
60 type {Definition of base types}
61 {$IFDEF ArbExtended}
62 ArbFloat = extended;
63 {$ELSE}
64 ArbFloat = double;
65 {$ENDIF}
66 ArbInt = LONGINT;
67 ArbString = AnsiString;
68
69 Float8Arb =ARRAY[0..7] OF BYTE;
70 Float10Arb =ARRAY[0..9] OF BYTE;
71
72 CONST {Some constants for the variables below, in binary formats.}
73 {$IFNDEF ArbExtended}
74 {First for REAL/Double}
75 TC1 : Float8Arb = ($00,$00,$00,$00,$00,$00,$B0,$3C);
76 TC2 : Float8Arb = ($FF,$FF,$FF,$FF,$FF,$FF,$EF,$7F);
77 TC3 : Float8Arb = ($00,$00,$00,$00,$01,$00,$10,$00);
78 TC5 : Float8Arb = ($EF,$39,$FA,$FE,$42,$2E,$86,$40);
79 TC6 : Float8Arb = ($D6,$BC,$FA,$BC,$2B,$23,$86,$C0);
80 {$ENDIF}
81
82 {For Extended}
83 {$IFDEF ArbExtended}
84 TC1 : Float10Arb = (0,0,$00,$00,$00,$00,0,128,192,63); {Eps}
85 TC2 : Float10Arb = ($FF,$FF,$FF,$FF,$FF,$FF,$FF,$D6,$FE,127); {9.99188560553925115E+4931}
86 TC3 : Float10Arb = (1,0,0,0,0,0,0,0,0,0); {3.64519953188247460E-4951}
87 TC5 : Float10Arb = (18,25,219,91,61,101,113,177,12,64); {1.13563488668777920E+0004}
88 TC6 : Float10Arb = (108,115,3,170,182,56,27,178,12,192); {-1.13988053843083006E+0004}
89 {$ENDIF}
90 { numdig is the number of useful (safe) decimal places of an "ArbFloat"
91 for display.
92 minform is the number of decimal places shown by the rtls
93 write(x:ArbFloat)
94 maxform is the maximal number of decimal positions
95 }
96
97 numdig = 25;
98 minform = 10;
99 maxform = 26;
100
101 var
102 macheps : ArbFloat absolute TC1; { macheps = r - 1, with r
103 the smallest ArbFloat > 1}
104 giant : ArbFloat absolute TC2; { the largest ArbFloat}
105 midget : ArbFloat absolute TC3; { the smallest positive ArbFloat}
106 LnGiant : ArbFloat absolute TC5; {ln of giant}
107 LnMidget : ArbFloat absolute TC6; {ln of midget}
108
109 {Copied from Det. Needs ArbExtended conditional}
110 const { og = 8^-maxexp, og�>=midget,
111 bg = 8^maxexp, bg�<=giant
112
113 midget and giant are defined in typ.pas}
114
115 {$IFDEF ArbExtended}
116 ogx: Float10Arb = (51,158,223,249,51,243,4,181,224,31);
117 bgx: Float10Arb = (108,119,117,92,70,38,155,234,254,95);
118 maxexpx : ArbInt = 2740;
119 {$ELSE}
120 ogx: Float8Arb= (84, 254, 32, 128, 32, 0, 0, 32);
121 bgx: Float8Arb= (149, 255, 255, 255, 255, 255, 239, 95);
122 maxexpx : ArbInt = 170;
123 {$ENDIF}
124
125 var
126 og : ArbFloat absolute ogx;
127 bg : ArbFloat absolute bgx;
128 MaxExp : ArbInt absolute maxexpx;
129
130
131 {Like standard EXP(), but for very small values (near lowest possible
132 ArbFloat this version returns 0}
expnull133 Function exp(x: ArbFloat): ArbFloat;
134
135 type
136 Complex = object
137 { Crude complex record. For me an example of
138 useless OOP, specially if you have operator overloading
139 }
140 xreal, imag : ArbFloat;
141 procedure Init (r, i: ArbFloat);
142 procedure Add (c: complex);
143 procedure Sub (c: complex);
Inpnull144 function Inp(z:complex):ArbFloat;
145 procedure Conjugate;
146 procedure Scale(s: ArbFloat);
Normnull147 Function Norm : ArbFloat;
Sizenull148 Function Size : ArbFloat;
Renull149 Function Re : ArbFloat;
150 procedure Unary;
Imnull151 Function Im : ArbFloat;
Argnull152 Function Arg : ArbFloat;
153 procedure MinC(c: complex);
154 procedure MaxC(c: complex);
155 Procedure TransF(var t: complex);
156
157 end;
158
159 vector = object
160 i, j, k: ArbFloat;
161 procedure Init (vii, vjj, vkk: ArbFloat);
162 procedure Unary;
163 procedure Add (c: vector);
164 procedure Sub (c: vector);
Vinull165 function Vi : ArbFloat;
Vjnull166 function Vj : ArbFloat;
Vknull167 function Vk : ArbFloat;
Normnull168 function Norm : ArbFloat;
Norm8null169 Function Norm8 : ArbFloat;
Sizenull170 function Size : ArbFloat;
InProdnull171 function InProd(c: vector): ArbFloat;
172 procedure Uitprod(c: vector; var e: vector);
173 procedure Scale(s: ArbFloat);
174 procedure DScale(s: ArbFloat);
175 procedure Normalize;
176 procedure Rotate(calfa, salfa: ArbFloat; axe: vector);
177 procedure Show(p,q: ArbInt);
178 end;
179
180 transformorg = record offset: complex; ss, sc: real end;
181 transform = record
182 offsetx, offsety, scalex, scaley: ArbFloat
183 end;
184
185 {Standard Functions used in NumLib}
186 rfunc1r = Function(x : ArbFloat): ArbFloat;
187 rfunc1rn = Function(x : ArbFloat): ArbFloat is nested;
188 rfunc2r = Function(x, y : ArbFloat): ArbFloat;
189
190 {Complex version}
191 rfunc1z = Function(z: complex): ArbFloat;
192
193 {Special Functions}
194 oderk1n = procedure(x: ArbFloat; var y, f: ArbFloat);
195 roofnrfunc = procedure(var x, fx: ArbFloat; var deff: boolean);
196
197
198 {Maximal n x m dimensions of matrix.
199 +/- highestelement*SIZEOF(elementtype) is
200 minimal size of matrix.}
201 const
202 highestfloatelement = High(ArbInt) div SizeOf(ArbFloat);
203 highestptrelement = High(ArbInt) div SizeOf(Pointer);
204 highestintelement = High(ArbInt) div SizeOf(ArbInt);
205 highestboolelement = High(ArbInt) div SizeOf(boolean);
206 highestcomplexelement = High(ArbInt) div SizeOf(complex);
207 highestvectorelement = High(ArbInt) div SizeOf(vector);
208
209
210 type
211 {Definition of matrix types in NumLib. First some vectors.
212 The high boundery is a maximal number only. Vectors can be smaller, but
213 not bigger. The difference is the starting number}
214 arfloat0 = array[0..highestfloatelement-1] of ArbFloat;
215 arfloat1 = array[1..highestfloatelement] of ArbFloat;
216 arfloat2 = array[2..highestfloatelement+1] of ArbFloat;
217 arfloat_1 = array[-1..highestfloatelement-2] of ArbFloat;
218
219 {A matrix is an array of floats}
220 ar2dr = array[0..highestptrelement-1] of ^arfloat0;
221 ar2dr1 = array[1..highestptrelement] of ^arfloat1;
222
223 {Matrices can get big, so we mosttimes allocate them on the heap.}
224 par2dr1 = ^ar2dr1;
225
226 {Integer vectors}
227 arint0 = array[0..highestintelement-1] of ArbInt;
228 arint1 = array[1..highestintelement] of ArbInt;
229
230 {Boolean (true/false) vectors}
231 arbool1 = array[1..highestboolelement] of boolean;
232
233 {Complex vectors}
234 arcomp0 = array[0..highestcomplexelement-1] of complex;
235 arcomp1 = array[1..highestcomplexelement] of complex;
236 arvect0 = array[0..highestvectorelement-1] of vector;
237 vectors = array[1..highestvectorelement] of vector;
238
239 parcomp = ^arcomp1;
240
241 {(de) Allocate mxn matrix to A}
242 procedure AllocateAr2dr(m, n: integer; var a: par2dr1);
243 procedure DeAllocateAr2dr(m, n: integer; var a: par2dr1);
244
245 {(de) allocate below-left triangle matrix for (de)convolution
246 (a 3x3 matrix looks like this
247
248 x
249 x x
250 x x x)
251 }
252 procedure AllocateL2dr(n: integer; var a: par2dr1);
253 procedure DeAllocateL2dr(n: integer; var a: par2dr1);
254
255 {Get the Re and Im parts of a complex type}
Renull256 Function Re(z: complex): ArbFloat;
Imnull257 Function Im(z: complex): ArbFloat;
258
259 { Creates a string from a floatingpoint value}
R2Snull260 Function R2S(x: ArbFloat; p, q: integer): string;
261
262 {Calculate inproduct of V1 and V2, which are vectors with N elements;
263 I1 and I2 are the SIZEOF the datatypes of V1 and V2
264 MvdV: Change this to "V1,V2:array of ArbFloat and forget the i1 and i2
265 parameters?}
266
Inprodnull267 Function Inprod(var V1, V2; n, i1, i2: ArbInt): ArbFloat;
268
269 {Return certain special machine constants.(macheps=1, Nan=7)}
MachCnstnull270 Function MachCnst(n: ArbInt): ArbFloat;
271
dllversionnull272 function dllversion:LONGINT;
273
274 implementation
275
MachCnstnull276 Function MachCnst(n: ArbInt): ArbFloat;
277 begin
278 case n of
279 1: MachCnst := macheps;
280 2: MachCnst := giant;
281 3: MachCnst := midget;
282 4: MachCnst := infinity;
283 5: MachCnst := LnGiant;
284 6: MachCnst := LnMidget;
285 7: MachCnst := Nan;
286 end
287 end;
288
289 { Are used in many of the example programs}
Renull290 Function Re(z: complex): ArbFloat;
291 begin
292 Re := z.xreal
293 end;
294
Imnull295 Function Im(z: complex): ArbFloat;
296 begin
297 Im := z.imag
298 end;
299
300 {Kind of Sysutils.TrimRight and TrimLeft called after eachother}
301 procedure Compress(var s: string);
302 var i, j: LONGINT;
303 begin
304 j := length(s);
305 while (j>0) and (s[j]=' ') do dec(j);
306 i := 1;
307 while (i<=j) and (s[i]=' ') do Inc(i);
308 s := copy(s, i, j+1-i)
309 end;
310
R2Snull311 Function R2S(x: ArbFloat; p, q: integer): string;
312 var s: string;
313 i, j, k: integer;
314 begin
315 if q=-1 then
316 begin
317 Str(x:p, s);
318 i := Pos('E', s)-1; k := i+1;
319 j := i+3; while (j<length(s)) and (s[j]='0') do inc(j);
320 while s[i]='0' do dec(i); if s[i]='.' then dec(i);
321 if s[j]='0' then s := copy(s,1,i) else
322 if s[k]='-' then
323 s := copy(s, 1, i)+'E-'+Copy(s, j, length(s)+1-j)
324 else
325 s := copy(s, 1, i)+'E'+Copy(s, j, length(s)+1-j)
326 end
327 else
328 Str(x:p:q, s);
329 Compress(s);
330 R2S := s
331 end;
332
333 procedure AllocateAr2dr(m, n: integer; var a: par2dr1);
334 var i: integer;
335 begin
336 GetMem(a, m*SizeOf(pointer));
337 for i:=1 to m do GetMem(a^[i], n*SizeOf(ArbFloat))
338 end;
339
340 procedure DeAllocateAr2dr(m, n: integer; var a: par2dr1);
341 var i: integer;
342 begin
343 for i:=m downto 1 do FreeMem(a^[i], n*SizeOf(ArbFloat));
344 FreeMem(a, m*SizeOf(pointer));
345 a := Nil
346 end;
347
348 procedure AllocateL2dr(n: integer; var a: par2dr1);
349 var i: integer;
350 begin
351 GetMem(a, n*SizeOf(pointer));
352 for i:=1 to n do GetMem(a^[i], i*SizeOf(ArbFloat))
353 end;
354
355 procedure DeAllocateL2dr(n: integer; var a: par2dr1);
356 var i: integer;
357 begin
358 for i:=n downto 1 do FreeMem(a^[i], i*SizeOf(ArbFloat));
359 FreeMem(a, n*SizeOf(pointer));
360 a := Nil
361 end;
362
363 procedure Complex.Init(r, i: ArbFloat);
364 begin
365 xreal:= r;
366 imag := i
367 end;
368
369 procedure Complex.Conjugate;
370 begin
371 imag := -imag
372 end;
373
Complex.Inpnull374 function Complex.Inp(z:complex):ArbFloat;
375 begin
376 Inp := xreal*z.xreal + imag*z.imag
377 end;
378
379 procedure Complex.MinC(c: complex);
380 begin if c.xreal<xreal then xreal := c.xreal;
381 if c.imag<imag then imag := c.imag
382 end;
383
384 procedure Complex.Maxc(c: complex);
385 begin if c.xreal>xreal then xreal := c.xreal;
386 if c.imag>imag then imag := c.imag
387 end;
388
389 procedure Complex.Add(c: complex);
390 begin
391 xreal := xreal + c.xreal; imag := imag + c.imag
392 end;
393
394 procedure Complex.Sub(c: complex);
395 begin
396 xreal := xreal - c.xreal; imag := imag - c.imag
397 end;
398
Complex.Normnull399 Function Complex.Norm: ArbFloat;
400 begin
401 Norm := Sqr(xreal) + Sqr(imag)
402 end;
403
Complex.Sizenull404 Function Complex.Size: ArbFloat;
405 begin
406 Size := Sqrt(Norm)
407 end;
408
Complex.Renull409 Function Complex.Re: ArbFloat;
410 begin
411 Re := xreal;
412 end;
413
Complex.Imnull414 Function Complex.Im: ArbFloat;
415 begin
416 Im := imag
417 end;
418
419 Procedure Complex.TransF(var t: complex);
420 var w: complex;
421 tt: transformorg absolute t;
422 begin
423 w := Self; Conjugate;
424 with tt do
425 begin
426 w.scale(ss);
427 scale(sc);
428 Add(offset)
429 end;
430 Add(w)
431 end;
432
433
434 procedure Complex.Unary;
435 begin
436 xreal := -xreal;
437 imag := -imag
438 end;
439
440 procedure Complex.Scale(s:ArbFloat);
441 begin
442 xreal := xreal*s; imag := imag*s
443 end;
444
Complex.Argnull445 Function Complex.Arg: ArbFloat;
446 begin
447 if xreal=0 then
448 if imag>0 then Arg := 0.5*pi else
449 if imag=0 then Arg := 0 else Arg := -0.5*pi else
450 if xReal>0 then Arg := ArcTan(imag/xReal)
451 else if imag>=0 then Arg := ArcTan(imag/xReal) + pi
452 else Arg := ArcTan(imag/xReal) - pi
453 end;
454
expnull455 Function exp(x: ArbFloat): ArbFloat;
456 begin
457 if x<LnMidget then exp := 0 else exp := system.exp(x)
458 end;
459
460 { procedure berekent: v1 = v1 + r*v2 i1 en i2 geven de
461 increments in bytes voor v1 en v2 }
462
Inprodnull463 Function Inprod(var V1, V2; n, i1, i2: ArbInt): ArbFloat;
464
465 VAR i: LONGINT;
466 p1, p2: ^ArbFloat;
467 s: ArbFloat;
468 begin
469 IF I1 <>SIZEOF(ArbFloat) THEN
470 BEGIN
471 WRITELN('1 Something went probably wrong while porting!');
472 HALT;
473 END;
474 p1 := @v1; p2 := @v2; s := 0;
475 for i:=1 to n do
476 begin
477 s := s + p1^*p2^;
478 Inc(ptrint(p1), i1);
479 Inc(ptrint(p2), i2)
480 end;
481 Inprod := s
482 end;
483
484 procedure Vector.Init(vii, vjj, vkk: ArbFloat);
485 begin
486 i := vii; j := vjj; k := vkk
487 end;
488
489 procedure Vector.Unary;
490 begin i := -i; j := -j; k := -k end;
491
492 procedure Vector.Add(c: vector);
493 begin
494 i := i + c.i; j := j + c.j; k := k + c.k
495 end;
496
497 procedure Vector.Sub(c: vector);
498 begin
499 i := i - c.i; j := j - c.j; k := k - c.k
500 end;
501
Vector.Vinull502 function Vector.Vi : ArbFloat; begin Vi := i end;
503
Vector.Vjnull504 function Vector.Vj : ArbFloat; begin Vj := j end;
505
Vector.Vknull506 function Vector.Vk : ArbFloat; begin Vk := k end;
507
Vector.Normnull508 function Vector.Norm:ArbFloat;
509 begin
510 Norm := Sqr(i) + Sqr(j) + Sqr(k)
511 end;
512
Vector.Norm8null513 function Vector.Norm8:ArbFloat;
514 var r: ArbFloat;
515 begin
516 r := abs(i);
517 if abs(j)>r then r := abs(j);
518 if abs(k)>r then r := abs(k);
519 Norm8 := r
520 end;
521
Vector.Sizenull522 function Vector.Size: ArbFloat;
523 begin
524 Size := Sqrt(Norm)
525 end;
526
Vector.InProdnull527 function Vector.InProd(c: vector): ArbFloat;
528 begin
529 InProd := i*c.i + j*c.j + k*c.k
530 end;
531
532 procedure Vector.Uitprod(c: vector; var e: vector);
533 begin
534 e.i := j*c.k - k*c.j;
535 e.j := k*c.i - i*c.k;
536 e.k := i*c.j - j*c.i
537 end;
538
539 procedure Vector.Scale(s: ArbFloat);
540 begin
541 i := i*s; j := j*s; k := k*s
542 end;
543
544 procedure Vector.DScale(s: ArbFloat);
545 begin
546 i := i/s; j := j/s; k := k/s
547 end;
548
549 procedure Vector.Normalize;
550 begin
551 DScale(Size)
552 end;
553
554 procedure Vector.Show(p,q:ArbInt);
555 begin writeln(i:p:q, 'I', j:p:q, 'J', k:p:q, 'K') end;
556
557 procedure Vector.Rotate(calfa, salfa: arbfloat; axe: vector);
558 var qv : vector;
559 begin
560 Uitprod(axe, qv); qv.scale(salfa);
561 axe.scale((1-calfa)*Inprod(axe));
562 scale(calfa); sub(qv); add(axe)
563 end;
564
dllversionnull565 function dllversion:LONGINT;
566
567 BEGIN
568 dllversion:=numlib_version;
569 END;
570
571
572 END.
573