Browse code at this revision
Parent(s): [3178ec]
| added | pmath/uibeta.p |
| added | pmath/uexpdist.p |
| added | pmath/uibtdist.p |
| added | pmath/umath.p |
| added | pmath/upoidist.p |
| added | pmath/Makefile |
| added | pmath/ugamdist.p |
| added | pmath/uigmdist.p |
| added | pmath/uigamma.p |
| removed | CERNlib.def.m2cc |
| removed | gamit.f |
| copied | SpezFunkC.def.m2cc -> pmath/upolev.p |
| copied | SpezFunkF77.def.m2cc -> pmath/types.inc |
| copied | TestSpezFunkt1.mod -> pmath/ugamma.p |
| copied | TstSpFun4pas.p -> pmath/TstSpFun4pas.p |
| copied | gamit.def.m2cc -> pmath/uminmax.p |
| copied | libmath.def.m2cc -> pmath/utypes.p |
CERNlib.def.m2cc
File was removed.
gamit.f
File was removed.
SpezFunkC.def.m2cc to pmath/upolev.p
--- a/SpezFunkC.def.m2cc +++ b/pmath/upolev.p @@ -1,33 +1,64 @@ -#ifdef __XDS__ -DEFINITION MODULE ["C"] SpezFunkC; -#endif -#ifdef __GM2__ -DEFINITION MODULE FOR "C" SpezFunkC; (* GNU M2 *) -#endif - - (*------------------------------------------------------------------------*) - (* interface to some C based functions *) - (* *) - (* Letzte Veraenderung: *) - (* *) - (* 16.05.18, MRi: Erstellen der ersten Version *) - (*------------------------------------------------------------------------*) - - (* $Id: SpezFunkC.def.m2cc,v 1.1 2018/06/21 08:01:13 mriedl Exp $ *) - -#ifdef __XDS__ -TYPE DOUBLE = LONGREAL; -#endif -#ifdef __GM2__ -TYPE DOUBLE = REAL; -#endif - - -PROCEDURE betai(a : DOUBLE; - b : DOUBLE; - x : DOUBLE) : DOUBLE; - -PROCEDURE gammq(a : DOUBLE; - x : DOUBLE) : DOUBLE; - -END SpezFunkC. +{ ****************************************************************** + Polynomial evaluations for special functions. + Translated from C code in Cephes library (http://www.moshier.net) + ****************************************************************** } + +unit upolev; + +interface + +uses + utypes; + +type + TabCoef = array[0..9] of Float; + +function PolEvl(var X : Float; Coef : TabCoef; N : Integer) : Float; + +function P1Evl(var X : Float; Coef : TabCoef; N : Integer) : Float; + +implementation + + function PolEvl(var X : Float; Coef : TabCoef; N : Integer) : Float; +{ ------------------------------------------------------------------ + Evaluates polynomial of degree N: + + 2 N + y = C + C x + C x +...+ C x + 0 1 2 N + + Coefficients are stored in reverse order: + + Coef[0] = C , ..., Coef[N] = C + N 0 + + The function P1Evl() assumes that Coef[N] = 1.0 and is + omitted from the array. Its calling arguments are + otherwise the same as PolEvl(). + ------------------------------------------------------------------ } + var + Ans : Float; + I : Integer; + begin + Ans := Coef[0]; + for I := 1 to N do + Ans := Ans * X + Coef[I]; + PolEvl := Ans; + end; + + function P1Evl(var X : Float; Coef : TabCoef; N : Integer) : Float; +{ ------------------------------------------------------------------ + Evaluate polynomial when coefficient of X is 1.0. + Otherwise same as PolEvl. + ------------------------------------------------------------------ } + var + Ans : Float; + I : Integer; + begin + Ans := X + Coef[0]; + for I := 1 to N - 1 do + Ans := Ans * X + Coef[I]; + P1Evl := Ans; + end; + +end.
SpezFunkF77.def.m2cc to pmath/types.inc
--- a/SpezFunkF77.def.m2cc +++ b/pmath/types.inc @@ -1,166 +1,305 @@ -#ifdef __XDS__ -DEFINITION MODULE ["C"] SpezFunkF77; -#endif -#ifdef __GM2__ -DEFINITION MODULE FOR "C" libmath; (* GNU M2 *) -#endif - - (*------------------------------------------------------------------------*) - (* interface to CERN besj0|besj1 function *) - (* *) - (* Letzte Veraenderung: *) - (* *) - (* 13.02.18, MRi: Erstellen der ersten Version *) - (* 15.05.18, MRi: Ergaenzen weiterer Funktionen in den letzten Tagen ... *) - (*------------------------------------------------------------------------*) - (* [1] Luke, Y.L.; "Mathematical functions and their approximations", *) - (* Academic Press, New York (US), pp 322â324 (1975) *) - (*------------------------------------------------------------------------*) - - (* $Id: SpezFunkF77.def.m2cc,v 1.1 2018/06/21 08:08:19 mriedl Exp mriedl $ *) - -#ifdef __XDS__ -TYPE DOUBLEPRECISION = LONGREAL; - DOUBLECOMPLEX = LONGCOMPLEX; - INTEGER4 = INTEGER; -#endif -#ifdef __GM2__ -TYPE DOUBLEPRECISION = REAL; - DOUBLECOMPLEX = COMPLEX; - INTEGER4 = INTEGER; -#endif - -PROCEDURE besselj0_(VAR x : DOUBLEPRECISION) : DOUBLEPRECISION; - - (*----------------------------------------------------------------*) - (* Interface to Fortran 2003 intrinsic Bessel function *) - (*----------------------------------------------------------------*) - -PROCEDURE dbesj0_(VAR x : DOUBLEPRECISION) : DOUBLEPRECISION; -PROCEDURE dbesj1_(VAR x : DOUBLEPRECISION) : DOUBLEPRECISION; -PROCEDURE dbesy0_(VAR x : DOUBLEPRECISION) : DOUBLEPRECISION; -PROCEDURE dbesy1_(VAR x : DOUBLEPRECISION) : DOUBLEPRECISION; - - (*----------------------------------------------------------------*) - (* CERN DBESJ0(x) with ENTRY DBESJ1(x), DBESY0(x), DBESY1(x) *) - (*----------------------------------------------------------------*) - -PROCEDURE dbessj0_(VAR x : DOUBLEPRECISION) : DOUBLEPRECISION; - - (*----------------------------------------------------------------*) - (* Reworked CERN version, only Bessel J0 *) - (*----------------------------------------------------------------*) - -PROCEDURE bessj0_(VAR x : DOUBLEPRECISION) : DOUBLEPRECISION; - - (*----------------------------------------------------------------*) - (* Single precision expansion (Abrahmowitz/Stegun) *) - (*----------------------------------------------------------------*) - -PROCEDURE bsslj_(VAR a : DOUBLECOMPLEX; - VAR In : INTEGER4; - VAR w : DOUBLECOMPLEX); - - (*----------------------------------------------------------------*) - (* NSWC SUBROUTINE BSSLJ(A,In,W) *) - (*----------------------------------------------------------------*) - -PROCEDURE rbesy_(VAR dnu : DOUBLEPRECISION; - VAR n : INTEGER4; - VAR Vm : DOUBLEPRECISION; - VAR By : ARRAY OF DOUBLEPRECISION; - VAR ier : INTEGER4); - - (*----------------------------------------------------------------*) - (* THIS SUBROUTINE COMPUTES THE BESSEL FUNCTION 'Y' OF REAL *) - (* ORDER DNU, DNU+1, . . ., DNU+N AND REAL ARGUMENT 'X'. *) - (*----------------------------------------------------------------*) - -PROCEDURE besj1_(VAR x : DOUBLEPRECISION) : DOUBLEPRECISION; -PROCEDURE besy1_(VAR x : DOUBLEPRECISION) : DOUBLEPRECISION; - - (*----------------------------------------------------------------*) - (* Cody, W.J. Version, Argonne National Laboratory *) - (*----------------------------------------------------------------*) - -PROCEDURE hankelh01_(VAR a : DOUBLEPRECISION) : DOUBLECOMPLEX; - - - (*----------------------------------------------------------------*) - (* Interface to CH12N of Shanjie Zhang and Jianming Jin *) - (*----------------------------------------------------------------*) - -PROCEDURE chae01_(VAR Z : DOUBLEPRECISION; - VAR RelErr : DOUBLEPRECISION; - VAR H01 : DOUBLECOMPLEX); - - (*----------------------------------------------------------------*) - (* Hankels asymptotic expansion *) - (*----------------------------------------------------------------*) - -PROCEDURE normprb_(VAR Z : DOUBLEPRECISION; - VAR P : DOUBLEPRECISION; - VAR Q : DOUBLEPRECISION; - VAR PDF : DOUBLEPRECISION); - - (*----------------------------------------------------------------*) - (* Normal distribution probabilities accurate to 1.e-15. *) - (* Z = no. of standard deviations from the mean. *) - (* P, Q = probabilities to the left & right of Z. P + Q = 1. *) - (* PDF = the probability density. *) - (* *) - (* Based upon algorithm 5666 for the error function, from: *) - (* Hart, J.F. et al, 'Computer Approximations', Wiley 1968 *) - (* *) - (* Programmer: Alan Miller *) - (*----------------------------------------------------------------*) - -PROCEDURE dlgama_(VAR X : DOUBLEPRECISION) : DOUBLEPRECISION; - - (*----------------------------------------------------------------*) - (* This routine calculates the LOG(GAMMA) function for a positive *) - (* real argument X. *) - (*----------------------------------------------------------------*) -(* -PROCEDURE betain_(VAR x : DOUBLEPRECISION; - VAR p : DOUBLEPRECISION; - VAR q : DOUBLEPRECISION; - VAR beta : DOUBLEPRECISION; - VAR ifault : INTEGER4) : DOUBLEPRECISION; - *) - -PROCEDURE cgamma_(VAR z : DOUBLECOMPLEX) : DOUBLECOMPLEX; - - (*----------------------------------------------------------------*) - (* CERN library 2006 *) - (*----------------------------------------------------------------*) - -PROCEDURE clngamma_(VAR z : DOUBLECOMPLEX) : DOUBLECOMPLEX; - - (*----------------------------------------------------------------*) - (* CERN library 2006 *) - (*----------------------------------------------------------------*) - -PROCEDURE cdgamma_(VAR z : DOUBLECOMPLEX; - VAR l : INTEGER4) : DOUBLECOMPLEX; - - (*----------------------------------------------------------------*) - (* Complex gamma function in double precision *) - (*----------------------------------------------------------------*) - -PROCEDURE zdgamma_(VAR z : DOUBLECOMPLEX) : DOUBLECOMPLEX; - - (*----------------------------------------------------------------*) - (* Complex gamma function in double precision *) - (* *) - (* The source of the original subroutine is *) - (* *) - (* http://momonga.t.u-tokyo.ac.jp/-ooura/gamerf.html *) - (* *) - (* Takuya Ooura *) - (* Research Institute for Mathematical Sciences *) - (* Kyoto University *) - (* Kyoto 606-01 Japan *) - (*----------------------------------------------------------------*) - -END SpezFunkF77. +{ ------------------------------------------------------------------ + Floating point type (Default = Double) + ------------------------------------------------------------------ } + +{$IFDEF SINGLEREAL} + type Float = Single; +{$ELSE} +{$IFDEF EXTENDEDREAL} + type Float = Extended; +{$ELSE} + {$DEFINE DOUBLEREAL} + type Float = Double; +{$ENDIF} +{$ENDIF} + +{ ------------------------------------------------------------------ + Mathematical constants + ------------------------------------------------------------------ } + +const + Pi = 3.14159265358979323846; { Pi } + Ln2 = 0.69314718055994530942; { Ln(2) } + Ln10 = 2.30258509299404568402; { Ln(10) } + LnPi = 1.14472988584940017414; { Ln(Pi) } + InvLn2 = 1.44269504088896340736; { 1/Ln(2) } + InvLn10 = 0.43429448190325182765; { 1/Ln(10) } + TwoPi = 6.28318530717958647693; { 2*Pi } + PiDiv2 = 1.57079632679489661923; { Pi/2 } + SqrtPi = 1.77245385090551602730; { Sqrt(Pi) } + Sqrt2Pi = 2.50662827463100050242; { Sqrt(2*Pi) } + InvSqrt2Pi = 0.39894228040143267794; { 1/Sqrt(2*Pi) } + LnSqrt2Pi = 0.91893853320467274178; { Ln(Sqrt(2*Pi)) } + Ln2PiDiv2 = 0.91893853320467274178; { Ln(2*Pi)/2 } + Sqrt2 = 1.41421356237309504880; { Sqrt(2) } + Sqrt2Div2 = 0.70710678118654752440; { Sqrt(2)/2 } + Gold = 1.61803398874989484821; { Golden Mean = (1 + Sqrt(5))/2 } + CGold = 0.38196601125010515179; { 2 - GOLD } + +{ ------------------------------------------------------------------ + Machine-dependent constants + ------------------------------------------------------------------ } + +{$IFDEF SINGLEREAL} +const + MachEp = 1.192093E-7; { Floating point precision: 2^(-23) } + MaxNum = 3.402823E+38; { Max. floating point number: 2^128 } + MinNum = 1.175495E-38; { Min. floating point number: 2^(-126) } + MaxLog = 88.72283; { Max. argument for Exp = Ln(MaxNum) } + MinLog = -87.33655; { Min. argument for Exp = Ln(MinNum) } + MaxFac = 33; { Max. argument for Factorial } + MaxGam = 34.648; { Max. argument for Gamma } + MaxLgm = 1.0383E+36; { Max. argument for LnGamma } +{$ENDIF} + +{$IFDEF DOUBLEREAL} +const + MachEp = 2.220446049250313E-16; { 2^(-52) } + MaxNum = 1.797693134862315E+308; { 2^1024 } + MinNum = 2.225073858507202E-308; { 2^(-1022) } + MaxLog = 709.7827128933840; + MinLog = -708.3964185322641; + MaxFac = 170; + MaxGam = 171.624376956302; + MaxLgm = 2.556348E+305; +{$ENDIF} + +{$IFDEF EXTENDEDREAL} +const + MachEp = 1.08420217248550444E-19; { 2^(-63) } + MaxNum = 5.9486574767861588254E+4931; { 2^16383 } + MinNum = 6.7242062862241870125E-4932; { 2^(-16381) } + MaxLog = 11355.830259113584004; + MinLog = -11354.443964752464114; + + MaxFac = 1754; + MaxGam = 1755.455; + MaxLgm = 1.04848146839019521E+4928; +{$ENDIF} + +{ ------------------------------------------------------------------ + Error codes for mathematical functions + ------------------------------------------------------------------ } + +const + FOk = 0; { No error } + FDomain = - 1; { Argument domain error } + FSing = - 2; { Function singularity } + FOverflow = - 3; { Overflow range error } + FUnderflow = - 4; { Underflow range error } + FTLoss = - 5; { Total loss of precision } + FPLoss = - 6; { Partial loss of precision } + +{ ------------------------------------------------------------------ + Error codes for matrix computations + ------------------------------------------------------------------ } + +const + MatOk = 0; { No error } + MatNonConv = -1; { Non-convergence } + MatSing = -2; { Quasi-singular matrix } + MatErrDim = -3; { Non-compatible dimensions } + MatNotPD = -4; { Matrix not positive definite } + +{ ------------------------------------------------------------------ + Error codes for optimization and nonlinear equations + ------------------------------------------------------------------ } + +const + OptOk = 0; { No error } + OptNonConv = -1; { Non-convergence } + OptSing = -2; { Quasi-singular hessian matrix } + OptBigLambda = -5; { Too high Marquardt parameter } + +{ ------------------------------------------------------------------ + Error codes for nonlinear regression + ------------------------------------------------------------------ } + +const + NLMaxPar = -6; { Max. number of parameters exceeded } + NLNullPar = -7; { Initial parameter equal to zero } + +{ ------------------------------------------------------------------ + Complex numbers + ------------------------------------------------------------------ } + +type Complex = record + X, Y : Float; +end; + +{ ------------------------------------------------------------------ + Vectors and matrices. + ------------------------------------------------------------------ } + +const { Maximal array size } +{$IFDEF _16BIT} + MaxSize = 65536; { 64 kilobytes : 2^16 } +{$ELSE} + MaxSize = 2147483648; { 2 gigabytes : 2^31 } +{$ENDIF} + + FltSize = SizeOf(Float); + CompSize = SizeOf(Complex); + IntSize = SizeOf(Integer); + BoolSize = SizeOf(Boolean); + PtrSize = SizeOf(Pointer); + +{$IFDEF __GPC__} + StrSize = 255; +{$ELSE} + StrSize = SizeOf(String); +{$ENDIF} + + MAX_FLT = Trunc(MaxSize / FltSize) - 2; + MAX_COMP = Trunc(MaxSize / CompSize) - 2; + MAX_INT = Trunc(MaxSize / IntSize) - 2; + MAX_BOOL = Trunc(MaxSize / BoolSize) - 2; + MAX_STR = Trunc(MaxSize / StrSize) - 2; + MAX_VEC = Trunc(MaxSize / PtrSize) - 2; + +type + TVector = array[0..MAX_FLT] of Float; + TIntVector = array[0..MAX_INT] of Integer; + TCompVector = array[0..MAX_COMP] of Complex; + TBoolVector = array[0..MAX_BOOL] of Boolean; + TStrVector = array[0..MAX_STR] of String; + + PVector = ^TVector; + PIntVector = ^TIntVector; + PBoolVector = ^TBoolVector; + PCompVector = ^TCompVector; + PStrVector = ^TStrVector; + +type + TMatrix = array[0..MAX_VEC] of PVector; + TIntMatrix = array[0..MAX_VEC] of PIntVector; + TBoolMatrix = array[0..MAX_VEC] of PBoolVector; + TCompMatrix = array[0..MAX_VEC] of PCompVector; + TStrMatrix = array[0..MAX_VEC] of PStrVector; + + PMatrix = ^TMatrix; + PIntMatrix = ^TIntMatrix; + PBoolMatrix = ^TBoolMatrix; + PCompMatrix = ^TCompMatrix; + PStrMatrix = ^TStrMatrix; + +{ ------------------------------------------------------------------ + Functional types + ------------------------------------------------------------------ } + +{ Function of one variable } +type TFunc = function(X : Float) : Float; + +{ Function of several variables } +type TFuncNVar = function(X : PVector) : Float; + +{ Nonlinear equation system } +type TEquations = procedure(X, F : PVector); + +{ Differential equation system } +type TDiffEqs = procedure(X : Float; Y, Yp : PVector); + +{ Jacobian } +type TJacobian = procedure(X : PVector; D : PMatrix); + +{ Gradient } +type TGradient = procedure(X, G : PVector); + +{ Hessian and Gradient } +type THessGrad = procedure(X, G : PVector; H : PMatrix); + +{ ------------------------------------------------------------------ + Random number generators + ------------------------------------------------------------------ } + +type RNG_Type = + (RNG_MWC, { Multiply-With-Carry } + RNG_MT, { Mersenne Twister } + RNG_UVAG); { Universal Virtual Array Generator } + +type + RNG_IntType = + {$IFDEF _16BIT}LongInt{$ELSE}Integer{$ENDIF}; { 32-bit, signed } + + RNG_LongType = + {$IFDEF _16BIT}LongInt{$ELSE}Cardinal{$ENDIF}; { 32-bit, signed or unsigned } + +{ ------------------------------------------------------------------ + Statistics + ------------------------------------------------------------------ } + +type StatClass = record { Statistical class } + Inf : Float; { Lower bound } + Sup : Float; { Upper bound } + N : Integer; { Number of values } + F : Float; { Frequency } + D : Float; { Density } +end; + +const + MAX_CLS = 1000; { Max. number of statistical classes } + +type + TStatClassVector = array[0..MAX_CLS] of StatClass; + PStatClassVector = ^TStatClassVector; + +{ ------------------------------------------------------------------ + Curve fit + ------------------------------------------------------------------ } + +type + TRegMode = (OLS, WLS); { Regression mode } + +type + TRegTest = record { Test of regression } + Vr : Float; { Residual variance } + R2 : Float; { Coefficient of determination } + R2a : Float; { Adjusted coeff. of determination } + F : Float; { Variance ratio (explained/residual) } + Nu1, Nu2 : Integer; { Degrees of freedom } + end; + +{ Optimization algorithms for nonlinear regression } +type + TOptAlgo = ( + NL_MARQ, { Marquardt algorithm } + NL_SIMP, { Simplex algorithm } + NL_BFGS, { BFGS algorithm } + NL_SA, { Simulated annealing } + NL_GA); { Genetic algorithm } + +{ Regression function } +type + TRegFunc = function(X : Float; B : PVector) : Float; + +{ Procedure to compute the derivatives of the regression function + with respect to the regression parameters } +type + TDerivProc = procedure(X, Y : Float; B, D : PVector); + +{ Variable of the integrated Michaelis equation: + Time, Substrate conc., Enzyme conc. } +type + TMintVar = (Var_T, Var_S, Var_E); + +{ ------------------------------------------------------------------ + Graphics + ------------------------------------------------------------------ } + +type + Str30 = String[30]; + TScale = (LinScale, LogScale); + TGrid = (NoGrid, HorizGrid, VertiGrid, BothGrid); + +{ ------------------------------------------------------------------ + Math parser + ------------------------------------------------------------------ } + +const + MaxArg = 26; { Max number of arguments for a function } + +type + TArgC = 1..MaxArg; + +type + TWrapper = function(ArgC : TArgC; ArgV : PVector) : Float;
TestSpezFunkt1.mod to pmath/ugamma.p
--- a/TestSpezFunkt1.mod +++ b/pmath/ugamma.p @@ -1,412 +1,344 @@ -MODULE TestSpezFunkt1; - - (*------------------------------------------------------------------------*) - (* Modul um einige Prozeduren aus SpezFunkt1 zu testen. *) - (* *) - (* Letzte Veraenderung *) - (* *) - (* 03.10.16, MRi: Erstellend der ersten Version *) - (* 04.10.16, MRi: Ersten Version von CubicEq *) - (* 12.10.16, MRi: Einfuehren von Sort und SolveP3 *) - (*------------------------------------------------------------------------*) - - (* - * $Id: TestSpezFunkt1.mod,v 1.1 2018/06/21 08:08:38 mriedl Exp mriedl $ - *) - -FROM LongComplexMath IMPORT scalarMult,abs; - IMPORT LongMath,LMathLib; -FROM LongMath IMPORT pi,cos,arccos; -FROM LMathLib IMPORT sqrt3; -FROM SpezFunkt1 IMPORT QuadGl,KubGl,CubicEq; -FROM RandomLib IMPORT Zufall; - IMPORT TFormAus; - IMPORT TIO; - -CONST WLn = TIO.WrLn; - WSt = TIO.WrStr; - WrLR = TIO.WrLngReal; - -<* IF DEFINED(debug) THEN *> -PROCEDURE WrCmplx(x : LONGCOMPLEX); - -BEGIN - WrLR(RE(x),12,-4); WrLR(IM(x),12,-4); WSt("*i"); -END WrCmplx; - -PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr; - -<* END *> - - -PROCEDURE QuickSortCmplxRC(VAR A : ARRAY OF LONGCOMPLEX; - n : INTEGER); - - PROCEDURE sort(L, R: INTEGER); - - VAR i,j : INTEGER; - w : LONGCOMPLEX; - x : LONGREAL; - BEGIN - i := L; j := R; - x := abs(A[(L+R) DIV 2]); - REPEAT - WHILE (abs(A[i]) < x) DO INC(i); END; - WHILE (x < abs(A[j])) DO DEC(j); END; - IF (i <= j) THEN - w := A[i]; A[i] := A[j]; A[j] := w; INC(i); DEC(j); - END; - UNTIL (i > j); - IF (L < j) THEN sort(L,j); END; - IF (i < R) THEN sort(i,R); END; - END sort; -BEGIN - sort(0,n-1); -END QuickSortCmplxRC; - -PROCEDURE Sort(VAR X1,X2,X3 : LONGCOMPLEX); - - VAR X : ARRAY [1..3] OF LONGCOMPLEX; -BEGIN - X[1]:=X1; - X[2]:=X2; - X[3]:=X3; - QuickSortCmplxRC(X,3); - X1:=X[1]; - X2:=X[2]; - X3:=X[3]; -END Sort; - -PROCEDURE SolveP3( a,b,c,d : LONGREAL; - VAR X1,X2,X3 : LONGCOMPLEX); - - (* Berechnet die Wurzeln der kubischen Gleichung *) - (* a*x**3 + b*x**2 + c*x + d *) - - CONST eps = LMathLib.MachEps; - - VAR P,Q,Q2,P3,A,B,Phi : LONGREAL; - Cp,Px,a0,a2,a3,re,im : LONGREAL; -BEGIN - a0 := a; - a := b / a0; - b := c / a0; - c := d / a0; - - a2 := a*a; - P := (a2 - 3.0*b) / 9.0; - Q := (a*(2.0*a2 - 9.0*b) + 27.0*c) / 54.0; - Q2 := Q*Q; - P3 := P*P*P; - a3 := a / 3.0; - - IF (Q2 < P3) THEN - Cp := Q / LongMath.sqrt(P3); - IF( Cp < -1.0) THEN Cp:= -1.0; ELSIF( Cp > 1.0) THEN Cp:= 1.0; END; - Phi := arccos(Cp); - Px := -2.0*LongMath.sqrt(P); - X1 := CMPLX(Px*cos( Phi / 3.0) - a3, 0.0); - X2 := CMPLX(Px*cos((Phi + 2.0*pi) / 3.0) - a3, 0.0); - X3 := CMPLX(Px*cos((Phi - 2.0*pi) / 3.0) - a3, 0.0); - ELSE - A := - LongMath.power(ABS(Q) + LongMath.sqrt(Q2 - P3),1./3.); - IF (Q < 0.0) THEN A := -A; END; - IF (ABS(A) < LMathLib.Small) THEN B:=0.0; ELSE B := P / A; END; - X1 := CMPLX( (A + B) - a3, 0.0); - re := -0.5* (A + B) - a3; - im := 0.5*sqrt3*(A - B); - IF (ABS(im) < eps*ABS(re)) THEN (* Zwei identische reelle Wurzeln *) - X2:=CMPLX(re,0.0); - X3:=CMPLX(re,0.0); - ELSE (* Ein paar konjugiert komplexe Wurzeln *) - X2 := CMPLX(re, im); - X3 := CMPLX(re,-im); - END; - END; -END SolveP3; -(* -PROCEDURE KubGl( A,B,C,D : LONGREAL; - VAR X1,X2,X3 : LONGCOMPLEX); - -<* IF DEFINED(debug) THEN *> - PROCEDURE CubicRedF(P,Q : LONGREAL; - Y : LONGCOMPLEX) : LONGCOMPLEX; - BEGIN - RETURN Y*Y*Y + scalarMult(3.0*P,Y) + CMPLX(2.0*Q,0.0); - END CubicRedF; -<* END *> - - PROCEDURE Root3(x : LONGREAL) : LONGREAL; - - VAR sign : LONGREAL; - BEGIN - sign:=1.0; IF (x < 0.0) THEN sign:=-1.0; END; - RETURN sign*LongMath.exp(LongMath.ln(ABS(x))/3.0); - END Root3; - - CONST sqrt3 = 1.7320508075688772; - W1 = CMPLX(-0.5, 0.5*sqrt3); - W2 = CMPLX(-0.5,-0.5*sqrt3); - zero = CMPLX(0.,0.); - - VAR P,Q,Z,T,Phi,P3 : LONGREAL; - Z1,Z2,Cp,Px,A0,A3 : LONGREAL; - V,U : LONGCOMPLEX; - Y1,Y2,Y3 : LONGCOMPLEX; -BEGIN - A0 := A; - A := B / A0; - B := C / A0; - C := D / A0; - - P := - sqr(A)/9.0 + B/3.0; - Q := A*(sqr(A)/27.0 - B/6.0) + C/2.0; - - P3 := sqr(P)*P; - Z := sqr(Q) + P3; - A3 := A / 3.0; - - IF (ABS(Z) < 1.0E-15) THEN - T := Root3(-Q); - Y1 := CMPLX(2.0*T - 0.0, 0.0); - Y2 := CMPLX( - T - 0.0, 0.0); - Y3 := CMPLX( - T - 0.0, 0.0); - ELSIF (Z > 0.0) THEN - Z1 := - Q + LongMath.sqrt(Z); - Z2 := - Q - LongMath.sqrt(Z); - U := CMPLX(Root3(Z1),0.0); - V := CMPLX(Root3(Z2),0.0); - Y1 := U + V; - Y2 := W1*U + W2*V; - Y3 := W2*U + W1*V; - ELSIF (Z < 0.0) THEN - Cp := - Q / LongMath.sqrt(-P3); - Phi := LongMath.arccos(Cp); - Px :=-2.0*LongMath.sqrt(-P); - Y1 := CMPLX(-Px*LongMath.cos( Phi /3.0), 0.0); - Y2 := CMPLX( Px*LongMath.cos((Phi + LongMath.pi)/3.0), 0.0); - Y3 := CMPLX( Px*LongMath.cos((Phi - LongMath.pi)/3.0), 0.0); - ELSE - Y1:=zero; Y2:=zero; Y3:=zero; - END; - X1 := Y1 - CMPLX(A3,0.0); - X2 := Y2 - CMPLX(A3,0.0); - X3 := Y3 - CMPLX(A3,0.0); - <* IF DEFINED(debug) THEN *> - Sort(Y1,Y2,Y3); - WLn; - WSt(" Y1 = "); WrCmplx(Y1); WLn; - WSt(" Y2 = "); WrCmplx(Y2); WLn; - WSt(" Y3 = "); WrCmplx(Y3); WLn; - WLn; - WSt("CubicRedF(Y1) = "); WrCmplx(CubicRedF(P,Q,Y1)); WLn; - WSt("CubicRedF(Y2) = "); WrCmplx(CubicRedF(P,Q,Y2)); WLn; - WSt("CubicRedF(Y3) = "); WrCmplx(CubicRedF(P,Q,Y3)); WLn; - WLn; - <* END *> -END KubGl; -*) -PROCEDURE CubicF(A,B,C,D : LONGCOMPLEX; - X : LONGCOMPLEX) : LONGCOMPLEX; -BEGIN - RETURN A*X*X*X + B*X*X + C*X + D; -END CubicF; - -VAR A,B,C,D : LONGREAL; - -PROCEDURE F2x(x : LONGCOMPLEX) : LONGCOMPLEX; - -BEGIN - RETURN scalarMult(A,x*x) + scalarMult(B,x) + CMPLX(C,0.0); -END F2x; - -PROCEDURE F3x(x : LONGCOMPLEX) : LONGCOMPLEX; - -BEGIN - RETURN scalarMult(A,x*x*x) + scalarMult(B,x*x) + - scalarMult(C,x) + CMPLX(D,0.0); -END F3x; - -VAR X1,X2,X3 : LONGCOMPLEX; - r1,r2,r3 : LONGCOMPLEX; - Ac,Bc,Cc,Dc : LONGCOMPLEX; - X1w,X2w,X3w : LONGCOMPLEX; - dr1,dr2,dr3 : LONGREAL; - bsp : CARDINAL; -BEGIN - TIO.ExpAnz:=2; - - FOR bsp:=1 TO 200 DO - - TIO.WrCharN("=",78); TIO.WrLn; - - IF (bsp = 1) THEN -(* x(x+1)(x-1) = (x**2 + x)(x-1) = x**3 - x**2 + x**2 -x = x**3 - x*) - A := 1.0; - B := 0.0; - C := -1.0; - D := 0.0; - X1w := CMPLX( 1., 0.0); - X2w := CMPLX( 0., 0.0); - X3w := CMPLX(-1., 0.0); - ELSIF (bsp = 2) THEN -(* x(x+1)(x-1)+1 = (x**2 + x)(x-1) + 1 *) - A := 1.0; - B := 0.0; - C := -1.0; - D := 1.0; - X1w := CMPLX(0., 0.0); - X2w := CMPLX(0., 0.0); - X3w := CMPLX(LongMath.power(2.,1./3.)-LongMath.power(4.,1./3.)-1.0,0.0); - - ELSIF (bsp = 3) THEN -(* x(x+1)(x-1)+(2/9)sqrt(3) *) - A := 1.0; - B := 0.0; - C := -1.0; - D := 2.0/9.0*LongMath.sqrt(3.0); - - ELSIF (bsp = 4) THEN - - A := 1.0; - B := 0.0; - C := -1.0; - D := 1.0; - X1w := CMPLX(0., 0.0); - X2w := CMPLX(0., 0.0); - X3w := CMPLX(LongMath.power(2.,1./3.)-LongMath.power(4.,1./3.)-1.0,0.0); - - ELSIF (bsp = 5) THEN - - A := 1.0; - B := 0.0; - C := 1.0; - D := 0.0; - X1w := CMPLX(0.0, 0.0); - X2w := CMPLX(0.0, 1.0); - X3w := CMPLX(0.0,-1.0); - - ELSIF (bsp = 6) THEN - - A := 3.0; - B := -8.0; - C := -11.0; - D := 10.0; - X1w := CMPLX(1.0 + LongMath.sqrt(6.), 0.0); - X2w := CMPLX(1.0 - LongMath.sqrt(6.), 0.0); - X3w := CMPLX(2.0 / 3.0 , 0.0); - - ELSIF (bsp >= 7) THEN - - WLn; WSt("Zufallszahlen ... "); WLn; - - A := 20.0*(0.5 - Zufall()); - B := 20.0*(0.5 - Zufall()); - C := 20.0*(0.5 - Zufall()); - D := 20.0*(0.5 - Zufall()); - - END; - - - TFormAus.Write5("/,X1,F7.3,S,F7.3,S,F7.3,/",A,"*x**2 + ",B,"*x + ",C); - - QuadGl(A,B,C,X1,X2); - TFormAus.Write4("/,X1,S,F14.7,F14.7,S,/","Wurzel 1 = ",RE(X1),IM(X1),"*i"); - TFormAus.Write4("X1,S,F14.7,F14.7,S,//","Wurzel 2 = ",RE(X2),IM(X2),"*i"); - r1 := F2x(X1); - r2 := F2x(X2); - TFormAus.Write4("X1,S,F14.7,F14.7,S,/", "F2x(W1) = ",RE(r1),IM(r1),"*i"); - TFormAus.Write4("X1,S,F14.7,F14.7,S,/", "F2x(W2) = ",RE(r2),IM(r2),"*i"); - - WLn; TIO.WrCharN("-",78); WLn; - - TFormAus.Write4("/,X1,F7.3,S,F7.3,S",A,"x**3 + ",B,"*x**2 + "); - TFormAus.Write3("F7.3,S,F7.3,//",C,"*x + ",D); - - KubGl(A,B,C,D,X1,X2,X3); - - Sort(X1,X2,X3); - TFormAus.Write4("X1,S,F14.7,F14.7,S,/", "Wurzel 1 = ",RE(X1),IM(X1),"*i"); - TFormAus.Write4("X1,S,F14.7,F14.7,S,/", "Wurzel 2 = ",RE(X2),IM(X2),"*i"); - TFormAus.Write4("X1,S,F14.7,F14.7,S,//","Wurzel 3 = ",RE(X3),IM(X3),"*i"); - - Ac:=CMPLX( A, 0.0); - Bc:=CMPLX( B, 0.0); - Cc:=CMPLX( C, 0.0); - Dc:=CMPLX( D, 0.0); - - CubicEq(Ac,Bc,Cc,Dc,X1w,X2w,X3w); - Sort(X1w,X2w,X3w); - - WLn; - WSt("X1w = "); WrLR(RE(X1w),12,7); WrLR(IM(X1w),12,7); WSt("*i"); WLn; - WSt("X2w = "); WrLR(RE(X2w),12,7); WrLR(IM(X2w),12,7); WSt("*i"); WLn; - WSt("X3w = "); WrLR(RE(X3w),12,7); WrLR(IM(X3w),12,7); WSt("*i"); WLn; - WLn; - - dr1 := RE(X1w) - RE(X1); - dr2 := RE(X2w) - RE(X2); - dr3 := RE(X3w) - RE(X3); - - WLn; - WSt("dr1 = "); WrLR(dr1,12,6); WLn; - WSt("dr2 = "); WrLR(dr2,12,6); WLn; - WSt("dr3 = "); WrLR(dr3,12,6); WLn; - WLn; - - r1 := F3x(X1); - r2 := F3x(X2); - r3 := F3x(X3); - TFormAus.Write4("X1,S,F14.7,F14.7,S,/", "F3x(W1) = ",RE(r1),IM(r1),"*i"); - TFormAus.Write4("X1,S,F14.7,F14.7,S,/", "F3x(W2) = ",RE(r2),IM(r2),"*i"); - TFormAus.Write4("X1,S,F14.7,F14.7,S,//","F3x(W3) = ",RE(r3),IM(r3),"*i"); - - IF (abs(r1 + r2 + r3) < 1.0E-12) THEN - WSt("Loesung richtig"); WLn; - ELSE - WSt("Loesung falsch !!!"); WLn; - END; - - CubicEq(Ac,Bc,Cc,Dc,X1,X2,X3); - Sort(X1,X2,X3); - - WLn; - WSt("X1 = "); WrLR(RE(X1),12,7); WrLR(IM(X1),12,7); WSt("*i"); WLn; - WSt("X2 = "); WrLR(RE(X2),12,7); WrLR(IM(X2),12,7); WSt("*i"); WLn; - WSt("X3 = "); WrLR(RE(X3),12,7); WrLR(IM(X3),12,7); WSt("*i"); WLn; - WLn; - - r1 := CubicF(Ac,Bc,Cc,Dc,X1); - r2 := CubicF(Ac,Bc,Cc,Dc,X2); - r3 := CubicF(Ac,Bc,Cc,Dc,X3); - - WLn; - WSt("r1 = "); WrLR(RE(r1),12,-3); WrLR(IM(r1),12,-3); WSt("*i"); WLn; - WSt("r2 = "); WrLR(RE(r2),12,-3); WrLR(IM(r2),12,-3); WSt("*i"); WLn; - WSt("r3 = "); WrLR(RE(r3),12,-3); WrLR(IM(r3),12,-3); WSt("*i"); WLn; - WLn; - - SolveP3(A,B,C,D,X1,X2,X3); - Sort(X1,X2,X3); - - WLn; - WSt("X1 = "); WrLR(RE(X1),12,7); WrLR(IM(X1),12,7); WSt("*i"); WLn; - WSt("X2 = "); WrLR(RE(X2),12,7); WrLR(IM(X2),12,7); WSt("*i"); WLn; - WSt("X3 = "); WrLR(RE(X3),12,7); WrLR(IM(X3),12,7); WSt("*i"); WLn; - WLn; - - r1 := F3x(X1); - r2 := F3x(X2); - r3 := F3x(X3); - TFormAus.Write4("X1,S,F14.7,F14.7,S,/", "F3x(X1) = ",RE(r1),IM(r1),"*i"); - TFormAus.Write4("X1,S,F14.7,F14.7,S,/", "F3x(X2) = ",RE(r2),IM(r2),"*i"); - TFormAus.Write4("X1,S,F14.7,F14.7,S,//","F3x(X3) = ",RE(r3),IM(r3),"*i"); - - IF (abs(r1 + r2 + r3) < 1.0E-12) THEN - WSt("Loesung richtig"); WLn; - ELSE - WSt("Loesung falsch !!!"); WLn; - END; - END; - -END TestSpezFunkt1. +{ ****************************************************************** + Gamma function and related functions. + Translated from C code in Cephes library (http://www.moshier.net) + ****************************************************************** } + +unit ugamma; + +interface + +uses + utypes, upolev; + +function SgnGamma(X : Float) : Integer; +{ Sign of Gamma function } + +function Stirling(X : Float) : Float; +{ Stirling's formula for the Gamma function } + +function StirLog(X : Float) : Float; +{ Approximate Ln(Gamma) by Stirling's formula, for X >= 13 } + +function Gamma(X : Float) : Float; +{ Gamma function } + +function LnGamma(X : Float) : Float; +{ Logarithm of Gamma function } + +implementation + + function SgnGamma(X : Float) : Integer; + begin + if X > 0.0 then + SgnGamma := 1 + else if Odd(Trunc(Abs(X))) then + SgnGamma := 1 + else + SgnGamma := - 1; + end; + + function Stirling(X : Float) : Float; + { Stirling's formula for the gamma function + Gamma(x) = Sqrt(2*Pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) + where P(x) is a polynomial } + const + STIR : TabCoef = ( + 7.147391378143610789273E-4, + - 2.363848809501759061727E-5, + - 5.950237554056330156018E-4, + 6.989332260623193171870E-5, + 7.840334842744753003862E-4, + - 2.294719747873185405699E-4, + - 2.681327161876304418288E-3, + 3.472222222230075327854E-3, + 8.333333333333331800504E-2, + 0); + + var + W, P : Float; + begin + W := 1.0 / X; + if X > 1024.0 then + begin + P := 6.97281375836585777429E-5 * W + 7.84039221720066627474E-4; + P := P * W - 2.29472093621399176955E-4; + P := P * W - 2.68132716049382716049E-3; + P := P * W + 3.47222222222222222222E-3; + P := P * W + 8.33333333333333333333E-2; + end + else + P := PolEvl(W, STIR, 8); + Stirling := Sqrt2Pi * Exp((X - 0.5) * Ln(X) - X) * (1.0 + W * P); + end; + + function GamSmall(X1, Z : Float) : Float; + { Gamma function for small values of the argument } + const + S : TabCoef = ( + - 1.193945051381510095614E-3, + 7.220599478036909672331E-3, + - 9.622023360406271645744E-3, + - 4.219773360705915470089E-2, + 1.665386113720805206758E-1, + - 4.200263503403344054473E-2, + - 6.558780715202540684668E-1, + 5.772156649015328608253E-1, + 1.000000000000000000000E0, + 0); + + SN : TabCoef = ( + 1.133374167243894382010E-3, + 7.220837261893170325704E-3, + 9.621911155035976733706E-3, + - 4.219773343731191721664E-2, + - 1.665386113944413519335E-1, + - 4.200263503402112910504E-2, + 6.558780715202536547116E-1, + 5.772156649015328608727E-1, + - 1.000000000000000000000E0, + 0); + + var + P : Float; + begin + if X1 = 0.0 then + begin + GamSmall := DefaultVal(FSing, MaxNum); + Exit; + end; + if X1 < 0.0 then + begin + X1 := - X1; + P := PolEvl(X1, SN, 8); + end + else + P := PolEvl(X1, S, 8); + GamSmall := Z / (X1 * P); + end; + + function StirLog(X : Float) : Float; + { Approximate Ln(Gamma) by Stirling's formula, for X >= 13 } + const + P : TabCoef = ( + 4.885026142432270781165E-3, + - 1.880801938119376907179E-3, + 8.412723297322498080632E-4, + - 5.952345851765688514613E-4, + 7.936507795855070755671E-4, + - 2.777777777750349603440E-3, + 8.333333333333331447505E-2, + 0, 0, 0); + + var + Q, W : Float; + begin + Q := Ln(X) * (X - 0.5) - X; + Q := Q + LnSqrt2Pi; + if X > 1.0E+10 then + StirLog := Q + else + begin + W := 1.0 / Sqr(X); + StirLog := Q + PolEvl(W, P, 6) / X; + end; + end; + + function Gamma(X : Float) : Float; + const + P : TabCoef = ( + 4.212760487471622013093E-5, + 4.542931960608009155600E-4, + 4.092666828394035500949E-3, + 2.385363243461108252554E-2, + 1.113062816019361559013E-1, + 3.629515436640239168939E-1, + 8.378004301573126728826E-1, + 1.000000000000000000009E0, + 0, 0); + + Q : TabCoef = ( + - 1.397148517476170440917E-5, + 2.346584059160635244282E-4, + - 1.237799246653152231188E-3, + - 7.955933682494738320586E-4, + 2.773706565840072979165E-2, + - 4.633887671244534213831E-2, + - 2.243510905670329164562E-1, + 4.150160950588455434583E-1, + 9.999999999999999999908E-1, + 0); + + var + SgnGam, N : Integer; + A, X1, Z : Float; + begin + SetErrCode(FOk); + SgnGam := SgnGamma(X); + + if (X = 0.0) or ((X < 0.0) and (Frac(X) = 0.0)) then + begin + Gamma := DefaultVal(FSing, SgnGam * MaxNum); + Exit; + end; + + if X > MaxGam then + begin + Gamma := DefaultVal(FOverflow, MaxNum); + Exit; + end; + + A := Abs(X); + if A > 13.0 then + begin + if X < 0.0 then + begin + N := Trunc(A); + Z := A - N; + if Z > 0.5 then + begin + N := N + 1; + Z := A - N; + end; + Z := Abs(A * Sin(Pi * Z)) * Stirling(A); + if Z <= Pi / MaxNum then + begin + Gamma := DefaultVal(FOverflow, SgnGam * MaxNum); + Exit; + end; + Z := PI / Z; + end + else + Z := Stirling(X); + Gamma := SgnGam * Z; + end + else + begin + Z := 1.0; + X1 := X; + while X1 >= 3.0 do + begin + X1 := X1 - 1.0; + Z := Z * X1; + end; + while X1 < - 0.03125 do + begin + Z := Z / X1; + X1 := X1 + 1.0; + end; + if X1 <= 0.03125 then + Gamma := GamSmall(X1, Z) + else + begin + while X1 < 2.0 do + begin + Z := Z / X1; + X1 := X1 + 1.0; + end; + if (X1 = 2.0) or (X1 = 3.0) then + Gamma := Z + else + begin + X1 := X1 - 2.0; + Gamma := Z * PolEvl(X1, P, 7) / PolEvl(X1, Q, 8); + end; + end; + end; + end; + + function LnGamma(X : Float) : Float; + const + P : TabCoef = ( + - 2.163690827643812857640E3, + - 8.723871522843511459790E4, + - 1.104326814691464261197E6, + - 6.111225012005214299996E6, + - 1.625568062543700591014E7, + - 2.003937418103815175475E7, + - 8.875666783650703802159E6, + 0, 0, 0); + + Q : TabCoef = ( + - 5.139481484435370143617E2, + - 3.403570840534304670537E4, + - 6.227441164066219501697E5, + - 4.814940379411882186630E6, + - 1.785433287045078156959E7, + - 3.138646407656182662088E7, + - 2.099336717757895876142E7, + 0, 0, 0); + + var + N : Integer; + A, X1, Z : Float; + begin + SetErrCode(FOk); + + if (X = 0.0) or ((X < 0.0) and (Frac(X) = 0.0)) then + begin + LnGamma := DefaultVal(FSing, MaxNum); + Exit; + end; + + if X > MaxLgm then + begin + LnGamma := DefaultVal(FOverflow, MaxNum); + Exit; + end; + + A := Abs(X); + if A > 34.0 then + begin + if X < 0.0 then + begin + N := Trunc(A); + Z := A - N; + if Z > 0.5 then + begin + N := N + 1; + Z := N - A; + end; + Z := A * Sin(Pi * Z); + if Z = 0.0 then + begin + LnGamma := DefaultVal(FOverflow, MaxNum); + Exit; + end; + Z := LnPi - Ln(Z) - StirLog(A); + end + else + Z := StirLog(X); + LnGamma := Z; + end + else if X < 13.0 then + begin + Z := 1.0; + X1 := X; + while X1 >= 3 do + begin + X1 := X1 - 1.0; + Z := Z * X1; + end; + while X1 < 2.0 do + begin + if Abs(X1) <= 0.03125 then + begin + LnGamma := Ln(Abs(GamSmall(X1, Z))); + Exit; + end; + Z := Z / X1; + X1 := X1 + 1.0; + end; + if Z < 0.0 then Z := - Z; + if X1 = 2.0 then + LnGamma := Ln(Z) + else + begin + X1 := X1 - 2.0; + LnGamma := X1 * PolEvl(X1, P, 6) / P1Evl(X1, Q, 7) + Ln(Z); + end; + end + else + LnGamma := StirLog(X); + end; + +end.
TstSpFun4pas.p to pmath/TstSpFun4pas.p
File was renamed.
gamit.def.m2cc to pmath/uminmax.p
--- a/gamit.def.m2cc +++ b/pmath/uminmax.p @@ -1,59 +1,98 @@ -#ifdef __XDS__ -DEFINITION MODULE ["C"] gamit; (* XDS *) -#endif -#ifdef __GM2__ -DEFINITION MODULE FOR "C" gamit; (* GNU M2 *) -#endif - - (*-----------------------------------------------------------------------*) - (* Schnittstelle zu gamit.f (wird in BoysInt2Lib genutzt) *) - (* *) - (* 13.02.17, MRi: Erstellen der ersten Version *) - (*-----------------------------------------------------------------------*) - - (* $Id *) - -#ifdef __XDS__ -TYPE REAL8 = LONGREAL; -#endif -#ifdef __GM2__ -TYPE REAL8 = REAL; -#endif - -PROCEDURE dgamit_(VAR a,x : REAL8) : REAL8; - - (*---------------------------------------------------------------*) - (* Calculate tricomis form of the incomplete gamma function *) - (*---------------------------------------------------------------*) - -PROCEDURE d9gmit_(VAR a,x,algap1,sgngam,alx : REAL8) : REAL8; - - (*---------------------------------------------------------------*) - (* compute tricomi's incomplete gamma function for small x *) - (*---------------------------------------------------------------*) - -PROCEDURE dgamr_(VAR x : REAL8) : REAL8; - - (*---------------------------------------------------------------*) - (* calculates the double precision reciprocal of the complete *) - (* gamma function for double precision argument x *) - (*---------------------------------------------------------------*) - -PROCEDURE dlngam_(VAR x : REAL8) : REAL8; - - (*---------------------------------------------------------------*) - (* calculates the double precision logarithm of the absolute *) - (* value of the gamma function for double precision argument x *) - (*---------------------------------------------------------------*) - -PROCEDURE dlgams_(VAR x,dlgam,sgngam : REAL8); - - (*---------------------------------------------------------------*) - (* dlgams(x,dlgam,sgngam) calculates the double precision *) - (* natural logarithm of the absolute value of the gamma function *) - (* for double precision argument x and stores the result in *) - (* double precision argument dlgam where sgngam defines the sign *) - (* of dlgam. *) - (*---------------------------------------------------------------*) - -END gamit. +{ ****************************************************************** + Minimum, maximum, sign and exchange + ****************************************************************** } + +unit uminmax; + +interface + +uses + utypes; + +function FMin(X, Y : Float) : Float; { Minimum of 2 reals } +function FMax(X, Y : Float) : Float; { Maximum of 2 reals } +function IMin(X, Y : Integer) : Integer; { Minimum of 2 integers } +function IMax(X, Y : Integer) : Integer; { Maximum of 2 integers } +function Sgn(X : Float) : Integer; { Sign (returns 1 if X = 0) } +function Sgn0(X : Float) : Integer; { Sign (returns 0 if X = 0) } +function DSgn(A, B : Float) : Float; { Sgn(B) * |A| } + +procedure FSwap(var X, Y : Float); { Exchange 2 reals } +procedure ISwap(var X, Y : Integer); { Exchange 2 integers } + +implementation + + function FMin(X, Y : Float) : Float; + begin + if X <= Y then + FMin := X + else + FMin := Y; + end; + + function FMax(X, Y : Float) : Float; + begin + if X >= Y then + FMax := X + else + FMax := Y; + end; + + function IMin(X, Y : Integer) : Integer; + begin + if X <= Y then + IMin := X + else + IMin := Y; + end; + + function IMax(X, Y : Integer) : Integer; + begin + if X >= Y then + IMax := X + else + IMax := Y; + end; + + function Sgn(X : Float) : Integer; + begin + if X >= 0.0 then + Sgn := 1 + else + Sgn := - 1; + end; + + function Sgn0(X : Float) : Integer; + begin + if X > 0.0 then + Sgn0 := 1 + else if X = 0.0 then + Sgn0 := 0 + else + Sgn0 := - 1; + end; + + function DSgn(A, B : Float) : Float; + begin + if B < 0.0 then DSgn := - Abs(A) else DSgn := Abs(A) + end; + + procedure FSwap(var X, Y : Float); + var + Temp : Float; + begin + Temp := X; + X := Y; + Y := Temp; + end; + + procedure ISwap(var X, Y : Integer); + var + Temp : Integer; + begin + Temp := X; + X := Y; + Y := Temp; + end; + +end.
libmath.def.m2cc to pmath/utypes.p
--- a/libmath.def.m2cc +++ b/pmath/utypes.p @@ -1,311 +1,541 @@ -#ifdef __XDS__ -DEFINITION MODULE ["C"] libmath; (* XDS *) -#endif -#ifdef __GM2__ -DEFINITION MODULE FOR "C" libmath; (* GNU M2 *) -#endif - - (*-----------------------------------------------------------------------*) - (* Interface to (some) C libmath functions and procedures *) - (*-----------------------------------------------------------------------*) - (* Letzte Veraenderung: *) - (* *) - (* 13.02.18, MRi: Erstellen der ersten Version *) - (* 13.05.18, MRi: Ergaenzen der Bessel-Funktionen *) - (* 16.05.18, MRi: Weitere Ergaenzungen *) - (*-----------------------------------------------------------------------*) - - (* $Id: libmath.def.m2cc,v 1.2 2018/06/21 08:08:19 mriedl Exp mriedl $ *) - -IMPORT SYSTEM; - -#ifdef __XDS__ -TYPE DOUBLE = LONGREAL; - INT8 = SYSTEM.INT8; - INTEGER32 = SYSTEM.int; -#endif -#ifdef __GM2__ -TYPE DOUBLE = REAL; - INT8 = SYSTEM.INTEGER8; - INTEGER32 = SYSTEM.INTEGER32; -#endif - -TYPE exception = RECORD - type : INTEGER32; - name : POINTER TO INT8; - arg1 : DOUBLE; - arg2 : DOUBLE; - retval : DOUBLE; - END; - -PROCEDURE matherr(VAR exept : exception) : INTEGER32; - - (*----------------------------------------------------------------*) - (* See man matherr(3) *) - (*----------------------------------------------------------------*) - -PROCEDURE isnan(x : DOUBLE) : INTEGER32; - - (*----------------------------------------------------------------*) - (* returns a nonzero value if (fpclassify(x) == FP_NAN *) - (*----------------------------------------------------------------*) - -PROCEDURE isinf(x : DOUBLE) : INTEGER32; - - (*----------------------------------------------------------------*) - (* returns 1 if x is positive and -1 if x is negative infinity. *) - (*----------------------------------------------------------------*) - -(*==========================================================================*) - -PROCEDURE acos(x : DOUBLE) : DOUBLE; - -PROCEDURE asin(x : DOUBLE) : DOUBLE; - -PROCEDURE atan(x : DOUBLE) : DOUBLE; - -PROCEDURE atan2(x : DOUBLE; - y : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The atan2 function calculates the principal value of the *) - (* arc tangent of y/x, using the signs of the two arguments to *) - (* determine the quadrant of the result. *) - (*----------------------------------------------------------------*) - -PROCEDURE cos(x : DOUBLE) : DOUBLE; - -PROCEDURE sin(x : DOUBLE) : DOUBLE; - -PROCEDURE tan(x : DOUBLE) : DOUBLE; - -PROCEDURE cosh(x : DOUBLE) : DOUBLE; - -PROCEDURE sinh(x : DOUBLE) : DOUBLE; - -PROCEDURE tanh(x : DOUBLE) : DOUBLE; - -PROCEDURE acosh(x : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The acosh function calculates the inverse hyperbolic cosine *) - (* of x, that is the value whose hyperbolic cosine is x. *) - (*----------------------------------------------------------------*) - -PROCEDURE asinh(x : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The asinh function calculates the inverse hyperbolic sine *) - (* of x, that is the value whose hyperbolic sine is x. *) - (*----------------------------------------------------------------*) - -PROCEDURE atanh(x : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The atanh function calculates the inverse hyperbolic tangent *) - (* of x, that is the value whose hyperbolic tangent is x. *) - (*----------------------------------------------------------------*) - -PROCEDURE log(x : DOUBLE) : DOUBLE; - -PROCEDURE log10(x : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The log10() function returns the base 10 logarithm of x. *) - (*----------------------------------------------------------------*) - -PROCEDURE log1p(x : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* log1p(x) returns a value equivalent to log (1 + x) *) - (* *) - (* It is computed in a way that is accurate even if the value *) - (* of x is near zero. *) - (*----------------------------------------------------------------*) - -PROCEDURE exp(x : DOUBLE) : DOUBLE; - -PROCEDURE expm1(x : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* expm1(x) returns a value equivalent to exp(x) - 1 *) - (* *) - (* It is computed in a way that is accurate even if the value *) - (* of x is near zero���a case where exp(x) - 1 would be *) - (* inaccurate due to subtraction of two numbers that are nearly *) - (* equal. *) - (*----------------------------------------------------------------*) - -PROCEDURE pow(x : DOUBLE; - y : DOUBLE) : DOUBLE; - -PROCEDURE sqrt(x : DOUBLE) : DOUBLE; - -PROCEDURE cbrt(x : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The cbrt function returns the (real) cube root of x. This *) - (* function cannot fail, every representable real value has a *) - (* representable real cube root. *) - (*----------------------------------------------------------------*) - -PROCEDURE fabs(x : DOUBLE) : DOUBLE; - -PROCEDURE hypot(x : DOUBLE; - y : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The hypot function returns sqrt(x*x+y*y). This is the length *) - (* of the hypotenuse of a right-angled triangle with sides of *) - (* length x and y, or the distance of the point (x,y) from the *) - (* origin. *) - (* The calculation is performed without undue overflow or *) - (* underflow during the intermediate steps of the calculation. *) - (*----------------------------------------------------------------*) - -PROCEDURE copysign(x : DOUBLE; - y : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The copysign functions return a value whose absolute value *) - (* matches that of x, but whose sign bit matches that of y. *) - (*----------------------------------------------------------------*) - -(*==========================================================================*) - -PROCEDURE erf(x : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The erf function returns the error function of x, defined as *) - (* erf(x) = 2/sqrt(pi)* integral from 0 to x of exp(-t*t) dt *) - (*----------------------------------------------------------------*) - -PROCEDURE erfc(x : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The erfc function returns the complementary error function *) - (* of x, that is, 1.0 - erf(x). *) - (*----------------------------------------------------------------*) - -VAR signgam : INTEGER32; - -PROCEDURE lgamma(x: DOUBLE): DOUBLE; - - (*---------------------------------------------------------------*) - (* Returns the natural logarithm of the absolute value of the *) - (* Gamma function. The sign of lgamma is provided in the global *) - (* variable signgam. *) - (*---------------------------------------------------------------*) - -PROCEDURE tgamma(x: DOUBLE): DOUBLE; - - (*---------------------------------------------------------------*) - (* Returns Gamma(x) *) - (*---------------------------------------------------------------*) - -(*==========================================================================*) - -PROCEDURE j0(x : DOUBLE) : DOUBLE; - - (*---------------------------------------------------------------*) - (* Bessel functions of x of the first kind of order 0. *) - (*---------------------------------------------------------------*) - -PROCEDURE j1(x : DOUBLE) : DOUBLE; - - (*---------------------------------------------------------------*) - (* Bessel functions of x of the first kind of order 1. *) - (*---------------------------------------------------------------*) - -PROCEDURE jn(n : INTEGER; - x : DOUBLE) : DOUBLE; - - (*---------------------------------------------------------------*) - (* Bessel functions of x of the first kind of order n. *) - (*---------------------------------------------------------------*) - -PROCEDURE y0(x : DOUBLE) : DOUBLE; - - (*---------------------------------------------------------------*) - (* Bessel functions of x of the second kind of order 0. *) - (*---------------------------------------------------------------*) - -PROCEDURE y1(x : DOUBLE) : DOUBLE; - - (*---------------------------------------------------------------*) - (* Bessel functions of x of the second kind of order 1. *) - (*---------------------------------------------------------------*) - -(*==========================================================================*) - -PROCEDURE ceil(x : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* These function returns the smallest integral value that is not *) - (* less than x. For example, ceil(0.5) is 1.0, and ceil(-0.5) *) - (* is 0.0. *) - (*----------------------------------------------------------------*) - -PROCEDURE floor(x : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* These function returns the largest integral value that is not *) - (* greater than x. For example, floor(0.5) is 0.0, and *) - (* floor(-0.5) is -1.0. *) - (*----------------------------------------------------------------*) - -PROCEDURE modf( x : DOUBLE; - VAR iptr : DOUBLE): DOUBLE; - - (*----------------------------------------------------------------*) - (* The modf function breaks the argument x into an integral part *) - (* and a fractional part, each of which has the same sign as x. *) - (* The integral part is stored in the location pointed to by iptr *) - (*----------------------------------------------------------------*) - -PROCEDURE fmod(x : DOUBLE; - y : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The fmod function computes the floating-point remainder of *) - (* dividing x by y. The return value is x - n*y, where n is the *) - (* quotient of x / y, rounded toward zero to an integer. *) - (*----------------------------------------------------------------*) - -PROCEDURE logb(x : DOUBLE) : DOUBLE; - - (*----------------------------------------------------------------*) - (* These functions extract the exponent from the internal *) - (* floating-point representation of x and return it as a *) - (* floating-point value. The integer constant FLT_RADIX, *) - (* defined in <float.h>, indicates the radix used for *) - (* the systems floating-point representation. If FLT_RADIX is 2, *) - (* logb(x) is equal to floor(log2(x)), except that it is *) - (* probably faster. If x is subnormal, logb returns the exponent *) - (* x would have if it were normalized. *) - (*----------------------------------------------------------------*) - -PROCEDURE ilogb(x : DOUBLE) : INTEGER32; - - (*----------------------------------------------------------------*) - (* These functions return the exponent part of their argument as *) - (* a signed integer. *) - (*----------------------------------------------------------------*) - -PROCEDURE frexp( x : DOUBLE; - VAR exp : INTEGER32) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The frexp function is used to split the number x into a *) - (* normalized fraction and an exponent which is stored in exp. *) - (*----------------------------------------------------------------*) - -PROCEDURE ldexp(x : DOUBLE; - exp : INTEGER32) : DOUBLE; - - (*----------------------------------------------------------------*) - (* The ldexp function returns the result of multiplying the *) - (* floating-point number x by 2 raised to the power exp. *) - (*----------------------------------------------------------------*) - -END libmath. +{ ****************************************************************** + Types and constants - Error handling - Dynamic arrays + ****************************************************************** + The default real type is DOUBLE (8-byte real). + Other types may be selected by defining the symbols: + + SINGLEREAL (Single precision, 4 bytes) + EXTENDEDREAL (Extended precision, 10 bytes) + ****************************************************************** } + +unit utypes; + +interface + +{$i types.inc} + +{ ------------------------------------------------------------------ + Error handling + ------------------------------------------------------------------ } + +procedure SetErrCode(ErrCode : Integer); +{ Sets the error code } + +function DefaultVal(ErrCode : Integer; DefVal : Float) : Float; +{ Sets error code and default function value } + +function MathErr : Integer; +{ Returns the error code } + +{ ------------------------------------------------------------------ + Dynamic arrays + ------------------------------------------------------------------ } + +procedure SetAutoInit(AutoInit : Boolean); +{ Sets the auto-initialization of arrays } + +procedure DimVector(var V : PVector; Ub : Integer); +{ Creates floating point vector V[0..Ub] } + +procedure DimIntVector(var V : PIntVector; Ub : Integer); +{ Creates integer vector V[0..Ub] } + +procedure DimCompVector(var V : PCompVector; Ub : Integer); +{ Creates complex vector V[0..Ub] } + +procedure DimBoolVector(var V : PBoolVector; Ub : Integer); +{ Creates boolean vector V[0..Ub] } + +procedure DimStrVector(var V : PStrVector; Ub : Integer); +{ Creates string vector V[0..Ub] } + +procedure DimMatrix(var A : PMatrix; Ub1, Ub2 : Integer); +{ Creates floating point matrix A[0..Ub1, 0..Ub2] } + +procedure DimIntMatrix(var A : PIntMatrix; Ub1, Ub2 : Integer); +{ Creates integer matrix A[0..Ub1, 0..Ub2] } + +procedure DimCompMatrix(var A : PCompMatrix; Ub1, Ub2 : Integer); +{ Creates complex matrix A[0..Ub1, 0..Ub2] } + +procedure DimBoolMatrix(var A : PBoolMatrix; Ub1, Ub2 : Integer); +{ Creates boolean matrix A[0..Ub1, 0..Ub2] } + +procedure DimStrMatrix(var A : PStrMatrix; Ub1, Ub2 : Integer); +{ Creates string matrix A[0..Ub1, 0..Ub2] } + +procedure DelVector(var V : PVector; Ub : Integer); +{ Deletes floating point vector V[0..Ub] } + +procedure DelIntVector(var V : PIntVector; Ub : Integer); +{ Deletes integer vector V[0..Ub] } + +procedure DelCompVector(var V : PCompVector; Ub : Integer); +{ Deletes complex vector V[0..Ub] } + +procedure DelBoolVector(var V : PBoolVector; Ub : Integer); +{ Deletes boolean vector V[0..Ub] } + +procedure DelStrVector(var V : PStrVector; Ub : Integer); +{ Deletes string vector V[0..Ub] } + +procedure DelMatrix(var A : PMatrix; Ub1, Ub2 : Integer); +{ Deletes floating point matrix A[0..Ub1, 0..Ub2] } + +procedure DelIntMatrix(var A : PIntMatrix; Ub1, Ub2 : Integer); +{ Deletes integer matrix A[0..Ub1, 0..Ub2] } + +procedure DelCompMatrix(var A : PCompMatrix; Ub1, Ub2 : Integer); +{ Deletes complex matrix A[0..Ub1, 0..Ub2] } + +procedure DelBoolMatrix(var A : PBoolMatrix; Ub1, Ub2 : Integer); +{ Deletes boolean matrix A[0..Ub1, 0..Ub2] } + +procedure DelStrMatrix(var A : PStrMatrix; Ub1, Ub2 : Integer); +{ Deletes string matrix A[0..Ub1, 0..Ub2] } + +implementation + +{$IFDEF _16BIT}const{$ELSE}var{$ENDIF} + gAutoInit : Boolean = True; + +var + gErrCode : Integer; + +procedure SetErrCode(ErrCode : Integer); +begin + gErrCode := ErrCode; +end; + +function DefaultVal(ErrCode : Integer; DefVal : Float) : Float; +begin + SetErrCode(ErrCode); + DefaultVal := DefVal; +end; + +function MathErr : Integer; +begin + MathErr := gErrCode; +end; + +procedure SetAutoInit(AutoInit : Boolean); +begin + gAutoInit := AutoInit; +end; + +procedure DimVector(var V : PVector; Ub : Integer); +var + I : Integer; +begin + { Check bounds } + if (Ub < 0) or (Ub > MAX_FLT) then + begin + V := nil; + Exit; + end; + + { Allocate vector } + GetMem(V, (Ub + 1) * FltSize); + if V = nil then Exit; + + { Initialize vector } + if gAutoInit then + for I := 0 to Ub do + V^[I] := 0.0; +end; + +procedure DimIntVector(var V : PIntVector; Ub : Integer); +var + I : Integer; +begin + { Check bounds } + if (Ub < 0) or (Ub > MAX_INT) then + begin + V := nil; + Exit; + end; + + { Allocate vector } + GetMem(V, (Ub + 1) * IntSize); + if V = nil then Exit; + + { Initialize vector } + if gAutoInit then + for I := 0 to Ub do + V^[I] := 0; +end; + +procedure DimCompVector(var V : PCompVector; Ub : Integer); +var + I : Integer; +begin + { Check bounds } + if (Ub < 0) or (Ub > MAX_COMP) then + begin + V := nil; + Exit; + end; + + { Allocate vector } + GetMem(V, (Ub + 1) * CompSize); + if V = nil then Exit; + + { Initialize vector } + if gAutoInit then + for I := 0 to Ub do + begin + V^[I].X := 0.0; + V^[I].Y := 0.0; + end; +end; + +procedure DimBoolVector(var V : PBoolVector; Ub : Integer); +var + I : Integer; +begin + { Check bounds } + if (Ub < 0) or (Ub > MAX_BOOL) then + begin + V := nil; + Exit; + end; + + { Allocate vector } + GetMem(V, (Ub + 1) * BoolSize); + if V = nil then Exit; + + { Initialize vector } + if gAutoInit then + for I := 0 to Ub do + V^[I] := False; +end; + +procedure DimStrVector(var V : PStrVector; Ub : Integer); +var + I : Integer; +begin + { Check bounds } + if (Ub < 0) or (Ub > MAX_STR) then + begin + V := nil; + Exit; + end; + + { Allocate vector } + GetMem(V, (Ub + 1) * StrSize); + if V = nil then Exit; + + { Initialize vector } + if gAutoInit then + for I := 0 to Ub do + V^[I] := ''; +end; + +procedure DimMatrix(var A : PMatrix; Ub1, Ub2 : Integer); +var + I, J : Integer; + RowSize : Word; +begin + if (Ub1 < 0) or (Ub1 > MAX_VEC) or (Ub2 < 0) or (Ub2 > MAX_FLT) then + begin + A := nil; + Exit; + end; + + { Allocate matrix } + GetMem(A, (Ub1 + 1) * PtrSize); + if A = nil then Exit; + + { Size of a row } + RowSize := (Ub2 + 1) * FltSize; + + { Allocate each row } + for I := 0 to Ub1 do + begin + GetMem(A^[I], RowSize); + if A^[I] = nil then + begin + A := nil; + Exit; + end; + end; + + { Initialize matrix } + if gAutoInit then + for I := 0 to Ub1 do + for J := 0 to Ub2 do + A^[I]^[J] := 0.0; +end; + +procedure DimIntMatrix(var A : PIntMatrix; Ub1, Ub2 : Integer); +var + I, J : Integer; + RowSize : Word; +begin + { Check bounds } + if (Ub1 < 0) or (Ub1 > MAX_VEC) or (Ub2 < 0) or (Ub2 > MAX_INT) then + begin + A := nil; + Exit; + end; + + { Allocate matrix } + GetMem(A, (Ub1 + 1) * PtrSize); + if A = nil then Exit; + + { Size of a row } + RowSize := (Ub2 + 1) * IntSize; + + { Allocate each row } + for I := 0 to Ub1 do + begin + GetMem(A^[I], RowSize); + if A^[I] = nil then + begin + A := nil; + Exit; + end; + end; + + { Initialize matrix } + if gAutoInit then + for I := 0 to Ub1 do + for J := 0 to Ub2 do + A^[I]^[J] := 0; +end; + +procedure DimCompMatrix(var A : PCompMatrix; Ub1, Ub2 : Integer); +var + I, J : Integer; + RowSize : Word; +begin + { Check bounds } + if (Ub1 < 0) or (Ub1 > MAX_VEC) or (Ub2 < 0) or (Ub2 > MAX_COMP) then + begin + A := nil; + Exit; + end; + + { Allocate matrix } + GetMem(A, (Ub1 + 1) * PtrSize); + if A = nil then Exit; + + { Size of a row } + RowSize := (Ub2 + 1) * CompSize; + + { Allocate each row } + for I := 0 to Ub1 do + begin + GetMem(A^[I], RowSize); + if A^[I] = nil then + begin + A := nil; + Exit; + end; + end; + + { Initialize matrix } + if gAutoInit then + for I := 0 to Ub1 do + for J := 0 to Ub2 do + begin + A^[I]^[J].X := 0.0; + A^[I]^[J].Y := 0.0; + end; +end; + +procedure DimBoolMatrix(var A : PBoolMatrix; Ub1, Ub2 : Integer); +var + I, J : Integer; + RowSize : Word; +begin + { Check bounds } + if (Ub1 < 0) or (Ub1 > MAX_VEC) or (Ub2 < 0) or (Ub2 > MAX_BOOL) then + begin + A := nil; + Exit; + end; + + { Allocate matrix } + GetMem(A, (Ub1 + 1) * PtrSize); + if A = nil then Exit; + + { Size of a row } + RowSize := (Ub2 + 1) * BoolSize; + + { Allocate each row } + for I := 0 to Ub1 do + begin + GetMem(A^[I], RowSize); + if A^[I] = nil then + begin + A := nil; + Exit; + end; + end; + + { Initialize matrix } + if gAutoInit then + for I := 0 to Ub1 do + for J := 0 to Ub2 do + A^[I]^[J] := False; +end; + +procedure DimStrMatrix(var A : PStrMatrix; Ub1, Ub2 : Integer); +var + I, J : Integer; + RowSize : Word; +begin + { Check bounds } + if (Ub1 < 0) or (Ub1 > MAX_VEC) or (Ub2 < 0) or (Ub2 > MAX_STR) then + begin + A := nil; + Exit; + end; + + { Allocate matrix } + GetMem(A, (Ub1 + 1) * PtrSize); + if A = nil then Exit; + + { Size of a row } + RowSize := (Ub2 + 1) * StrSize; + + { Allocate each row } + for I := 0 to Ub1 do + begin + GetMem(A^[I], RowSize); + if A^[I] = nil then + begin + A := nil; + Exit; + end; + end; + + { Initialize matrix } + if gAutoInit then + for I := 0 to Ub1 do + for J := 0 to Ub2 do + A^[I]^[J] := ''; +end; + +procedure DelVector(var V : PVector; Ub : Integer); +begin + if V <> nil then + begin + FreeMem(V, (Ub + 1) * FltSize); + V := nil; + end; +end; + +procedure DelIntVector(var V : PIntVector; Ub : Integer); +begin + if V <> nil then + begin + FreeMem(V, (Ub + 1) * IntSize); + V := nil; + end; +end; + +procedure DelCompVector(var V : PCompVector; Ub : Integer); +begin + if V <> nil then + begin + FreeMem(V, (Ub + 1) * CompSize); + V := nil; + end; +end; + +procedure DelBoolVector(var V : PBoolVector; Ub : Integer); +begin + if V <> nil then + begin + FreeMem(V, (Ub + 1) * BoolSize); + V := nil; + end; +end; + +procedure DelStrVector(var V : PStrVector; Ub : Integer); +begin + if V <> nil then + begin + FreeMem(V, (Ub + 1) * StrSize); + V := nil; + end; +end; + +procedure DelMatrix(var A : PMatrix; Ub1, Ub2 : Integer); +var + I : Integer; + RowSize : Word; +begin + if A <> nil then + begin + RowSize := (Ub2 + 1) * FltSize; + for I := Ub1 downto 0 do + FreeMem(A^[I], RowSize); + FreeMem(A, (Ub1 + 1) * PtrSize); + A := nil; + end; +end; + +procedure DelIntMatrix(var A : PIntMatrix; Ub1, Ub2 : Integer); +var + I : Integer; + RowSize : Word; +begin + if A <> nil then + begin + RowSize := (Ub2 + 1) * IntSize; + for I := Ub1 downto 0 do + FreeMem(A^[I], RowSize); + FreeMem(A, (Ub1 + 1) * PtrSize); + A := nil; + end; +end; + +procedure DelCompMatrix(var A : PCompMatrix; Ub1, Ub2 : Integer); +var + I : Integer; + RowSize : Word; +begin + if A <> nil then + begin + RowSize := (Ub2 + 1) * CompSize; + for I := Ub1 downto 0 do + FreeMem(A^[I], RowSize); + FreeMem(A, (Ub1 + 1) * PtrSize); + A := nil; + end; +end; + +procedure DelBoolMatrix(var A : PBoolMatrix; Ub1, Ub2 : Integer); +var + I : Integer; + RowSize : Word; +begin + if A <> nil then + begin + RowSize := (Ub2 + 1) * BoolSize; + for I := Ub1 downto 0 do + FreeMem(A^[I], RowSize); + FreeMem(A, (Ub1 + 1) * PtrSize); + A := nil; + end; +end; + +procedure DelStrMatrix(var A : PStrMatrix; Ub1, Ub2 : Integer); +var + I : Integer; + RowSize : Word; +begin + if A <> nil then + begin + RowSize := (Ub2 + 1) * StrSize; + for I := Ub1 downto 0 do + FreeMem(A^[I], RowSize); + FreeMem(A, (Ub1 + 1) * PtrSize); + A := nil; + end; +end; + +end.