ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng)/ProfitPress-MegaCDROM2.B6I/PROG/PASCAL/TPMATH.ZIP/MATH.PAS

Pascal/Delphi Source File | 1989-12-12 | 40.5 KB | 1,341 lines

Unit Math;
Interface

(*
 ╔════════════════════════════════════════════════════════════════════════════╗
 ║ Standard Mathematical Functions.        Version 1.0                        ║
 ║                                                                            ║
 ║ This code has been placed in the public domain.                            ║
 ║ All non-trivial functions that do not reside permanently inside my head    ║
 ║ (e.g. Bessel Functions) were obtain from:                                  ║
 ║                                                                            ║
 ║   _CRC_Standard_Mathematical_Tables_, CRC Press, 25th Edition, William H.  ║
 ║   Beyer, editor.                                                           ║
 ║                                                                            ║
 ║   HP-15C Owner's Manual, Feb. 1984, Hewlett-Packard Company.               ║
 ║                                                                            ║
 ║   _Probability,_Random_Variables,_and_Stochastic_Processes_, Athanasios    ║
 ║   Papoulis, Second Edition, 1984, McGraw-Hill                              ║
 ║                                                                            ║
 ║                                                                            ║
 ║ All code has been tested, though not rigorusly, and is correct to the best ║
 ║ of my knowledge. However, I make no guarentees. If you find any errors or  ║
 ║ make any changes/enhancements, I would appreciate it if you inform me. My  ║
 ║ address is:                                                                ║
 ║                                                                            ║
 ║ Scott Bostater                                                             ║
 ║ Georgia Tech Research Institute                                            ║
 ║ 7220 Richardson Road                                                       ║
 ║ Symrna, GA 30080                                                           ║
 ║ (404)-528-7782    (please, only in an emergency, my boss to work for HIM)  ║
 ║                                                                            ║
 ║ e-mail:                                                                    ║
 ║ ---------------------------------------                                    ║
 ║ Internet: [e-mail]                                            ║
 ║ UUCP:     ...!{allegra,amd,hplabs,ut-ngp}!gatech!edu!bb16                  ║
 ║ Bitnet:   BBOSTATE@GTRI01                                                  ║
 ║                                                                            ║
 ║                                                                            ║
 ║ Special Thanks to Mike Baden and Stan West for thier contributions         ║
 ║                                                                            ║
 ╚════════════════════════════════════════════════════════════════════════════╝
*)

{ Delcare Standard Math types}
Type
  {$IFOPT N+}
    Float = Double;
  {$Else}
    Float = Real;
  {$ENDIF}
  Complex = Record
              Real: Float;
              Imag: Float;
            End;

{*
 * Run of the mill math functions that I always need
 *}

Function  Log                 (x:Float):Float;                  { Log base 10 }
Function  Power               (x,y:Float):Float;
Function  Atan360             (x,y:Float):Float;
Function  Tan                 (x:Float):Float;
Function  ArcSin              (x:Float):Float;
Function  ArcCos              (x:Float):Float;
Function  HypSin              (x:Float):Float;
Function  HypCos              (x:Float):Float;

Function  Inverse_tangent_deg (X,Y:Float):Float;
Function  Inverse_tangent_rad (X,Y:Float):Float;
Function  Inverse_sin_deg     (Y,R:Float):Float;
Function  Inverse_sin_rad     (Y,R:Float):Float;
Function  Inverse_cos_deg     (X,R:Float):Float;
Function  Inverse_cos_rad     (X,R:Float):Float;

{*
 * Complex math
 *}

Procedure Complex_add         (A,B:Complex; Var C:Complex);
Procedure Complex_subtract    (A,B:Complex; Var C:Complex);
Procedure Complex_multiply    (A,B:Complex; Var C:Complex);
Procedure Complex_divide      (A,B:Complex; Var C:Complex);
Function  Magnitude           (A:Complex):Float;
Procedure complex_conjugate   (A:Complex; var c:complex);


{*
 * Bessel Functions, (I hate Bessel Functions!)
 *}

function  I0                  (x:Float):Float;   { modified Bessel of Order 0 }
function  I1                  (x:Float):Float;   { modified Bessel of Order 1 }
function  I2                  (x:Float):Float;   { modified Bessel of Order 2 }
function  J0                  (x:Float):Float;   { Bessel function of Order 0 }
function  J1                  (x:Float):Float;   { Bessel function of Order 1 }

{*
 * Family of Unit Step functions
 *
 * U(0,t)  = 1 (t=0), = 0 (else)
 * U(-1,t) = 1 (t>0), 1/2 (t=0), 0 (t<0)     [integral of U(0,t)]
 * U(-2,t) = t (t>=0), 0 (t<0)               [integral of U(-1,t)]
 * U(-3,t) = ½t² (t>=0), 0 (t<0)             [integral of U(-2,t)]
 *}

Function  U                   (index:Integer; T:Float):Float;

{*
 * Voltage and Power Decibel Functions
 *}

Function  db10                (x:Float):Float;
Function  db20                (x:Float):Float;
Function  undb10              (X:Float):Float;
Function  undb20              (x:Float):Float;

{*
 * Hex output Functions  (TP really should have had these built in)
 *}

type
  String9 = String[9];
  String4 = String[4];
  String2 = String[2];

Function  HexLong             (a:LongInt):String9;
Function  HexWord             (a:Word):   String4;
Function  HexInt              (a:Integer):String4;
Function  HexByte             (a:Byte):   String2;

