home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / turbopas / tpmath.arc / MATH.PAS next >
Pascal/Delphi Source File  |  1989-12-12  |  41KB  |  1,341 lines
  1. Unit Math;
  2. Interface
  3.  
  4. (*
  5.  ╔════════════════════════════════════════════════════════════════════════════╗
  6.  ║ Standard Mathematical Functions.        Version 1.0                        ║
  7.  ║                                                                            ║
  8.  ║ This code has been placed in the public domain.                            ║
  9.  ║ All non-trivial functions that do not reside permanently inside my head    ║
  10.  ║ (e.g. Bessel Functions) were obtain from:                                  ║
  11.  ║                                                                            ║
  12.  ║   _CRC_Standard_Mathematical_Tables_, CRC Press, 25th Edition, William H.  ║
  13.  ║   Beyer, editor.                                                           ║
  14.  ║                                                                            ║
  15.  ║   HP-15C Owner's Manual, Feb. 1984, Hewlett-Packard Company.               ║
  16.  ║                                                                            ║
  17.  ║   _Probability,_Random_Variables,_and_Stochastic_Processes_, Athanasios    ║
  18.  ║   Papoulis, Second Edition, 1984, McGraw-Hill                              ║
  19.  ║                                                                            ║
  20.  ║                                                                            ║
  21.  ║ All code has been tested, though not rigorusly, and is correct to the best ║
  22.  ║ of my knowledge. However, I make no guarentees. If you find any errors or  ║
  23.  ║ make any changes/enhancements, I would appreciate it if you inform me. My  ║
  24.  ║ address is:                                                                ║
  25.  ║                                                                            ║
  26.  ║ Scott Bostater                                                             ║
  27.  ║ Georgia Tech Research Institute                                            ║
  28.  ║ 7220 Richardson Road                                                       ║
  29.  ║ Symrna, GA 30080                                                           ║
  30.  ║ (404)-528-7782    (please, only in an emergency, my boss to work for HIM)  ║
  31.  ║                                                                            ║
  32.  ║ e-mail:                                                                    ║
  33.  ║ ---------------------------------------                                    ║
  34.  ║ Internet: bb16@prism.gatech.edu                                            ║
  35.  ║ UUCP:     ...!{allegra,amd,hplabs,ut-ngp}!gatech!edu!bb16                  ║
  36.  ║ Bitnet:   BBOSTATE@GTRI01                                                  ║
  37.  ║                                                                            ║
  38.  ║                                                                            ║
  39.  ║ Special Thanks to Mike Baden and Stan West for thier contributions         ║
  40.  ║                                                                            ║
  41.  ╚════════════════════════════════════════════════════════════════════════════╝
  42. *)
  43.  
  44. { Delcare Standard Math types}
  45. Type
  46.   {$IFOPT N+}
  47.     Float = Double;
  48.   {$Else}
  49.     Float = Real;
  50.   {$ENDIF}
  51.   Complex = Record
  52.               Real: Float;
  53.               Imag: Float;
  54.             End;
  55.  
  56. {*
  57.  * Run of the mill math functions that I always need
  58.  *}
  59.  
  60. Function  Log                 (x:Float):Float;                  { Log base 10 }
  61. Function  Power               (x,y:Float):Float;
  62. Function  Atan360             (x,y:Float):Float;
  63. Function  Tan                 (x:Float):Float;
  64. Function  ArcSin              (x:Float):Float;
  65. Function  ArcCos              (x:Float):Float;
  66. Function  HypSin              (x:Float):Float;
  67. Function  HypCos              (x:Float):Float;
  68.  
  69. Function  Inverse_tangent_deg (X,Y:Float):Float;
  70. Function  Inverse_tangent_rad (X,Y:Float):Float;
  71. Function  Inverse_sin_deg     (Y,R:Float):Float;
  72. Function  Inverse_sin_rad     (Y,R:Float):Float;
  73. Function  Inverse_cos_deg     (X,R:Float):Float;
  74. Function  Inverse_cos_rad     (X,R:Float):Float;
  75.  
  76. {*
  77.  * Complex math
  78.  *}
  79.  
  80. Procedure Complex_add         (A,B:Complex; Var C:Complex);
  81. Procedure Complex_subtract    (A,B:Complex; Var C:Complex);
  82. Procedure Complex_multiply    (A,B:Complex; Var C:Complex);
  83. Procedure Complex_divide      (A,B:Complex; Var C:Complex);
  84. Function  Magnitude           (A:Complex):Float;
  85. Procedure complex_conjugate   (A:Complex; var c:complex);
  86.  
  87.  
  88. {*
  89.  * Bessel Functions, (I hate Bessel Functions!)
  90.  *}
  91.  
  92. function  I0                  (x:Float):Float;   { modified Bessel of Order 0 }
  93. function  I1                  (x:Float):Float;   { modified Bessel of Order 1 }
  94. function  I2                  (x:Float):Float;   { modified Bessel of Order 2 }
  95. function  J0                  (x:Float):Float;   { Bessel function of Order 0 }
  96. function  J1                  (x:Float):Float;   { Bessel function of Order 1 }
  97.  
  98. {*
  99.  * Family of Unit Step functions
  100.  *
  101.  * U(0,t)  = 1 (t=0), = 0 (else)
  102.  * U(-1,t) = 1 (t>0), 1/2 (t=0), 0 (t<0)     [integral of U(0,t)]
  103.  * U(-2,t) = t (t>=0), 0 (t<0)               [integral of U(-1,t)]
  104.  * U(-3,t) = ½t² (t>=0), 0 (t<0)             [integral of U(-2,t)]
  105.  *}
  106.  
  107. Function  U                   (index:Integer; T:Float):Float;
  108.  
  109. {*
  110.  * Voltage and Power Decibel Functions
  111.  *}
  112.  
  113. Function  db10                (x:Float):Float;
  114. Function  db20                (x:Float):Float;
  115. Function  undb10              (X:Float):Float;
  116. Function  undb20              (x:Float):Float;
  117.  
  118. {*
  119.  * Hex output Functions  (TP really should have had these built in)
  120.  *}
  121.  
  122. type
  123.   String9 = String[9];
  124.   String4 = String[4];
  125.   String2 = String[2];
  126.  
  127. Function  HexLong             (a:LongInt):String9;
  128. Function  HexWord             (a:Word):   String4;
  129. Function  HexInt              (a:Integer):String4;
  130. Function  HexByte             (a:Byte):   String2;
  131.  
  132. {*
  133.  * If you have an 80387, These commands are _much_ better than the standard
  134.  * Sin/Cos functions that TP provides
  135.  *}
  136.  
  137. {$Ifdef P387 }
  138. Function Cos( x: Float): Float;
  139. Function Sin( x: Float): Float;
  140. {$EndIf}
  141.  
  142. {*
  143.  * Statistical Routines
  144.  *}
  145. type
  146.   StatsArray = Array[1..2048] of Float;
  147.   StatRecord = Record
  148.                  Nx:    0..2048;                    { Number of x points      }
  149.                  Ny:    0..2048;                    { Number of y points      }
  150.                  x:     ^StatsArray;                { x data                  }
  151.                  y:     ^StatsArray;                { y data                  }
  152.                  MeanX: Float;                      { X average               }
  153.                  MeanY: Float;                      { Y average               }
  154.                  StdX:  Float;                      { Standard Deviation of x }
  155.                  StdY:  Float;                      { Standard Deviation of y }
  156.                  r:     Float;                      { correlation coefficient }
  157.                End; { Record }
  158.  
  159. Procedure GetStatMemory(  var a: StatRecord);
  160. Procedure FreeStatMemory( var a: StatRecord);
  161. Procedure ComputeStats(   var a: StatRecord);       { mean, std dev, cor.coef }
  162. Function  Gauss( mean, StdDev: Float):Float;        { sample from Gauss dist.}
  163.  
  164. {*
  165.  * Matirx Manipulations
  166.  *}
  167.  
  168. Const
  169.   Arraysize = 32;
  170.   Arraysize2= 64;
  171.  
  172. Type
  173.   Vector   = array[0..ArraySize]  of Real;
  174.   Vector2  = array[0..ArraySize2] of real;
  175.   Matrix   = array[0..ArraySize] of Vector;
  176.   Matrix2  = array[0..ArraySize2] of Vector2;
  177.   Matrix2ptr = ^Matrix2;
  178.  
  179. Procedure MatrixMultiply         (    Dimen:             Integer;
  180.                                   Var InputMatrix1:      Matrix;
  181.                                   Var InputMatrix2:      Matrix;
  182.                                   Var ProductMatrix:     Matrix);
  183. Procedure MatrixAdd              (    Dimen:             Integer;
  184.                                   Var InputMatrix1:      Matrix;
  185.                                   Var InputMatrix2:      Matrix;
  186.                                   Var AddedMatrix:       Matrix);
  187. Procedure MatrixSubtract         (    Dimen:             Integer;
  188.                                   Var InputMatrix1:      Matrix;
  189.                                   Var InputMatrix2:      Matrix;
  190.                                   Var Matrix1Minus2:     Matrix);
  191. Procedure MatrixScalarMultiply   (    Dimen:             Integer;
  192.                                   Var InputMatrix:       Matrix;
  193.                                       Scalar:            Real;
  194.                                   Var ProductMatrix:     Matrix);
  195. Procedure BigMatrixScalarMultiply(    Dimen:             Integer;
  196.                                   Var InputMatrix:       Matrix2;
  197.                                       Scalar:            Real;
  198.                                   Var ProductMatrix:     Matrix2);
  199. Procedure BigPtrMatrixScalarMultiply(    Dimen:             Integer;
  200.                                   Var InputMatrix:       Matrix2ptr;
  201.                                       Scalar:            Real;
  202.                                   Var ProductMatrix:     Matrix2ptr);
  203. Procedure ComplexMatrixInverse   (    Dimen:             Integer;
  204.                                   Var RealInputMatrix:   Matrix;
  205.                                   Var ImagInputMatrix:   Matrix;
  206.                                   Var RealInverseMatrix: Matrix;
  207.                                   Var ImagInverseMatrix: Matrix;
  208.                                   Var Error:             Byte);
  209. Procedure MatrixInverse          (    Dimen:             Integer;
  210.                                   Var Input:             Matrix;
  211.                                   Var Output:            Matrix;
  212.                                   Var Error:             Byte);
  213. Procedure BigMatrixInverse       (    Dimen:             Integer;
  214.                                   Var Input:             Matrix2;
  215.                                   Var Output:            Matrix2;
  216.                                   Var Error:             Byte);
  217. Procedure BigPtrMatrixInverse    (    Dimen:             Integer;
  218.                                   Var Input:             Matrix2Ptr;
  219.                                   Var Output:            Matrix2Ptr;
  220.                                   Var Error:             Byte);
  221.  
  222.  
  223. Implementation
  224. const
  225.   Minus2: Integer = -2;
  226.  
  227.  
  228.  
  229. {$IFDef P387 }                                    { If 80387 available, load }
  230. Function Cos( x: Float): Float; external;         { in 80387 specific        }
  231. Function Sin( x: Float): Float; external;         { functions}
  232. Function Tan( x: Float): Float; external;
  233. {$L trig387.obj }
  234.  
  235. {$ELSE}
  236. {$IFOPT N+}
  237. Function Tan( X: Float): Float; external;         { 8087 or 80287 Presend     }
  238.                                                   { Load in Assembly routine  }
  239. {$L trig.obj }                                    { for fast Tan(x)           }
  240.  
  241. {$Else}
  242. Function Tan( X: Float): Float;                   { No Coprocessor present.   }
  243. Begin                                             { Out of luck. Do it the old}
  244.   Tan := Sin(x)/Cos(x);                           { fashioned way             }
  245. End;
  246. {$ENDIF}
  247. {$EndIf}
  248.  
  249.  
  250. {═════════════════════════════════════════════════════════════════════════════}
  251.  
  252. Function U( index: Integer; T: Float): Float;
  253. var
  254.   P: ShortInt;
  255.  
  256. Begin
  257.   P := -1-index;
  258.   U := 0;              { default }
  259.  
  260.   If T = 0 Then
  261.    Begin
  262.     If index = 0 Then U := 1
  263.    End
  264.   Else If T > 0 Then
  265.     Case P of
  266.      0: U := 1;
  267.      1: U := T;
  268.      2: U := Sqr(t)/2;
  269.      3..127: U := Power(t,P)/P;
  270.     End; { Case }
  271. End; { Function U }
  272.  
  273. {═════════════════════════════════════════════════════════════════════════════}
  274.  
  275. Function Magnitude (A:Complex):Float;
  276. Begin
  277.   Magnitude := Sqrt(Sqr(A.Real) + Sqr(A.Imag));
  278. End;
  279.  
  280. {═════════════════════════════════════════════════════════════════════════════}
  281.  
  282. Procedure complex_conjugate  (A:Complex; var c:complex);
  283. Begin
  284.   C.Real := A.Real;
  285.   C.Imag :=-A.Imag;
  286. End;
  287.  
  288. {═════════════════════════════════════════════════════════════════════════════}
  289.  
  290. Function Log (X:Float):Float;
  291. Const ln10 = 2.302585093;
  292. Begin
  293.   Log := ln(X) / ln10;
  294. End;
  295.  
  296. {═════════════════════════════════════════════════════════════════════════════}
  297.  
  298. Procedure Complex_Add (A,B:Complex; Var C:Complex);
  299. { C = A + B }
  300. Begin
  301.   C.Real := A.Real + B.Real;
  302.   C.Imag := A.Imag + B.Imag;
  303. End;
  304.  
  305. {═════════════════════════════════════════════════════════════════════════════}
  306.  
  307. Procedure Complex_subtract (A,B:Complex; Var C:Complex);
  308. { C = A - B }
  309. Begin
  310.   C.Real := A.Real - B.Real;
  311.   C.Imag := A.Imag - B.Imag;
  312. End;
  313.  
  314. {═════════════════════════════════════════════════════════════════════════════}
  315.  
  316. Procedure Complex_multiply (A,B:Complex; Var C:Complex);
  317. { C = A * B }
  318. Begin
  319.   C.Real := A.Real * B.Real  -  A.Imag * B.Imag;
  320.   C.Imag := A.Real * B.Imag  +  A.Imag * B.Real;
  321. End;
  322.  
  323. {═════════════════════════════════════════════════════════════════════════════}
  324.  
  325. Procedure Complex_divide (A,B:Complex; Var C:Complex);
  326. { C = A / B }
  327. Var
  328.   Temp:Real;
  329. Begin
  330.   Temp   :=  B.Real * B.Real  +  B.Imag * B.Imag;
  331.   C.Real := (A.Real * B.Real  +  A.Imag * B.Imag) / Temp;
  332.   C.Imag := (A.Imag * B.Real  -  A.Real * B.Imag) / Temp;
  333. End;
  334.  
  335. {═════════════════════════════════════════════════════════════════════════════}
  336.  
  337. Function Inverse_tangent_deg(X,Y:Float):Float;
  338. { 0 <= result <= 360 }
  339. Var
  340.   Angle:Float;
  341. Begin
  342.   if X=0 then Angle := Pi/2.0 else angle := ArcTan(Y/X);
  343.   angle := angle*180.0/pi;
  344.   if (X<=0.0) and (Y<0) then angle := angle - 180.0;
  345.   if (X< 0.0) and (Y>0) then angle := angle + 180.0;
  346.   If angle < 0 then angle := angle+360.0;
  347.   Inverse_tangent_deg := angle;
  348. End;
  349.  
  350. {═════════════════════════════════════════════════════════════════════════════}
  351.  
  352. Function Inverse_tangent_rad(X,Y:Float):Float;
  353. { 0 <= result <= 2pi }
  354. Var
  355.   Angle:Float;
  356. Begin
  357.   if X=0 then Angle := Pi/2.0 else angle := ArcTan(Y/X);
  358.   if (X<=0.0) and (Y<0) then angle := angle - pi;
  359.   if (X< 0.0) and (Y>0) then angle := angle + pi;
  360.   If angle < 0 then angle := angle+2*pi;
  361.   Inverse_tangent_rad := angle;
  362. End;
  363.  
  364. {═════════════════════════════════════════════════════════════════════════════}
  365.  
  366. Function Inverse_sin_deg (Y,R:Float):Float;
  367. { -90 <= result <= 90 }
  368. Var
  369.   X:Float;
  370.   temp:Float;
  371. Begin
  372.   X := Sqrt(Sqr(R)-Sqr(Y));
  373.   Temp := Inverse_tangent_deg(X,Y);
  374.   If Temp > 90 then Temp := Temp - 360;
  375.   Inverse_sin_deg := Temp;
  376. End;
  377.  
  378. {═════════════════════════════════════════════════════════════════════════════}
  379.  
  380. Function Inverse_sin_rad (Y,R:Float):Float;
  381. { -90 <= result <= 90 }
  382. Var
  383.   X:Float;
  384.   temp:Float;
  385. Begin
  386.   X := Sqrt(Sqr(R)-Sqr(Y));
  387.   Temp := Inverse_tangent_rad(X,Y);
  388.   If Temp > 90 then Temp := Temp - 360;
  389.   Inverse_sin_rad := Temp;
  390. End;
  391.  
  392. {═════════════════════════════════════════════════════════════════════════════}
  393.  
  394. Function Inverse_cos_deg (X,R:Float):Float;
  395. { -90 <= result <= 90 }
  396. Var
  397.   Y:Float;
  398.   temp:Float;
  399. Begin
  400.   Y := Sqrt(Sqr(R)-Sqr(X));
  401.   Temp := Inverse_tangent_deg(X,Y);
  402.   If Temp > 90 then Temp := Temp - 360;
  403.   Inverse_cos_deg := Temp;
  404. End;
  405.  
  406. {═════════════════════════════════════════════════════════════════════════════}
  407.  
  408. Function Inverse_cos_rad (X,R:Float):Float;
  409. { -90 <= result <= 90 }
  410. Var
  411.   Y:Float;
  412.   temp:Float;
  413. Begin
  414.   Y := Sqrt(Sqr(R)-Sqr(X));
  415.   Temp := Inverse_tangent_rad(X,Y);
  416.   If Temp > 90 then Temp := Temp - 360;
  417.   Inverse_cos_rad := Temp;
  418. End;
  419.  
  420. {═════════════════════════════════════════════════════════════════════════════}
  421.  
  422. Function ArcSin( x: Float): Float;
  423. Begin
  424.   ArcSin := Arctan( x / Sqrt( 1-Sqr(x)));
  425. End;
  426.  
  427. {═════════════════════════════════════════════════════════════════════════════}
  428.  
  429. const
  430.   pi4 = 0.785398164;
  431.  
  432. Function ArcCos( x: Float): Float;
  433. Begin
  434.   If x > 1 Then
  435.     ArcCos := 0
  436.   Else if x < -1 Then
  437.     ArcCos := pi
  438.   Else
  439.     ArcCos := pi4 - Arctan( x / Sqrt( 1-Sqr(x)));
  440. End;
  441.  
  442. {═════════════════════════════════════════════════════════════════════════════}
  443.  
  444. Function HypSin( x: Float): Float;
  445. Begin
  446.   Hypsin := 0.5 * (exp(x) - exp(-x));
  447. End;
  448.  
  449. {═════════════════════════════════════════════════════════════════════════════}
  450.  
  451. Function HypCos( x: Float): Float;
  452. Begin
  453.   HypCos := 0.5 * (exp(x) + exp(-x));
  454. End;
  455.  
  456. {═════════════════════════════════════════════════════════════════════════════}
  457.  
  458. function I0(x: Float): Float;
  459. var
  460.   tnum,factor: word;
  461.   term,i0temp: Float;
  462. begin
  463.   term:=1;
  464.   i0temp:=term;
  465.   factor:=0;
  466.   tnum:=1;
  467.   repeat
  468.     factor:=factor+2+2*(tnum-1);
  469.     term:=term*sqr(x/2)/factor;
  470.     i0temp:=i0temp+term;
  471.     inc(tnum);
  472.   until (abs(term)<1e-15) or (tnum=0);      {no overflow}
  473.   i0:=i0temp;
  474. end; {I0}
  475.  
  476. {═════════════════════════════════════════════════════════════════════════════}
  477.  
  478. function I1(x: Float): Float;
  479. var
  480.   tnum,factor: word;
  481.   term,i1temp: Float;
  482. begin
  483.   term:=x/4;
  484.   i1temp:=term;
  485.   factor:=0;
  486.   tnum:=1;
  487.   repeat
  488.     factor:=factor+3+2*(tnum-1);
  489.     term:=term*sqr(x/2)/factor;
  490.     writeln(term);
  491.     i1temp:=i1temp+term;
  492.     inc(tnum);
  493.   until (abs(term)<1e-15) or (tnum=0);      {no overflow}
  494.   i1:=i1temp;
  495. end; {I0}
  496.  
  497. {═════════════════════════════════════════════════════════════════════════════}
  498.  
  499. function I2(x: Float): Float;
  500. var
  501.   tnum,factor: word;
  502.   term,i2temp: Float;
  503. begin
  504.   term:=sqr(x)/24;
  505.   i2temp:=term;
  506.   factor:=0;
  507.   tnum:=1;
  508.   repeat
  509.     factor:=factor+4+2*(tnum-1);
  510.     term:=term*sqr(x/2)/factor;
  511.     writeln(term);
  512.     i2temp:=i2temp+term;
  513.     inc(tnum);
  514.   until (abs(term)<1e-15) or (tnum=0);      {no overflow}
  515.   i2:=i2temp;
  516. end; {I0}
  517.  
  518. {═════════════════════════════════════════════════════════════════════════════}
  519.  
  520. function J0(x: Float): Float;
  521. var
  522.   tnum: word;
  523.   term,j0temp: Float;
  524. begin
  525.   j0temp:=1;
  526.   term:=1;
  527.   tnum:=1;
  528.   repeat
  529.     term:=term*-sqr(x)/(4*sqr(tnum));
  530.     j0temp:=j0temp+term;
  531.     inc(tnum);
  532.   until (abs(term)<1e-15) or (tnum=0);      {no overflow}
  533.   j0:=j0temp;
  534. end; {J0}
  535.  
  536. {═════════════════════════════════════════════════════════════════════════════}
  537.  
  538. function J1(x: Float): Float;
  539. var
  540.   tnum,factor: word;
  541.   term,j1temp: Float;
  542. begin
  543.   term:=x/2;
  544.   j1temp:=term;
  545.   factor:=0;
  546.   tnum:=1;
  547.   repeat
  548.     factor:=factor+2+2*(tnum-1);
  549.     term:=term*-sqr(x)/(4*factor);
  550.     j1temp:=j1temp+term;
  551.     inc(tnum);
  552.   until (abs(term)<1e-15) or (tnum=0);      {no overflow}
  553.   j1:=j1temp;
  554. end; {J1}
  555.  
  556. {═════════════════════════════════════════════════════════════════════════════}
  557.  
  558. {═════════════════════════════════════════════════════════════════════════════}
  559.  
  560. Function Power( x,y: Float): Float;
  561. Begin
  562.   If y = 0 Then
  563.     power := 1.0
  564.   Else if x = 0 Then
  565.     Power := 0.0
  566.   Else If x > 0 Then
  567.     Power := exp( y * ln(x))
  568.   Else if Trunc(y) mod 2 = 0 Then
  569.     Power := exp( y * ln(abs(x)))
  570.   Else
  571.     Power := -exp( y * ln(abs(x)));
  572. End;
  573.  
  574. {═════════════════════════════════════════════════════════════════════════════}
  575.  
  576. Function Atan360( x, y: Float): Float;
  577. Var
  578.   Angle: Float;
  579. Begin
  580.   if X=0 then Angle := Pi/2.0 else angle := ArcTan(Y/X);
  581.   angle := angle*180.0/pi;
  582.   if (X<=0.0) and (Y<0) then angle := angle - 180.0;
  583.   if (X< 0.0) and (Y>0) then angle := angle + 180.0;
  584.   If angle < 0 then angle := angle+360.0;
  585.   Atan360 := angle;
  586. End;
  587.  
  588. {═════════════════════════════════════════════════════════════════════════════}
  589.  
  590. {$IFDEF N+}       { Limits change depending on whter 6 byte or 8 byte real }
  591. const
  592.   ln10: Float  = 2.30258509299405E+0000;
  593.  
  594. Function db10( X: Float): Float;
  595. Begin
  596.   If x < 1.0e-50 Then
  597.     db10 := -500
  598.   Else
  599.     db10 := 10 * ln(x) / ln10;
  600. End; { Function db10 }
  601.  
  602. {═════════════════════════════════════════════════════════════════════════════}
  603.  
  604. Function db20( X: Float): Float;
  605. Begin
  606.   If x < 1.0e-50 Then
  607.     db20 := -1000
  608.   Else
  609.     db20 := 20 * ln(x) / ln10;
  610. End; { Function db20 }
  611.  
  612.  
  613. {═════════════════════════════════════════════════════════════════════════════}
  614. {$ELSE}         { using 6 byte reals !}
  615. const
  616.   ln10: Float = 2.30258509299405E+0000;
  617.  
  618. Function db10( X: Float): Float;
  619. Begin
  620.   If x < 1.0e-25 Then
  621.     db10 := -250
  622.   Else
  623.     db10 := 10 * ln(x) / ln10;
  624. End; { Function db10 }
  625.  
  626. {═════════════════════════════════════════════════════════════════════════════}
  627.  
  628. Function db20( X: Float): Float;
  629. Begin
  630.   If x < 1.0e-25 Then
  631.     db20 := -500
  632.   Else
  633.     db20 := 20 * ln(x) / ln10;
  634. End; { Function db20 }
  635.  
  636. {$ENDIF}
  637. {═════════════════════════════════════════════════════════════════════════════}
  638.  
  639. Function Undb10( X: Float): Float;
  640. Begin
  641.   Undb10 := Power(10, X/10);
  642. End; { Function db10 }
  643.  
  644. {═════════════════════════════════════════════════════════════════════════════}
  645.  
  646. Function Undb20( X: Float): Float;
  647. Begin
  648.   Undb20 := Power(10, X/20);
  649. End; { Function db20 }
  650.  
  651. {═════════════════════════════════════════════════════════════════════════════}
  652.  
  653. const
  654.   h: Array[0..15] of Char = ( '0', '1', '2', '3', '4', '5', '6', '7',
  655.                               '8', '9', 'A', 'B', 'C', 'D', 'E', 'F');
  656.  
  657. Function HexLong( a: LongInt): String9;
  658. Begin
  659.   HexLong := H[a shr 28] + H[(a shr 24) and $f] + H[ (a shr 20) and $f] +
  660.              H[ (a shr 16) and $f] + ' ' + H[(a shr 12) and $f] +
  661.              H[ (a shr 8) and $f] + H[ (a shr 4) and $f] + H[ a and $f];
  662. End;
  663.  
  664. {═════════════════════════════════════════════════════════════════════════════}
  665.  
  666. Function HexWord( a: Word): String4;
  667. Begin
  668.   HexWord := H[a shr 12] + H[ (a shr 8) and $f] + H[ (a shr 4) and $f] +
  669.              H[ a and $f];
  670. End;
  671.  
  672. {═════════════════════════════════════════════════════════════════════════════}
  673.  
  674. Function HexInt( a: Integer): String4;
  675. Begin
  676.   HexInt := H[a shr 12] + H[ (a shr 8) and $f] + H[ (a shr 4) and $f] +
  677.             H[ a and $f];
  678. End;
  679.  
  680. {═════════════════════════════════════════════════════════════════════════════}
  681.  
  682. Function HexByte( a: Byte): String2;
  683. Begin
  684.   HexByte := H[ a shr 4 and $f] + H[ a and $f];
  685. End;
  686.  
  687. {═════════════════════════════════════════════════════════════════════════════}
  688.  
  689. Procedure GetStatMemory( var a: StatRecord );
  690.  
  691. { get free memory for a.x^ and a.y^ for the number of points specified in
  692.   a.Nx and a.Ny.  It also clears out the previous contents of meanX...CovXY  }
  693.  
  694. Begin
  695.  
  696.   FillChar( a.meanX, Sizeof(Float)*5, 0);             { clear Stats data }
  697.  
  698.   If a.Nx > 0 Then
  699.    Begin
  700.     GetMem( a.x, Sizeof(Float)*a.Nx);
  701.     FillChar( a.x^[1], Sizeof(Float)*a.Nx, 0);
  702.    End; { If }
  703.  
  704.   If a.Ny > 0 Then
  705.    Begin
  706.     GetMem( a.y, Sizeof(Float)*a.Ny);
  707.     FillChar( a.y^[1], Sizeof(Float)*a.Ny, 0);
  708.    End; { If }
  709. End; { Procedure GetStatMemory }
  710.  
  711. {═════════════════════════════════════════════════════════════════════════════}
  712.  
  713. Procedure FreeStatMemory( var a: StatRecord);
  714.  
  715. { free memory reserved by GetStatMemory. This procedure uses FreeMem which is
  716.   incompatible with Mark/Release. Use Dispose to get rid of all your other
  717.   pointers                                                                    }
  718.  
  719. Begin
  720.   If a.Nx > 0 Then FreeMem( a.x, Sizeof(Float)*a.Nx);
  721.   If a.Ny > 0 Then FreeMem( a.y, Sizeof(Float)*a.Ny);
  722. End; { Procedure FreeStatMemory }
  723.  
  724. {═════════════════════════════════════════════════════════════════════════════}
  725.  
  726. Procedure ComputeStats( var a: StatRecord);
  727.  
  728. { These routines came straight out of My HP 15C Owner's Manual p.206 }
  729.  
  730. var
  731.   i:      1..2048;
  732.   Sum_x:  Float;
  733.   Sum_x2: Float;
  734.   Sum_y:  Float;
  735.   Sum_y2: Float;
  736.   Sum_xy: Float;
  737.   M:      Float;
  738.   N:      Float;
  739.   P:      Float;
  740.  
  741. Begin
  742.  
  743.   Sum_x  := 0;
  744.   Sum_x2 := 0;
  745.   Sum_y  := 0;
  746.   Sum_y2 := 0;
  747.   Sum_xy := 0;
  748.  
  749.   If (a.Nx = 0) and (a.Ny = 0) Then                { Nothing to do! }
  750.     Exit;
  751.  
  752.   If a.Ny = 0 Then        { X Data only }
  753.    Begin
  754.     For i := 1 to a.Ny do
  755.      Begin
  756.       Sum_x  := Sum_x  + a.x^[i];
  757.       Sum_x2 := sum_x2 + Sqr(a.x^[i]);
  758.      End;
  759.     a.MeanX := Sum_x / A.Nx;
  760.     a.StdX  := Sqrt((Sum_x2 - Sqr(Sum_x)/A.Nx)/Pred(a.Nx));
  761.    End
  762.   Else If a.Nx = 0 Then         { Y Data only }
  763.    Begin
  764.     For i := 1 to a.Ny do
  765.      Begin
  766.       Sum_y  := Sum_y  + a.y^[i];
  767.       Sum_y2 := sum_y2 + Sqr(a.y^[i]);
  768.      End;
  769.     a.MeanY := Sum_y / A.Ny;
  770.     a.StdY  := Sqrt((Sum_y2 - Sqr(Sum_y)/A.Ny)/Pred(a.Ny));
  771.    End
  772.   Else If A.Nx <> a.Ny Then
  773.    Begin
  774.     Writeln( 'Number of x points = ', a.Nx);
  775.     Writeln( 'Number of y points = ', a.Ny);
  776.     Writeln( 'Nx and Ny must be same value or zero!.');
  777.     Halt(1);
  778.    End
  779.   Else
  780.    Begin
  781.     For i := 1 to a.Ny do
  782.      Begin
  783.       Sum_x  := Sum_x  + a.x^[i];
  784.       Sum_x2 := sum_x2 + Sqr(a.x^[i]);
  785.       Sum_y  := Sum_y  + a.y^[i];
  786.       Sum_y2 := sum_y2 + Sqr(a.y^[i]);
  787.       Sum_xy := Sum_xy + a.x^[i] * a.y^[i]
  788.      End;
  789.     a.MeanX := Sum_x / A.Nx;
  790.     a.MeanY := Sum_y / A.Ny;
  791.     M       := a.Nx*Sum_x2 - Sqr(Sum_x);
  792.     N       := a.Nx*Sum_y2 - Sqr(Sum_y);
  793.     P       := A.Nx*Sum_xy - Sum_x * Sum_y;
  794.     a.StdX  := Sqrt( M/(Int(A.nx)*Int(Pred(a.Nx))));
  795.     a.StdY  := Sqrt( N/(Int(A.ny)*Int(Pred(a.Ny))));
  796.     a.r     := P / Sqrt(M*N);
  797.    End; { If }
  798. End;
  799.  
  800. {═════════════════════════════════════════════════════════════════════════════}
  801.  
  802. var
  803.   Odd:           Boolean;
  804.   Value_Save:    Float;
  805.  
  806. Function Gauss( mean, StdDev: Float): Float;
  807.  
  808. { This function returns a sample from a Gaussian distribution with the
  809.   specified mean and standard deviation. It does it from the taking 2
  810.   independant samples from a uniform distribution. Randomize must be called
  811.   before using this routine. (I've included it in the unit initialization
  812.   portion).
  813.                                                                             }
  814.  
  815. const
  816.   sqrt2: Float = 1.414213562373095150;
  817.   pi2:   Float = 6.283185307179586230;
  818.  
  819. var
  820.   a, b: Float;
  821.   temp1, temp2:  Float;
  822.  
  823. Begin
  824.   If Odd Then
  825.    Begin
  826.     Odd := False;
  827.     Gauss := Value_save * StdDev + mean;
  828.     Exit;
  829.    End; { if }
  830.  
  831.    a := Random;
  832.    b := random;
  833.  
  834.    Temp1 := pi2 * a;
  835.    Temp2 := Sqrt2 * Sqrt( -ln(b));
  836.  
  837.    Gauss := Cos( Temp1) * temp2 * STDDev + mean;
  838.    Value_save := Sin( Temp1) * temp2;
  839.    Odd := True;
  840. End;
  841.  
  842. {═════════════════════════════════════════════════════════════════════════════}
  843.  
  844. Procedure MatrixMultiply         (    Dimen:             Integer;
  845.                                   Var InputMatrix1:      Matrix;
  846.                                   Var InputMatrix2:      Matrix;
  847.                                   Var ProductMatrix:     Matrix);
  848. Var I,J,K:Integer;
  849. Begin
  850.   For I:=1 to Dimen do
  851.     For J:=1 to Dimen do
  852.     Begin
  853.       ProductMatrix[I][J] := 0;
  854.       For K:=1 to Dimen do
  855.         ProductMatrix[I][J] := ProductMatrix[I][J] +
  856.                                InputMatrix1[I][K] * InputMatrix2[K][J];
  857.     End;
  858. End;
  859.  
  860. {═════════════════════════════════════════════════════════════════════════════}
  861.  
  862. Procedure MatrixAdd              (    Dimen:             Integer;
  863.                                   Var InputMatrix1:      Matrix;
  864.                                   Var InputMatrix2:      Matrix;
  865.                                   Var AddedMatrix:       Matrix);
  866. Var I,J:Integer;
  867. Begin
  868.   For I:=1 to Dimen do
  869.     For J:=1 to Dimen do
  870.         AddedMatrix[I][J] := InputMatrix1[I][J] + InputMatrix2[I][J];
  871. End;
  872.  
  873. {═════════════════════════════════════════════════════════════════════════════}
  874.  
  875. Procedure MatrixSubtract         (    Dimen:             Integer;
  876.                                   Var InputMatrix1:      Matrix;
  877.                                   Var InputMatrix2:      Matrix;
  878.                                   Var Matrix1Minus2:     Matrix);
  879. Var I,J:Integer;
  880. Begin
  881.   For I:=1 to Dimen do
  882.     For J:=1 to Dimen do
  883.         Matrix1Minus2[I][J] := InputMatrix1[I][J] - InputMatrix2[I][J];
  884. End;
  885.  
  886. {═════════════════════════════════════════════════════════════════════════════}
  887.  
  888. Procedure MatrixScalarMultiply   (    Dimen:             Integer;
  889.                                   Var InputMatrix:       Matrix;
  890.                                       Scalar:            Real;
  891.                                   Var ProductMatrix:     Matrix);
  892. Var I,J:Integer;
  893. Begin
  894.   For I:=1 to Dimen do
  895.     For J:=1 to Dimen do
  896.         ProductMatrix[I][J] := Scalar * InputMatrix[I][J];
  897. End;
  898.  
  899. {═════════════════════════════════════════════════════════════════════════════}
  900.  
  901. Procedure BigMatrixScalarMultiply(    Dimen:             Integer;
  902.                                   Var InputMatrix:       Matrix2;
  903.                                       Scalar:            Real;
  904.                                   Var ProductMatrix:     Matrix2);
  905. Var I,J:Integer;
  906. Begin
  907.   For I:=1 to Dimen do
  908.     For J:=1 to Dimen do
  909.         ProductMatrix[I][J] := Scalar * InputMatrix[I][J];
  910. End;
  911.  
  912. {═════════════════════════════════════════════════════════════════════════════}
  913.  
  914. Procedure BigPtrMatrixScalarMultiply(    Dimen:             Integer;
  915.                                   Var InputMatrix:       Matrix2ptr;
  916.                                       Scalar:            Real;
  917.                                   Var ProductMatrix:     Matrix2ptr);
  918. Var I,J:Integer;
  919. Begin
  920.   For I:=1 to Dimen do
  921.     For J:=1 to Dimen do
  922.         ProductMatrix^[I][J] := Scalar * InputMatrix^[I][J];
  923. End;
  924.  
  925. {═════════════════════════════════════════════════════════════════════════════}
  926.  
  927. Procedure MatrixInverse          (    Dimen:             Integer;
  928.                                   Var Input:             Matrix;
  929.                                   Var Output:            Matrix;
  930.                                   Var Error:             Byte);
  931. {This program was copied from an old fortran routine of unknown origin.  It
  932.  is better than the Turbo Numerical methods inversion routine in that it
  933.  can handle larger arrays.}
  934. Var
  935.   Determ:       double;
  936.   I,J,K,L:      Integer;
  937.   Ipivot:       array[1..Arraysize]  of integer;
  938.   pivot:        array[1..Arraysize]  of Real;
  939.   Index2:       Array[1..Arraysize, 1..2] of integer;
  940.   Irow:         Integer;
  941.   Icol:         Integer;
  942.   Swap:         Real;
  943.   Amax:         Real;
  944.   L1:           Integer;
  945.   T:            Real;
  946.   Jrow, Jcol:   Integer;
  947.  
  948. Begin
  949.   { Initialization }
  950. { for i:=1 to dimen do                                     { debug }
  951. {   for j:=1 to dimen do                                   { debug }
  952. {     writeln('input[',i:2,',',j:2,'] = ', input[i][j]); { debug }
  953.   MatrixScalarMultiply ( Dimen, Input, 1, Output); {copies input to output}
  954. { for i:=1 to dimen do                                     { debug }
  955. {   for j:=1 to dimen do                                   { debug }
  956. {     writeln('output[',i:2,',',j:2,'] = ', output[i][j]); { debug }
  957.   Error := 0;
  958.   Determ := 1.0;
  959.   For I:=1 to Dimen do Ipivot[I] := 0;
  960.  
  961.   { Search for pivot element }
  962.   For I:=1 to Dimen do
  963.   Begin
  964.     Amax := 0.0;
  965.     icol := -1;
  966.     For J:=1 to Dimen do
  967.       If Ipivot[J] <> 1 then
  968.       Begin
  969.         For K:=1 to Dimen do
  970.         Begin
  971.           If Ipivot[K] -1 > 0 then exit;
  972.           if (ipivot[K] -1 < 0) and ((abs(amax)-abs(Output[J][K]))<0) then
  973.           Begin
  974.             Irow := J;
  975.             Icol := K;
  976.             Amax := Output[J][K];
  977.           End;
  978.         End;
  979.       End;
  980.     if icol < 0 then
  981.     begin
  982.       error := 1;
  983.       writeln('Zero Matrix?');
  984.       exit;
  985.     end;
  986.     Ipivot[Icol] := Ipivot[Icol]+1;
  987.  
  988.     { Interchange rows to put pivot element on diagonal }
  989.     If Irow <> Icol then
  990.     Begin
  991.       Determ := -Determ;
  992.       For L:= 1 to Dimen do
  993.       Begin
  994.         Swap := Output[Irow][L];
  995.         Output[Irow][L] := Output[Icol][L];
  996.         Output[Icol][L] := Swap;
  997.       End;
  998.     End;
  999.     Index2[I][1] := Irow;
  1000.     Index2[I][2] := Icol;
  1001.     Pivot[I] := Output[Icol][Icol];
  1002. {writeln('determ, pivot[',i,']: ', determ, pivot[i]); { debug }
  1003.     Determ := Determ * Pivot[I];
  1004.  
  1005.     { Divide pivot row by pivot element }
  1006.     Output[Icol][Icol] := 1.0;
  1007.     For L:=1 to Dimen do
  1008.       If pivot[I] <> 0.0 then
  1009.         Output[Icol][L] := Output[Icol][L] / Pivot[I]
  1010.       Else
  1011.       Begin
  1012.         Error := 1;
  1013.         Exit;
  1014.       End;
  1015.  
  1016.     { Reduce non-pivot rows }
  1017.     For L1 := 1 to Dimen do
  1018.       If l1 <> Icol then
  1019.       Begin
  1020.         T := Output[L1][Icol];
  1021.         Output[L1][Icol] := 0.0;
  1022.         For L:=1 to Dimen do
  1023.           Output[L1][L] := Output[L1][L] - Output[Icol][L]*T;
  1024.       End;
  1025.   End;
  1026.  
  1027.   { Interchange columns }
  1028.   For I:=1 to Dimen do
  1029.   Begin
  1030.     L := Dimen + 1 - I;
  1031.     If Index2[L,1] <> Index2[L,2] then
  1032.     Begin
  1033.       Jrow := Index2[L,1];
  1034.       Jcol := Index2[L,2];
  1035.       For K:=1 to dimen do
  1036.       Begin
  1037.         Swap := Output[K][Jrow];
  1038.         Output[K][Jrow] := Output[K][Jcol];
  1039.         Output[K][Jcol] := Swap;
  1040.       End;
  1041.     End;
  1042.   End;
  1043.  
  1044. End;
  1045.  
  1046.  
  1047. {═════════════════════════════════════════════════════════════════════════════}
  1048.  
  1049.  
  1050. Procedure BigMatrixInverse       (    Dimen:             Integer;
  1051.                                   Var Input:             Matrix2;
  1052.                                   Var Output:            Matrix2;
  1053.                                   Var Error:             Byte);
  1054. {This program was copied from an old fortran routine of unknown origin.  It
  1055.  is better than the Turbo Numerical methods inversion routine in that it
  1056.  can handle larger arrays.}
  1057. Var
  1058.   Determ:       Double;
  1059.   I,J,K,L:      Integer;
  1060.   Ipivot:       array[1..Arraysize2]  of integer;
  1061.   pivot:        array[1..Arraysize2]  of Real;
  1062.   Index2:       Array[1..Arraysize2, 1..2] of integer;
  1063.   Irow:         Integer;
  1064.   Icol:         Integer;
  1065.   Swap:         Real;
  1066.   Amax:         Real;
  1067.   L1:           Integer;
  1068.   T:            Real;
  1069.   Jrow, Jcol:   Integer;
  1070.  
  1071. Begin
  1072.   { Initialization }
  1073. { for i:=1 to dimen do                                     { debug }
  1074. {   for j:=1 to dimen do                                   { debug }
  1075. {     writeln('input[',i:2,',',j:2,'] = ', input[i][j]); { debug }
  1076.   BigMatrixScalarMultiply ( Dimen, Input, 1, Output); {copies input to output}
  1077. { for i:=1 to dimen do                                     { debug }
  1078. {   for j:=1 to dimen do                                   { debug }
  1079. {     writeln('output[',i:2,',',j:2,'] = ', output[i][j]); { debug }
  1080.   Error := 0;
  1081.   Determ := 1.0;
  1082.   For I:=1 to Dimen do Ipivot[I] := 0;
  1083.  
  1084.   { Search for pivot element }
  1085.   For I:=1 to Dimen do
  1086.   Begin
  1087.     Amax := 0.0;
  1088.     icol := -1;
  1089.     For J:=1 to Dimen do
  1090.       If Ipivot[J] <> 1 then
  1091.       Begin
  1092.         For K:=1 to Dimen do
  1093.         Begin
  1094.           If Ipivot[K] -1 > 0 then exit;
  1095.           if (ipivot[K] -1 < 0) and ((abs(amax)-abs(Output[J][K]))<0) then
  1096.           Begin
  1097.             Irow := J;
  1098.             Icol := K;
  1099.             Amax := Output[J][K];
  1100.           End;
  1101.         End;
  1102.       End;
  1103.     if icol < 0 then
  1104.     begin
  1105.       error := 1;
  1106.       writeln('Zero Matrix?');
  1107.       exit;
  1108.     end;
  1109.     Ipivot[Icol] := Ipivot[Icol]+1;
  1110.  
  1111.     { Interchange rows to put pivot element on diagonal }
  1112.     If Irow <> Icol then
  1113.     Begin
  1114.       Determ := -Determ;
  1115.       For L:= 1 to Dimen do
  1116.       Begin
  1117.         Swap := Output[Irow][L];
  1118.         Output[Irow][L] := Output[Icol][L];
  1119.         Output[Icol][L] := Swap;
  1120.       End;
  1121.     End;
  1122.     Index2[I][1] := Irow;
  1123.     Index2[I][2] := Icol;
  1124.     Pivot[I] := Output[Icol][Icol];
  1125. {writeln('determ, pivot[',i,']: ', determ, pivot[i]); { debug }
  1126.     Determ := Determ * Pivot[I];
  1127.  
  1128.     { Divide pivot row by pivot element }
  1129.     Output[Icol][Icol] := 1.0;
  1130.     For L:=1 to Dimen do
  1131.       If pivot[I] <> 0.0 then
  1132.         Output[Icol][L] := Output[Icol][L] / Pivot[I]
  1133.       Else
  1134.       Begin
  1135.         Error := 1;
  1136.         Exit;
  1137.       End;
  1138.  
  1139.     { Reduce non-pivot rows }
  1140.     For L1 := 1 to Dimen do
  1141.       If l1 <> Icol then
  1142.       Begin
  1143.         T := Output[L1][Icol];
  1144.         Output[L1][Icol] := 0.0;
  1145.         For L:=1 to Dimen do
  1146.           Output[L1][L] := Output[L1][L] - Output[Icol][L]*T;
  1147.       End;
  1148.   End;
  1149.  
  1150.   { Interchange columns }
  1151.   For I:=1 to Dimen do
  1152.   Begin
  1153.     L := Dimen + 1 - I;
  1154.     If Index2[L,1] <> Index2[L,2] then
  1155.     Begin
  1156.       Jrow := Index2[L,1];
  1157.       Jcol := Index2[L,2];
  1158.       For K:=1 to dimen do
  1159.       Begin
  1160.         Swap := Output[K][Jrow];
  1161.         Output[K][Jrow] := Output[K][Jcol];
  1162.         Output[K][Jcol] := Swap;
  1163.       End;
  1164.     End;
  1165.   End;
  1166.  
  1167. End;
  1168.  
  1169. {═════════════════════════════════════════════════════════════════════════════}
  1170.  
  1171.  
  1172. Procedure BigPtrMatrixInverse    (    Dimen:             Integer;
  1173.                                   Var Input:             Matrix2ptr;
  1174.                                   Var Output:            Matrix2ptr;
  1175.                                   Var Error:             Byte);
  1176. {This program was copied from an old fortran routine of unknown origin.  It
  1177.  is better than the Turbo Numerical methods inversion routine in that it
  1178.  can handle larger arrays.}
  1179. Var
  1180.   Determ:       Double;
  1181.   I,J,K,L:      Integer;
  1182.   Ipivot:       array[1..Arraysize2]  of integer;
  1183.   pivot:        array[1..Arraysize2]  of Real;
  1184.   Index2:       Array[1..Arraysize2, 1..2] of integer;
  1185.   Irow:         Integer;
  1186.   Icol:         Integer;
  1187.   Swap:         Real;
  1188.   Amax:         Real;
  1189.   L1:           Integer;
  1190.   T:            Real;
  1191.   Jrow, Jcol:   Integer;
  1192.  
  1193. Begin
  1194.   { Initialization }
  1195. { for i:=1 to dimen do                                     { debug }
  1196. {   for j:=1 to dimen do                                   { debug }
  1197. {     writeln('input[',i:2,',',j:2,'] = ', input[i][j]); { debug }
  1198.   BigPtrMatrixScalarMultiply ( Dimen, Input, 1, Output); {copies input to output}
  1199. { for i:=1 to dimen do                                     { debug }
  1200. {   for j:=1 to dimen do                                   { debug }
  1201. {     writeln('output[',i:2,',',j:2,'] = ', output[i][j]); { debug }
  1202.   Error := 0;
  1203.   Determ := 1.0;
  1204.   For I:=1 to Dimen do Ipivot[I] := 0;
  1205.  
  1206.   { Search for pivot element }
  1207.   For I:=1 to Dimen do
  1208.   Begin
  1209.     Amax := 0.0;
  1210.     icol := -1;
  1211.     For J:=1 to Dimen do
  1212.       If Ipivot[J] <> 1 then
  1213.       Begin
  1214.         For K:=1 to Dimen do
  1215.         Begin
  1216.           If Ipivot[K] -1 > 0 then exit;
  1217.           if (ipivot[K] -1 < 0) and ((abs(amax)-abs(Output^[J][K]))<0) then
  1218.           Begin
  1219.             Irow := J;
  1220.             Icol := K;
  1221.             Amax := Output^[J][K];
  1222.           End;
  1223.         End;
  1224.       End;
  1225.     if icol < 0 then
  1226.     begin
  1227.       error := 1;
  1228.       writeln('Zero Matrix?');
  1229.       exit;
  1230.     end;
  1231.     Ipivot[Icol] := Ipivot[Icol]+1;
  1232.  
  1233.     { Interchange rows to put pivot element on diagonal }
  1234.     If Irow <> Icol then
  1235.     Begin
  1236.       Determ := -Determ;
  1237.       For L:= 1 to Dimen do
  1238.       Begin
  1239.         Swap := Output^[Irow][L];
  1240.         Output^[Irow][L] := Output^[Icol][L];
  1241.         Output^[Icol][L] := Swap;
  1242.       End;
  1243.     End;
  1244.     Index2[I][1] := Irow;
  1245.     Index2[I][2] := Icol;
  1246.     Pivot[I] := Output^[Icol][Icol];
  1247. {writeln('determ, pivot[',i,']: ', determ, pivot[i]); { debug }
  1248.     Determ := Determ * Pivot[I];
  1249.  
  1250.     { Divide pivot row by pivot element }
  1251.     Output^[Icol][Icol] := 1.0;
  1252.     For L:=1 to Dimen do
  1253.       If pivot[I] <> 0.0 then
  1254.         Output^[Icol][L] := Output^[Icol][L] / Pivot[I]
  1255.       Else
  1256.       Begin
  1257.         Error := 1;
  1258.         Exit;
  1259.       End;
  1260.  
  1261.     { Reduce non-pivot rows }
  1262.     For L1 := 1 to Dimen do
  1263.       If l1 <> Icol then
  1264.       Begin
  1265.         T := Output^[L1][Icol];
  1266.         Output^[L1][Icol] := 0.0;
  1267.         For L:=1 to Dimen do
  1268.           Output^[L1][L] := Output^[L1][L] - Output^[Icol][L]*T;
  1269.       End;
  1270.   End;
  1271.  
  1272.   { Interchange columns }
  1273.   For I:=1 to Dimen do
  1274.   Begin
  1275.     L := Dimen + 1 - I;
  1276.     If Index2[L,1] <> Index2[L,2] then
  1277.     Begin
  1278.       Jrow := Index2[L,1];
  1279.       Jcol := Index2[L,2];
  1280.       For K:=1 to dimen do
  1281.       Begin
  1282.         Swap := Output^[K][Jrow];
  1283.         Output^[K][Jrow] := Output^[K][Jcol];
  1284.         Output^[K][Jcol] := Swap;
  1285.       End;
  1286.     End;
  1287.   End;
  1288.  
  1289. End;
  1290.  
  1291.  
  1292. {═════════════════════════════════════════════════════════════════════════════}
  1293.  
  1294. Procedure ComplexMatrixInverse   (    Dimen:             Integer;
  1295.                                   Var RealInputMatrix:   Matrix;
  1296.                                   Var ImagInputMatrix:   Matrix;
  1297.                                   Var RealInverseMatrix: Matrix;
  1298.                                   Var ImagInverseMatrix: Matrix;
  1299.                                   Var Error:             Byte);
  1300.  
  1301.  
  1302. Var
  1303.   Phase:    Complex;
  1304.   I,J,K:    Integer;
  1305.   big:      Matrix2ptr;
  1306.   bigout:   Matrix2ptr;
  1307.  
  1308. begin
  1309.   new(big);
  1310.   new(bigout);
  1311.   for i:=1 to Dimen do
  1312.     for j:=1 to Dimen do
  1313.     begin
  1314.       big^[i,j] := RealInputMatrix[i,j];
  1315.       big^[i+Dimen, j+Dimen] := RealInputMatrix[i,j];
  1316.     end;
  1317.   for i:=1 to Dimen do
  1318.     for j:=1 to Dimen do
  1319.     begin
  1320.       big^[i, j+Dimen] := ImagInputMatrix[i,j];
  1321.       big^[i+Dimen, j] := ImagInputMatrix[i,j];
  1322.     end;
  1323.     BigPtrMatrixInverse ( Dimen * 2, big, bigout, error);
  1324.     for i:=1 to Dimen do
  1325.       for j:=1 to Dimen do
  1326.       begin
  1327.         realinversematrix[i,j] :=  bigout^[i,j];
  1328.         imaginversematrix[i,j] := -bigout^[i,j+Dimen];
  1329.       end;
  1330.     dispose(bigout);
  1331.     dispose(big);
  1332.   end;
  1333.  
  1334. {═════════════════════════════════════════════════════════════════════════════}
  1335.  
  1336.  
  1337. Begin
  1338.   Randomize;
  1339.   odd := False;
  1340. End.
  1341.