{*
 * If you have an 80387, These commands are _much_ better than the standard
 * Sin/Cos functions that TP provides
 *}

{$Ifdef P387 }
Function Cos( x: Float): Float;
Function Sin( x: Float): Float;
{$EndIf}

{*
 * Statistical Routines
 *}
type
  StatsArray = Array[1..2048] of Float;
  StatRecord = Record
                 Nx:    0..2048;                    { Number of x points      }
                 Ny:    0..2048;                    { Number of y points      }
                 x:     ^StatsArray;                { x data                  }
                 y:     ^StatsArray;                { y data                  }
                 MeanX: Float;                      { X average               }
                 MeanY: Float;                      { Y average               }
                 StdX:  Float;                      { Standard Deviation of x }
                 StdY:  Float;                      { Standard Deviation of y }
                 r:     Float;                      { correlation coefficient }
               End; { Record }

Procedure GetStatMemory(  var a: StatRecord);
Procedure FreeStatMemory( var a: StatRecord);
Procedure ComputeStats(   var a: StatRecord);       { mean, std dev, cor.coef }
Function  Gauss( mean, StdDev: Float):Float;        { sample from Gauss dist.}

{*
 * Matirx Manipulations
 *}

Const
  Arraysize = 32;
  Arraysize2= 64;

Type
  Vector   = array[0..ArraySize]  of Real;
  Vector2  = array[0..ArraySize2] of real;
  Matrix   = array[0..ArraySize] of Vector;
  Matrix2  = array[0..ArraySize2] of Vector2;
  Matrix2ptr = ^Matrix2;

Procedure MatrixMultiply         (    Dimen:             Integer;
                                  Var InputMatrix1:      Matrix;
                                  Var InputMatrix2:      Matrix;
                                  Var ProductMatrix:     Matrix);
Procedure MatrixAdd              (    Dimen:             Integer;
                                  Var InputMatrix1:      Matrix;
                                  Var InputMatrix2:      Matrix;
                                  Var AddedMatrix:       Matrix);
Procedure MatrixSubtract         (    Dimen:             Integer;
                                  Var InputMatrix1:      Matrix;
                                  Var InputMatrix2:      Matrix;
                                  Var Matrix1Minus2:     Matrix);
Procedure MatrixScalarMultiply   (    Dimen:             Integer;
                                  Var InputMatrix:       Matrix;
                                      Scalar:            Real;
                                  Var ProductMatrix:     Matrix);
Procedure BigMatrixScalarMultiply(    Dimen:             Integer;
                                  Var InputMatrix:       Matrix2;
                                      Scalar:            Real;
                                  Var ProductMatrix:     Matrix2);
Procedure BigPtrMatrixScalarMultiply(    Dimen:             Integer;
                                  Var InputMatrix:       Matrix2ptr;
                                      Scalar:            Real;
                                  Var ProductMatrix:     Matrix2ptr);
Procedure ComplexMatrixInverse   (    Dimen:             Integer;
                                  Var RealInputMatrix:   Matrix;
                                  Var ImagInputMatrix:   Matrix;
                                  Var RealInverseMatrix: Matrix;
                                  Var ImagInverseMatrix: Matrix;
                                  Var Error:             Byte);
Procedure MatrixInverse          (    Dimen:             Integer;
                                  Var Input:             Matrix;
                                  Var Output:            Matrix;
                                  Var Error:             Byte);
Procedure BigMatrixInverse       (    Dimen:             Integer;
                                  Var Input:             Matrix2;
                                  Var Output:            Matrix2;
                                  Var Error:             Byte);
Procedure BigPtrMatrixInverse    (    Dimen:             Integer;
                                  Var Input:             Matrix2Ptr;
                                  Var Output:            Matrix2Ptr;
                                  Var Error:             Byte);


Implementation
const
  Minus2: Integer = -2;



{$IFDef P387 }                                    { If 80387 available, load }
Function Cos( x: Float): Float; external;         { in 80387 specific        }
Function Sin( x: Float): Float; external;         { functions}
Function Tan( x: Float): Float; external;
{$L trig387.obj }

{$ELSE}
{$IFOPT N+}
Function Tan( X: Float): Float; external;         { 8087 or 80287 Presend     }
                                                  { Load in Assembly routine  }
{$L trig.obj }                                    { for fast Tan(x)           }

{$Else}
Function Tan( X: Float): Float;                   { No Coprocessor present.   }
Begin                                             { Out of luck. Do it the old}
  Tan := Sin(x)/Cos(x);                           { fashioned way             }
End;
{$ENDIF}
{$EndIf}


{═════════════════════════════════════════════════════════════════════════════}

Function U( index: Integer; T: Float): Float;
var
  P: ShortInt;

Begin
  P := -1-index;
  U := 0;              { default }

  If T = 0 Then
   Begin
    If index = 0 Then U := 1
   End
  Else If T > 0 Then
    Case P of
     0: U := 1;
     1: U := T;
     2: U := Sqr(t)/2;
     3..127: U := Power(t,P)/P;
    End; { Case }
End; { Function U }

{═════════════════════════════════════════════════════════════════════════════}

Function Magnitude (A:Complex):Float;
Begin
  Magnitude := Sqrt(Sqr(A.Real) + Sqr(A.Imag));
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure complex_conjugate  (A:Complex; var c:complex);
Begin
  C.Real := A.Real;
  C.Imag :=-A.Imag;
End;

{═════════════════════════════════════════════════════════════════════════════}

Function Log (X:Float):Float;
Const ln10 = 2.302585093;
Begin
  Log := ln(X) / ln10;
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure Complex_Add (A,B:Complex; Var C:Complex);
{ C = A + B }
Begin
  C.Real := A.Real + B.Real;
  C.Imag := A.Imag + B.Imag;
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure Complex_subtract (A,B:Complex; Var C:Complex);
{ C = A - B }
Begin
  C.Real := A.Real - B.Real;
  C.Imag := A.Imag - B.Imag;
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure Complex_multiply (A,B:Complex; Var C:Complex);
{ C = A * B }
Begin
  C.Real := A.Real * B.Real  -  A.Imag * B.Imag;
  C.Imag := A.Real * B.Imag  +  A.Imag * B.Real;
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure Complex_divide (A,B:Complex; Var C:Complex);
{ C = A / B }
Var
  Temp:Real;
Begin
  Temp   :=  B.Real * B.Real  +  B.Imag * B.Imag;
  C.Real := (A.Real * B.Real  +  A.Imag * B.Imag) / Temp;
  C.Imag := (A.Imag * B.Real  -  A.Real * B.Imag) / Temp;
End;

{═════════════════════════════════════════════════════════════════════════════}

Function Inverse_tangent_deg(X,Y:Float):Float;
{ 0 <= result <= 360 }
Var
  Angle:Float;
Begin
  if X=0 then Angle := Pi/2.0 else angle := ArcTan(Y/X);
  angle := angle*180.0/pi;
  if (X<=0.0) and (Y<0) then angle := angle - 180.0;
  if (X< 0.0) and (Y>0) then angle := angle + 180.0;
  If angle < 0 then angle := angle+360.0;
  Inverse_tangent_deg := angle;
End;

{═════════════════════════════════════════════════════════════════════════════}

Function Inverse_tangent_rad(X,Y:Float):Float;
{ 0 <= result <= 2pi }
Var
  Angle:Float;
Begin
  if X=0 then Angle := Pi/2.0 else angle := ArcTan(Y/X);
  if (X<=0.0) and (Y<0) then angle := angle - pi;
  if (X< 0.0) and (Y>0) then angle := angle + pi;
  If angle < 0 then angle := angle+2*pi;
  Inverse_tangent_rad := angle;
End;

{═════════════════════════════════════════════════════════════════════════════}

Function Inverse_sin_deg (Y,R:Float):Float;
{ -90 <= result <= 90 }
Var
  X:Float;
  temp:Float;
Begin
  X := Sqrt(Sqr(R)-Sqr(Y));
  Temp := Inverse_tangent_deg(X,Y);
  If Temp > 90 then Temp := Temp - 360;
  Inverse_sin_deg := Temp;
End;

{═════════════════════════════════════════════════════════════════════════════}

Function Inverse_sin_rad (Y,R:Float):Float;
{ -90 <= result <= 90 }
Var
  X:Float;
  temp:Float;
Begin
  X := Sqrt(Sqr(R)-Sqr(Y));
  Temp := Inverse_tangent_rad(X,Y);
  If Temp > 90 then Temp := Temp - 360;
  Inverse_sin_rad := Temp;
End;

{═════════════════════════════════════════════════════════════════════════════}

Function Inverse_cos_deg (X,R:Float):Float;
{ -90 <= result <= 90 }
Var
  Y:Float;
  temp:Float;
Begin
  Y := Sqrt(Sqr(R)-Sqr(X));
  Temp := Inverse_tangent_deg(X,Y);
  If Temp > 90 then Temp := Temp - 360;
  Inverse_cos_deg := Temp;
End;

{═════════════════════════════════════════════════════════════════════════════}

Function Inverse_cos_rad (X,R:Float):Float;
{ -90 <= result <= 90 }
Var
  Y:Float;
  temp:Float;
Begin
  Y := Sqrt(Sqr(R)-Sqr(X));
  Temp := Inverse_tangent_rad(X,Y);
  If Temp > 90 then Temp := Temp - 360;
  Inverse_cos_rad := Temp;
End;

{═════════════════════════════════════════════════════════════════════════════}

Function ArcSin( x: Float): Float;
Begin
  ArcSin := Arctan( x / Sqrt( 1-Sqr(x)));
End;

{═════════════════════════════════════════════════════════════════════════════}

const
  pi4 = 0.785398164;

Function ArcCos( x: Float): Float;
Begin
  If x > 1 Then
    ArcCos := 0
  Else if x < -1 Then
    ArcCos := pi
  Else
    ArcCos := pi4 - Arctan( x / Sqrt( 1-Sqr(x)));
End;

{═════════════════════════════════════════════════════════════════════════════}

Function HypSin( x: Float): Float;
Begin
  Hypsin := 0.5 * (exp(x) - exp(-x));
End;

{═════════════════════════════════════════════════════════════════════════════}

Function HypCos( x: Float): Float;
Begin
  HypCos := 0.5 * (exp(x) + exp(-x));
End;

{═════════════════════════════════════════════════════════════════════════════}

function I0(x: Float): Float;
var
  tnum,factor: word;
  term,i0temp: Float;
begin
  term:=1;
  i0temp:=term;
  factor:=0;
  tnum:=1;
  repeat
    factor:=factor+2+2*(tnum-1);
    term:=term*sqr(x/2)/factor;
    i0temp:=i0temp+term;
    inc(tnum);
  until (abs(term)<1e-15) or (tnum=0);      {no overflow}
  i0:=i0temp;
end; {I0}

{═════════════════════════════════════════════════════════════════════════════}

function I1(x: Float): Float;
var
  tnum,factor: word;
  term,i1temp: Float;
begin
  term:=x/4;
  i1temp:=term;
  factor:=0;
  tnum:=1;
  repeat
    factor:=factor+3+2*(tnum-1);
    term:=term*sqr(x/2)/factor;
    writeln(term);
    i1temp:=i1temp+term;
    inc(tnum);
  until (abs(term)<1e-15) or (tnum=0);      {no overflow}
  i1:=i1temp;
end; {I0}

{═════════════════════════════════════════════════════════════════════════════}

function I2(x: Float): Float;
var
  tnum,factor: word;
  term,i2temp: Float;
begin
  term:=sqr(x)/24;
  i2temp:=term;
  factor:=0;
  tnum:=1;
  repeat
    factor:=factor+4+2*(tnum-1);
    term:=term*sqr(x/2)/factor;
    writeln(term);
    i2temp:=i2temp+term;
    inc(tnum);
  until (abs(term)<1e-15) or (tnum=0);      {no overflow}
  i2:=i2temp;
end; {I0}

{═════════════════════════════════════════════════════════════════════════════}

function J0(x: Float): Float;
var
  tnum: word;
  term,j0temp: Float;
begin
  j0temp:=1;
  term:=1;
  tnum:=1;
  repeat
    term:=term*-sqr(x)/(4*sqr(tnum));
    j0temp:=j0temp+term;
    inc(tnum);
  until (abs(term)<1e-15) or (tnum=0);      {no overflow}
  j0:=j0temp;
end; {J0}

{═════════════════════════════════════════════════════════════════════════════}

function J1(x: Float): Float;
var
  tnum,factor: word;
  term,j1temp: Float;
begin
  term:=x/2;
  j1temp:=term;
  factor:=0;
  tnum:=1;
  repeat
    factor:=factor+2+2*(tnum-1);
    term:=term*-sqr(x)/(4*factor);
    j1temp:=j1temp+term;
    inc(tnum);
  until (abs(term)<1e-15) or (tnum=0);      {no overflow}
  j1:=j1temp;
end; {J1}

{═════════════════════════════════════════════════════════════════════════════}

{═════════════════════════════════════════════════════════════════════════════}

Function Power( x,y: Float): Float;
Begin
  If y = 0 Then
    power := 1.0
  Else if x = 0 Then
    Power := 0.0
  Else If x > 0 Then
    Power := exp( y * ln(x))
  Else if Trunc(y) mod 2 = 0 Then
    Power := exp( y * ln(abs(x)))
  Else
    Power := -exp( y * ln(abs(x)));
End;

{═════════════════════════════════════════════════════════════════════════════}

Function Atan360( x, y: Float): Float;
Var
  Angle: Float;
Begin
  if X=0 then Angle := Pi/2.0 else angle := ArcTan(Y/X);
  angle := angle*180.0/pi;
  if (X<=0.0) and (Y<0) then angle := angle - 180.0;
  if (X< 0.0) and (Y>0) then angle := angle + 180.0;
  If angle < 0 then angle := angle+360.0;
  Atan360 := angle;
End;

{═════════════════════════════════════════════════════════════════════════════}

{$IFDEF N+}       { Limits change depending on whter 6 byte or 8 byte real }
const
  ln10: Float  = 2.30258509299405E+0000;

Function db10( X: Float): Float;
Begin
  If x < 1.0e-50 Then
    db10 := -500
  Else
    db10 := 10 * ln(x) / ln10;
End; { Function db10 }

{═════════════════════════════════════════════════════════════════════════════}

Function db20( X: Float): Float;
Begin
  If x < 1.0e-50 Then
    db20 := -1000
  Else
    db20 := 20 * ln(x) / ln10;
End; { Function db20 }


{═════════════════════════════════════════════════════════════════════════════}
{$ELSE}         { using 6 byte reals !}
const
  ln10: Float = 2.30258509299405E+0000;

Function db10( X: Float): Float;
Begin
  If x < 1.0e-25 Then
    db10 := -250
  Else
    db10 := 10 * ln(x) / ln10;
End; { Function db10 }

{═════════════════════════════════════════════════════════════════════════════}

Function db20( X: Float): Float;
Begin
  If x < 1.0e-25 Then
    db20 := -500
  Else
    db20 := 20 * ln(x) / ln10;
End; { Function db20 }

{$ENDIF}
{═════════════════════════════════════════════════════════════════════════════}

Function Undb10( X: Float): Float;
Begin
  Undb10 := Power(10, X/10);
End; { Function db10 }

{═════════════════════════════════════════════════════════════════════════════}

Function Undb20( X: Float): Float;
Begin
  Undb20 := Power(10, X/20);
End; { Function db20 }

{═════════════════════════════════════════════════════════════════════════════}

const
  h: Array[0..15] of Char = ( '0', '1', '2', '3', '4', '5', '6', '7',
                              '8', '9', 'A', 'B', 'C', 'D', 'E', 'F');

Function HexLong( a: LongInt): String9;
Begin
  HexLong := H[a shr 28] + H[(a shr 24) and $f] + H[ (a shr 20) and $f] +
             H[ (a shr 16) and $f] + ' ' + H[(a shr 12) and $f] +
             H[ (a shr 8) and $f] + H[ (a shr 4) and $f] + H[ a and $f];
End;

{═════════════════════════════════════════════════════════════════════════════}

Function HexWord( a: Word): String4;
Begin
  HexWord := H[a shr 12] + H[ (a shr 8) and $f] + H[ (a shr 4) and $f] +
             H[ a and $f];
End;

{═════════════════════════════════════════════════════════════════════════════}

Function HexInt( a: Integer): String4;
Begin
  HexInt := H[a shr 12] + H[ (a shr 8) and $f] + H[ (a shr 4) and $f] +
            H[ a and $f];
End;

{═════════════════════════════════════════════════════════════════════════════}

Function HexByte( a: Byte): String2;
Begin
  HexByte := H[ a shr 4 and $f] + H[ a and $f];
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure GetStatMemory( var a: StatRecord );

{ get free memory for a.x^ and a.y^ for the number of points specified in
  a.Nx and a.Ny.  It also clears out the previous contents of meanX...CovXY  }

Begin

  FillChar( a.meanX, Sizeof(Float)*5, 0);             { clear Stats data }

  If a.Nx > 0 Then
   Begin
    GetMem( a.x, Sizeof(Float)*a.Nx);
    FillChar( a.x^[1], Sizeof(Float)*a.Nx, 0);
   End; { If }

  If a.Ny > 0 Then
   Begin
    GetMem( a.y, Sizeof(Float)*a.Ny);
    FillChar( a.y^[1], Sizeof(Float)*a.Ny, 0);
   End; { If }
End; { Procedure GetStatMemory }

{═════════════════════════════════════════════════════════════════════════════}

Procedure FreeStatMemory( var a: StatRecord);

{ free memory reserved by GetStatMemory. This procedure uses FreeMem which is
  incompatible with Mark/Release. Use Dispose to get rid of all your other
  pointers                                                                    }

Begin
  If a.Nx > 0 Then FreeMem( a.x, Sizeof(Float)*a.Nx);
  If a.Ny > 0 Then FreeMem( a.y, Sizeof(Float)*a.Ny);
End; { Procedure FreeStatMemory }

{═════════════════════════════════════════════════════════════════════════════}

Procedure ComputeStats( var a: StatRecord);

{ These routines came straight out of My HP 15C Owner's Manual p.206 }

var
  i:      1..2048;
  Sum_x:  Float;
  Sum_x2: Float;
  Sum_y:  Float;
  Sum_y2: Float;
  Sum_xy: Float;
  M:      Float;
  N:      Float;
  P:      Float;

Begin

  Sum_x  := 0;
  Sum_x2 := 0;
  Sum_y  := 0;
  Sum_y2 := 0;
  Sum_xy := 0;

  If (a.Nx = 0) and (a.Ny = 0) Then                { Nothing to do! }
    Exit;

  If a.Ny = 0 Then        { X Data only }
   Begin
    For i := 1 to a.Ny do
     Begin
      Sum_x  := Sum_x  + a.x^[i];
      Sum_x2 := sum_x2 + Sqr(a.x^[i]);
     End;
    a.MeanX := Sum_x / A.Nx;
    a.StdX  := Sqrt((Sum_x2 - Sqr(Sum_x)/A.Nx)/Pred(a.Nx));
   End
  Else If a.Nx = 0 Then         { Y Data only }
   Begin
    For i := 1 to a.Ny do
     Begin
      Sum_y  := Sum_y  + a.y^[i];
      Sum_y2 := sum_y2 + Sqr(a.y^[i]);
     End;
    a.MeanY := Sum_y / A.Ny;
    a.StdY  := Sqrt((Sum_y2 - Sqr(Sum_y)/A.Ny)/Pred(a.Ny));
   End
  Else If A.Nx <> a.Ny Then
   Begin
    Writeln( 'Number of x points = ', a.Nx);
    Writeln( 'Number of y points = ', a.Ny);
    Writeln( 'Nx and Ny must be same value or zero!.');
    Halt(1);
   End
  Else
   Begin
    For i := 1 to a.Ny do
     Begin
      Sum_x  := Sum_x  + a.x^[i];
      Sum_x2 := sum_x2 + Sqr(a.x^[i]);
      Sum_y  := Sum_y  + a.y^[i];
      Sum_y2 := sum_y2 + Sqr(a.y^[i]);
      Sum_xy := Sum_xy + a.x^[i] * a.y^[i]
     End;
    a.MeanX := Sum_x / A.Nx;
    a.MeanY := Sum_y / A.Ny;
    M       := a.Nx*Sum_x2 - Sqr(Sum_x);
    N       := a.Nx*Sum_y2 - Sqr(Sum_y);
    P       := A.Nx*Sum_xy - Sum_x * Sum_y;
    a.StdX  := Sqrt( M/(Int(A.nx)*Int(Pred(a.Nx))));
    a.StdY  := Sqrt( N/(Int(A.ny)*Int(Pred(a.Ny))));
    a.r     := P / Sqrt(M*N);
   End; { If }
End;

{═════════════════════════════════════════════════════════════════════════════}

var
  Odd:           Boolean;
  Value_Save:    Float;

Function Gauss( mean, StdDev: Float): Float;

{ This function returns a sample from a Gaussian distribution with the
  specified mean and standard deviation. It does it from the taking 2
  independant samples from a uniform distribution. Randomize must be called
  before using this routine. (I've included it in the unit initialization
  portion).
                                                                            }

const
  sqrt2: Float = 1.414213562373095150;
  pi2:   Float = 6.283185307179586230;

var
  a, b: Float;
  temp1, temp2:  Float;

Begin
  If Odd Then
   Begin
    Odd := False;
    Gauss := Value_save * StdDev + mean;
    Exit;
   End; { if }

   a := Random;
   b := random;

   Temp1 := pi2 * a;
   Temp2 := Sqrt2 * Sqrt( -ln(b));

   Gauss := Cos( Temp1) * temp2 * STDDev + mean;
   Value_save := Sin( Temp1) * temp2;
   Odd := True;
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure MatrixMultiply         (    Dimen:             Integer;
                                  Var InputMatrix1:      Matrix;
                                  Var InputMatrix2:      Matrix;
                                  Var ProductMatrix:     Matrix);
Var I,J,K:Integer;
Begin
  For I:=1 to Dimen do
    For J:=1 to Dimen do
    Begin
      ProductMatrix[I][J] := 0;
      For K:=1 to Dimen do
        ProductMatrix[I][J] := ProductMatrix[I][J] +
                               InputMatrix1[I][K] * InputMatrix2[K][J];
    End;
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure MatrixAdd              (    Dimen:             Integer;
                                  Var InputMatrix1:      Matrix;
                                  Var InputMatrix2:      Matrix;
                                  Var AddedMatrix:       Matrix);
Var I,J:Integer;
Begin
  For I:=1 to Dimen do
    For J:=1 to Dimen do
        AddedMatrix[I][J] := InputMatrix1[I][J] + InputMatrix2[I][J];
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure MatrixSubtract         (    Dimen:             Integer;
                                  Var InputMatrix1:      Matrix;
                                  Var InputMatrix2:      Matrix;
                                  Var Matrix1Minus2:     Matrix);
Var I,J:Integer;
Begin
  For I:=1 to Dimen do
    For J:=1 to Dimen do
        Matrix1Minus2[I][J] := InputMatrix1[I][J] - InputMatrix2[I][J];
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure MatrixScalarMultiply   (    Dimen:             Integer;
                                  Var InputMatrix:       Matrix;
                                      Scalar:            Real;
                                  Var ProductMatrix:     Matrix);
Var I,J:Integer;
Begin
  For I:=1 to Dimen do
    For J:=1 to Dimen do
        ProductMatrix[I][J] := Scalar * InputMatrix[I][J];
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure BigMatrixScalarMultiply(    Dimen:             Integer;
                                  Var InputMatrix:       Matrix2;
                                      Scalar:            Real;
                                  Var ProductMatrix:     Matrix2);
Var I,J:Integer;
Begin
  For I:=1 to Dimen do
    For J:=1 to Dimen do
        ProductMatrix[I][J] := Scalar * InputMatrix[I][J];
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure BigPtrMatrixScalarMultiply(    Dimen:             Integer;
                                  Var InputMatrix:       Matrix2ptr;
                                      Scalar:            Real;
                                  Var ProductMatrix:     Matrix2ptr);
Var I,J:Integer;
Begin
  For I:=1 to Dimen do
    For J:=1 to Dimen do
        ProductMatrix^[I][J] := Scalar * InputMatrix^[I][J];
End;

{═════════════════════════════════════════════════════════════════════════════}

Procedure MatrixInverse          (    Dimen:             Integer;
                                  Var Input:             Matrix;
                                  Var Output:            Matrix;
                                  Var Error:             Byte);
{This program was copied from an old fortran routine of unknown origin.  It
 is better than the Turbo Numerical methods inversion routine in that it
 can handle larger arrays.}
Var
  Determ:       double;
  I,J,K,L:      Integer;
  Ipivot:       array[1..Arraysize]  of integer;
  pivot:        array[1..Arraysize]  of Real;
  Index2:       Array[1..Arraysize, 1..2] of integer;
  Irow:         Integer;
  Icol:         Integer;
  Swap:         Real;
  Amax:         Real;
  L1:           Integer;
  T:            Real;
  Jrow, Jcol:   Integer;

Begin
  { Initialization }
{ for i:=1 to dimen do                                     { debug }
{   for j:=1 to dimen do                                   { debug }
{     writeln('input[',i:2,',',j:2,'] = ', input[i][j]); { debug }
  MatrixScalarMultiply ( Dimen, Input, 1, Output); {copies input to output}
{ for i:=1 to dimen do                                     { debug }
{   for j:=1 to dimen do                                   { debug }
{     writeln('output[',i:2,',',j:2,'] = ', output[i][j]); { debug }
  Error := 0;
  Determ := 1.0;
  For I:=1 to Dimen do Ipivot[I] := 0;

  { Search for pivot element }
  For I:=1 to Dimen do
  Begin
    Amax := 0.0;
    icol := -1;
    For J:=1 to Dimen do
      If Ipivot[J] <> 1 then
      Begin
        For K:=1 to Dimen do
        Begin
          If Ipivot[K] -1 > 0 then exit;
          if (ipivot[K] -1 < 0) and ((abs(amax)-abs(Output[J][K]))<0) then
          Begin
            Irow := J;
            Icol := K;
            Amax := Output[J][K];
          End;
        End;
      End;
    if icol < 0 then
    begin
      error := 1;
      writeln('Zero Matrix?');
      exit;
    end;
    Ipivot[Icol] := Ipivot[Icol]+1;

    { Interchange rows to put pivot element on diagonal }
    If Irow <> Icol then
    Begin
      Determ := -Determ;
      For L:= 1 to Dimen do
      Begin
        Swap := Output[Irow][L];
        Output[Irow][L] := Output[Icol][L];
        Output[Icol][L] := Swap;
      End;
    End;
    Index2[I][1] := Irow;
    Index2[I][2] := Icol;
    Pivot[I] := Output[Icol][Icol];
{writeln('determ, pivot[',i,']: ', determ, pivot[i]); { debug }
    Determ := Determ * Pivot[I];

    { Divide pivot row by pivot element }
    Output[Icol][Icol] := 1.0;
    For L:=1 to Dimen do
      If pivot[I] <> 0.0 then
        Output[Icol][L] := Output[Icol][L] / Pivot[I]
      Else
      Begin
        Error := 1;
        Exit;
      End;

    { Reduce non-pivot rows }
    For L1 := 1 to Dimen do
      If l1 <> Icol then
      Begin
        T := Output[L1][Icol];
        Output[L1][Icol] := 0.0;
        For L:=1 to Dimen do
          Output[L1][L] := Output[L1][L] - Output[Icol][L]*T;
      End;
  End;

  { Interchange columns }
  For I:=1 to Dimen do
  Begin
    L := Dimen + 1 - I;
    If Index2[L,1] <> Index2[L,2] then
    Begin
      Jrow := Index2[L,1];
      Jcol := Index2[L,2];
      For K:=1 to dimen do
      Begin
        Swap := Output[K][Jrow];
        Output[K][Jrow] := Output[K][Jcol];
        Output[K][Jcol] := Swap;
      End;
    End;
  End;

End;


{═════════════════════════════════════════════════════════════════════════════}


Procedure BigMatrixInverse       (    Dimen:             Integer;
                                  Var Input:             Matrix2;
                                  Var Output:            Matrix2;
                                  Var Error:             Byte);
{This program was copied from an old fortran routine of unknown origin.  It
 is better than the Turbo Numerical methods inversion routine in that it
 can handle larger arrays.}
Var
  Determ:       Double;
  I,J,K,L:      Integer;
  Ipivot:       array[1..Arraysize2]  of integer;
  pivot:        array[1..Arraysize2]  of Real;
  Index2:       Array[1..Arraysize2, 1..2] of integer;
  Irow:         Integer;
  Icol:         Integer;
  Swap:         Real;
  Amax:         Real;
  L1:           Integer;
  T:            Real;
  Jrow, Jcol:   Integer;

Begin
  { Initialization }
{ for i:=1 to dimen do                                     { debug }
{   for j:=1 to dimen do                                   { debug }
{     writeln('input[',i:2,',',j:2,'] = ', input[i][j]); { debug }
  BigMatrixScalarMultiply ( Dimen, Input, 1, Output); {copies input to output}
{ for i:=1 to dimen do                                     { debug }
{   for j:=1 to dimen do                                   { debug }
{     writeln('output[',i:2,',',j:2,'] = ', output[i][j]); { debug }
  Error := 0;
  Determ := 1.0;
  For I:=1 to Dimen do Ipivot[I] := 0;

  { Search for pivot element }
  For I:=1 to Dimen do
  Begin
    Amax := 0.0;
    icol := -1;
    For J:=1 to Dimen do
      If Ipivot[J] <> 1 then
      Begin
        For K:=1 to Dimen do
        Begin
          If Ipivot[K] -1 > 0 then exit;
          if (ipivot[K] -1 < 0) and ((abs(amax)-abs(Output[J][K]))<0) then
          Begin
            Irow := J;
            Icol := K;
            Amax := Output[J][K];
          End;
        End;
      End;
    if icol < 0 then
    begin
      error := 1;
      writeln('Zero Matrix?');
      exit;
    end;
    Ipivot[Icol] := Ipivot[Icol]+1;

    { Interchange rows to put pivot element on diagonal }
    If Irow <> Icol then
    Begin
      Determ := -Determ;
      For L:= 1 to Dimen do
      Begin
        Swap := Output[Irow][L];
        Output[Irow][L] := Output[Icol][L];
        Output[Icol][L] := Swap;
      End;
    End;
    Index2[I][1] := Irow;
    Index2[I][2] := Icol;
    Pivot[I] := Output[Icol][Icol];
{writeln('determ, pivot[',i,']: ', determ, pivot[i]); { debug }
    Determ := Determ * Pivot[I];

    { Divide pivot row by pivot element }
    Output[Icol][Icol] := 1.0;
    For L:=1 to Dimen do
      If pivot[I] <> 0.0 then
        Output[Icol][L] := Output[Icol][L] / Pivot[I]
      Else
      Begin
        Error := 1;
        Exit;
      End;

    { Reduce non-pivot rows }
    For L1 := 1 to Dimen do
      If l1 <> Icol then
      Begin
        T := Output[L1][Icol];
        Output[L1][Icol] := 0.0;
        For L:=1 to Dimen do
          Output[L1][L] := Output[L1][L] - Output[Icol][L]*T;
      End;
  End;

  { Interchange columns }
  For I:=1 to Dimen do
  Begin
    L := Dimen + 1 - I;
    If Index2[L,1] <> Index2[L,2] then
    Begin
      Jrow := Index2[L,1];
      Jcol := Index2[L,2];
      For K:=1 to dimen do
      Begin
        Swap := Output[K][Jrow];
        Output[K][Jrow] := Output[K][Jcol];
        Output[K][Jcol] := Swap;
      End;
    End;
  End;

End;

{═════════════════════════════════════════════════════════════════════════════}


Procedure BigPtrMatrixInverse    (    Dimen:             Integer;
                                  Var Input:             Matrix2ptr;
                                  Var Output:            Matrix2ptr;
                                  Var Error:             Byte);
{This program was copied from an old fortran routine of unknown origin.  It
 is better than the Turbo Numerical methods inversion routine in that it
 can handle larger arrays.}
Var
  Determ:       Double;
  I,J,K,L:      Integer;
  Ipivot:       array[1..Arraysize2]  of integer;
  pivot:        array[1..Arraysize2]  of Real;
  Index2:       Array[1..Arraysize2, 1..2] of integer;
  Irow:         Integer;
  Icol:         Integer;
  Swap:         Real;
  Amax:         Real;
  L1:           Integer;
  T:            Real;
  Jrow, Jcol:   Integer;

Begin
  { Initialization }
{ for i:=1 to dimen do                                     { debug }
{   for j:=1 to dimen do                                   { debug }
{     writeln('input[',i:2,',',j:2,'] = ', input[i][j]); { debug }
  BigPtrMatrixScalarMultiply ( Dimen, Input, 1, Output); {copies input to output}
{ for i:=1 to dimen do                                     { debug }
{   for j:=1 to dimen do                                   { debug }
{     writeln('output[',i:2,',',j:2,'] = ', output[i][j]); { debug }
  Error := 0;
  Determ := 1.0;
  For I:=1 to Dimen do Ipivot[I] := 0;

  { Search for pivot element }
  For I:=1 to Dimen do
  Begin
    Amax := 0.0;
    icol := -1;
    For J:=1 to Dimen do
      If Ipivot[J] <> 1 then
      Begin
        For K:=1 to Dimen do
        Begin
          If Ipivot[K] -1 > 0 then exit;
          if (ipivot[K] -1 < 0) and ((abs(amax)-abs(Output^[J][K]))<0) then
          Begin
            Irow := J;
            Icol := K;
            Amax := Output^[J][K];
          End;
        End;
      End;
    if icol < 0 then
    begin
      error := 1;
      writeln('Zero Matrix?');
      exit;
    end;
    Ipivot[Icol] := Ipivot[Icol]+1;

    { Interchange rows to put pivot element on diagonal }
    If Irow <> Icol then
    Begin
      Determ := -Determ;
      For L:= 1 to Dimen do
      Begin
        Swap := Output^[Irow][L];
        Output^[Irow][L] := Output^[Icol][L];
        Output^[Icol][L] := Swap;
      End;
    End;
    Index2[I][1] := Irow;
    Index2[I][2] := Icol;
    Pivot[I] := Output^[Icol][Icol];
{writeln('determ, pivot[',i,']: ', determ, pivot[i]); { debug }
    Determ := Determ * Pivot[I];

    { Divide pivot row by pivot element }
    Output^[Icol][Icol] := 1.0;
    For L:=1 to Dimen do
      If pivot[I] <> 0.0 then
        Output^[Icol][L] := Output^[Icol][L] / Pivot[I]
      Else
      Begin
        Error := 1;
        Exit;
      End;

    { Reduce non-pivot rows }
    For L1 := 1 to Dimen do
      If l1 <> Icol then
      Begin
        T := Output^[L1][Icol];
        Output^[L1][Icol] := 0.0;
        For L:=1 to Dimen do
          Output^[L1][L] := Output^[L1][L] - Output^[Icol][L]*T;
      End;
  End;

  { Interchange columns }
  For I:=1 to Dimen do
  Begin
    L := Dimen + 1 - I;
    If Index2[L,1] <> Index2[L,2] then
    Begin
      Jrow := Index2[L,1];
      Jcol := Index2[L,2];
      For K:=1 to dimen do
      Begin
        Swap := Output^[K][Jrow];
        Output^[K][Jrow] := Output^[K][Jcol];
        Output^[K][Jcol] := Swap;
      End;
    End;
  End;

End;


{═════════════════════════════════════════════════════════════════════════════}

Procedure ComplexMatrixInverse   (    Dimen:             Integer;
                                  Var RealInputMatrix:   Matrix;
                                  Var ImagInputMatrix:   Matrix;
                                  Var RealInverseMatrix: Matrix;
                                  Var ImagInverseMatrix: Matrix;
                                  Var Error:             Byte);


Var
  Phase:    Complex;
  I,J,K:    Integer;
  big:      Matrix2ptr;
  bigout:   Matrix2ptr;

begin
  new(big);
  new(bigout);
  for i:=1 to Dimen do
    for j:=1 to Dimen do
    begin
      big^[i,j] := RealInputMatrix[i,j];
      big^[i+Dimen, j+Dimen] := RealInputMatrix[i,j];
    end;
  for i:=1 to Dimen do
    for j:=1 to Dimen do
    begin
      big^[i, j+Dimen] := ImagInputMatrix[i,j];
      big^[i+Dimen, j] := ImagInputMatrix[i,j];
    end;
    BigPtrMatrixInverse ( Dimen * 2, big, bigout, error);
    for i:=1 to Dimen do
      for j:=1 to Dimen do
      begin
        realinversematrix[i,j] :=  bigout^[i,j];
        imaginversematrix[i,j] := -bigout^[i,j+Dimen];
      end;
    dispose(bigout);
    dispose(big);
  end;

{═════════════════════════════════════════════════════════════════════════════}


Begin
  Randomize;
  odd := False;
End.