复数运算单元

UNIT ComplexOps; { see demo at end }

{This UNIT provides complex arithmetic and transcendental functions.

(C) Copyright 1990, 1992, Earl F. Glynn, Overland Park, KS. Compuserve 73257,3527.
All rights reserved. This program may be freely distributed only for
non-commercial use.

Some ideas in this UNIT were borrowed from "A Pascal Tool for Complex
Numbers", Journal of Pascal, Ada, & Modula-2, May/June 1985, pp. 23-29.
Many complex formulas were taken from Chapter 4, "Handbook of Mathematical
Functions" (Ninth Printing), Abramowitz and Stegun (editors), Dover, 1972.}

INTERFACE

TYPE
RealType = DOUBLE;
ComplexForm = (polar,rectangular);
Complex =
RECORD
CASE form: ComplexForm OF
rectangular: (x,y : RealType); {z = x + y*i}
polar : (r,theta: RealType); {z = r*CIS(theta)}
END; {where CIS(theta) = COS(theta) + SIN(theta)*i}
{ theta = -PI..PI (in canonical form)}

CONST
MaxTerm : BYTE = 35;
EpsilonSqr : RealType = 1.0E-20;
Infinity : RealType = 1.0E25; {virtual infinity}

{complex number definition/conversion/output}
PROCEDURE CConvert (VAR z: Complex; f: ComplexForm);
PROCEDURE CSet (VAR z: Complex; a,b: RealType; f: ComplexForm);
FUNCTION CStr (z: Complex; w,d: BYTE; f: ComplexForm): STRING;

{complex arithmetic}
PROCEDURE CAdd (VAR z: Complex; a,b: Complex); {z = a + b}
PROCEDURE CDiv (VAR z: Complex; a,b: Complex); {z = a / b}
PROCEDURE CMult (VAR z: Complex; a,b: Complex); {z = a * b}
PROCEDURE CSub (VAR z: Complex; a,b: Complex); {z = a - b}
PROCEDURE CNeg (VAR z: Complex; a : Complex); {z = -a }

{complex natural log, exponential}
PROCEDURE CLn (VAR fn : Complex; z: Complex); {fn = ln(z) }
PROCEDURE CExp (VAR z : Complex; a: Complex); {z = exp(a)}
PROCEDURE CPwr (VAR z : Complex; a,b: Complex); {z = a^b }

{complex trig functions}
PROCEDURE CCos (VAR z: Complex; a: Complex); {z = cos(a)}
PROCEDURE CSin (VAR z: Complex; a: Complex); {z = sin(a)}
PROCEDURE CTan (VAR z: Complex; a: Complex); {z = tan(a)}

PROCEDURE CSec (VAR z: Complex; a: Complex); {z = sec(a)}
PROCEDURE CCsc (VAR z: Complex; a: Complex); {z = csc(a)}
PROCEDURE CCot (VAR z: Complex; a: Complex); {z = cot(a)}

{complex hyperbolic functions}
PROCEDURE CCosh (VAR z: Complex; a: Complex); {z = cosh(a)}
PROCEDURE CSinh (VAR z: Complex; a: Complex); {z = sinh(a)}
PROCEDURE CTanh (VAR z: Complex; a: Complex); {z = tanh(a)}

PROCEDURE CSech (VAR z: Complex; a: Complex); {z = sech(a)}
PROCEDURE CCsch (VAR z: Complex; a: Complex); {z = csch(a)}
PROCEDURE CCoth (VAR z: Complex; a: Complex); {z = coth(a)}

{miscellaneous complex functions}
FUNCTION CAbs (z: Complex): RealType; {CAbs = |z|}
FUNCTION CAbsSqr (z: Complex): RealType; {CAbsSqr = |z|^2}
PROCEDURE CIntPwr (VAR z: Complex; a: Complex; n: INTEGER); {z = a^n}
PROCEDURE CRealPwr (VAR z: Complex; a: Complex; x: RealType); {z = a^x}
PROCEDURE CConjugate (VAR z: Complex; a: Complex); {z = a*}
PROCEDURE CSqrt (VAR z: Complex; a: Complex); {z = SQRT(a)}
PROCEDURE CRoot (VAR z: Complex; a: Complex; k,n: WORD); {z = a^(1/n)}

{complex Bessel functions of order zero}
PROCEDURE CI0 (VAR sum: Complex; z: Complex); {sum = I0(z)}
PROCEDURE CJ0 (VAR sum: Complex; z: Complex); {sum = J0(z)}

PROCEDURE CLnGamma (VAR z: Complex; a: Complex);
PROCEDURE CGamma (VAR z: Complex; a: Complex);

{treat "fuzz" of real numbers}
PROCEDURE CDeFuzz (VAR z: Complex);
FUNCTION DeFuzz (x: RealType): RealType;
PROCEDURE SetFuzz (value: RealType);

{miscellaneous}
FUNCTION FixAngle (theta: RealType): RealType; {-PI < theta <= PI}
FUNCTION Pwr (x,y: RealType): RealType; {Pwr = x^y}
FUNCTION Log10 (x: RealType): RealType;
FUNCTION LongMod (l1,l2: LongInt): LongInt;
FUNCTION Cosh (x: RealType): RealType;
FUNCTION Sinh (x: RealType): RealType;

IMPLEMENTATION

VAR
fuzz : RealType;
Cone : Complex;
Cinfinity: Complex;
Czero : Complex;
hLnPI : RealType;
hLn2PI : RealType;
ln2 : RealType;

{complex number definition/conversion/output}
PROCEDURE CConvert (VAR z: Complex; f: ComplexForm);
VAR a: Complex;
BEGIN
IF z.form = f
THEN CDeFuzz (z)
ELSE BEGIN
CASE z.form OF
polar: {polar-to-rectangular conversion}
BEGIN
a.form := rectangular;
a.x := z.r * COS(z.theta);
a.y := z.r * SIN(z.theta)
END;
rectangular: {rectangular-to-polar conversion}
BEGIN
a.form := polar;
IF DeFuzz(z.x) = 0.0
THEN BEGIN
IF DeFuzz(z.y) = 0.0
THEN BEGIN
a.r := 0.0;
a.theta := 0.0
END
ELSE
IF z.y > 0.0
THEN BEGIN
a.r := z.y;
a.theta := 0.5*PI
END
ELSE BEGIN
a.r := -z.y;
a.theta := -0.5*PI
END
END
ELSE BEGIN
a.r := CAbs(z);
a.theta := ARCTAN(z.y/z.x); {4th/1st quadrant -PI/2..PI/2}
IF z.x < 0.0 {2nd/3rd quadrants}
THEN
IF z.y >= 0.0
THEN a.theta := PI + a.theta {2nd quadrant: PI/2..PI}
ELSE a.theta := -PI + a.theta {3rd quadrant: -PI..-PI/2}
END
END;
END;
CDeFuzz (a);
z := a
END
END {CConvert};

PROCEDURE CSet (VAR z: Complex; a,b: RealType; f: ComplexForm);
BEGIN
z.form := f;
CASE f OF
polar:
BEGIN
z.r := a;
z.theta := b
END;
rectangular:
BEGIN
z.x := a;
z.y := b
END;
END
END {CSet};

FUNCTION CStr (z: Complex; w,d: BYTE; f: ComplexForm): STRING;
VAR s1,s2: STRING;
BEGIN
CConvert (z,f);
CASE f OF
polar:
BEGIN
Str (z.r:w:d, s1);
Str (z.theta:w:d, s2);
CStr := s1+'*CIS('+s2+')'
END;
rectangular:
BEGIN
Str (z.x:w:d, s1);
Str (ABS(z.y):w:d, s2);
IF z.y >= 0
THEN CStr := s1+' +'+s2+'i'
ELSE CStr := s1+' -'+s2+'i'
END
END
END {CStr};

{complex arithmetic}
PROCEDURE CAdd (VAR z: Complex; a,b: Complex); {z = a + b}
BEGIN {complex addition}
CConvert (a,rectangular);
CConvert (b,rectangular);
z.form := rectangular;
z.x := a.x + b.x; {real part}
z.y := a.y + b.y; {imaginary part}
END {CAdd};

PROCEDURE CDiv (VAR z: Complex; a,b: Complex); {z = a / b}
VAR temp: RealType;
BEGIN
CConvert (b,a.form); {arbitrarily convert one to type of other}
z.form := a.form;
CASE a.form OF
polar:
BEGIN
z.r := a.r / b.r;
z.theta := FixAngle(a.theta - b.theta)
END;
rectangular:
BEGIN
temp := SQR(b.x) + SQR(b.y);
z.x := (a.x*b.x + a.y*b.y) / temp;
z.y := (a.y*b.x - a.x*b.y) / temp
END
END
END {CDiv};

PROCEDURE CMult(VAR z: Complex; a,b: Complex); {z = a * b}
BEGIN
CConvert (b,a.form); {arbitrarily convert one to type of other}
z.form := a.form;
CASE a.form OF
polar:
BEGIN
z.r := a.r * b.r;
z.theta := FixAngle(a.theta + b.theta)
END;
rectangular:
BEGIN
z.x := a.x*b.x - a.y*b.y;
z.y := a.x*b.y + a.y*b.x
END
END
END {CMult};

PROCEDURE CSub (VAR z: Complex; a,b: Complex); {z = a - b}
BEGIN {complex subtraction}
CConvert (a,rectangular);
CConvert (b,rectangular);
z.form := rectangular;
z.x := a.x - b.x; {real part}
z.y := a.y - b.y; {imaginary part}
END {CSub};

PROCEDURE CNeg (VAR z: Complex; a : Complex); {z = -a }
BEGIN
z.form := a.form;
CASE a.form OF
polar:
BEGIN
z.r := a.r;
z.theta := FixAngle(a.theta + PI)
END;
rectangular:
BEGIN
z.x := -a.x;
z.y := -a.y
END
END
END {CNeg};
{complex natural log, exponential}
PROCEDURE CLn (VAR fn: Complex; z: Complex); {fn = ln(z)}
BEGIN {Abramowitz formula 4.1.2 on p. 67}
CConvert (z,polar);
fn.form := rectangular;
fn.x := LN(z.r);
fn.y := FixAngle(z.theta)
END {CLn}; {principal value only}

PROCEDURE CExp (VAR z : Complex; a: Complex); {z = exp(a)}
VAR
temp: RealType;
BEGIN {Euler's Formula: Abramowitz formula 4.3.47 on p. 74}
CConvert (a,rectangular);
temp := EXP(a.x);
CSet (z, temp*COS(a.y),temp*SIN(a.y), rectangular)
END {CExp};

PROCEDURE CPwr (VAR z : Complex; a,b: Complex); {z = a^b }
VAR
blna,lna: Complex;
BEGIN {Abramowitz formula 4.2.7 on p. 69}
CDeFuzz (a);
CDeFuzz (b);
IF CAbsSqr(a) = 0.0
THEN
IF (CAbsSqr(b) = 0.0)
THEN z := Cone {lim a^a = 1 as a -> 0}
ELSE z := Czero {0^b = 0, b <> 0}
ELSE BEGIN
CLn (lna,a);
CMult(blna,b,lna);
CExp (z, blna)
END
END {CPwr};
{complex trig functions}
PROCEDURE CCos (VAR z: Complex; a: Complex); {z = cos(a)}
BEGIN {Abramowitz formula 4.3.56 on p. 74}
CConvert (a,rectangular);
CSet (z, COS(a.x)*COSH(a.y), -SIN(a.x)*SINH(a.y), rectangular)
END {CCos};

PROCEDURE CSin (VAR z: Complex; a: Complex); {z = sin(a)}
BEGIN {Abramowitz formula 4.3.55 on p. 74}
CConvert (a,rectangular);
CSet (z, SIN(a.x)*COSH(a.y), COS(a.x)*SINH(a.y), rectangular)
END {CSin};

PROCEDURE CTan (VAR z: Complex; a: Complex); {z = tan(a)}
VAR
temp: RealType;
BEGIN {Abramowitz formula 4.3.57 on p. 74}
CConvert (a,rectangular);
temp := COS(2.0*a.x) + COSH(2.0*a.y);
IF DeFuzz(temp) <> 0.0
THEN BEGIN
CSet (z,SIN(2.0*a.x)/temp,SINH(2.0*a.y)/temp,rectangular)
END
ELSE z := Cinfinity
END {CTan};

PROCEDURE CSec (VAR z: Complex; a: Complex); {z = sec(a)}
VAR
temp: Complex;
BEGIN {Abramowitz formula 4.3.5 on p. 72}
CCos (temp, a);
IF DeFuzz( Cabs(temp) ) <> 0.0
THEN CDiv (z, Cone,temp)
ELSE z := Cinfinity
END {CSec};

PROCEDURE CCsc (VAR z: Complex; a: Complex); {z = csc(a)}
VAR
temp: Complex;
BEGIN {Abramowitz formula 4.3.4 on p. 72}
CSin (temp, a);
IF DeFuzz( Cabs(temp) ) <> 0.0
THEN CDiv (z, Cone,temp)
ELSE z := Cinfinity
END {CCsc};

PROCEDURE CCot (VAR z: Complex; a: Complex); {z = cot(a)}
VAR
temp: RealType;
BEGIN {Abramowitz formula 4.3.58 on p. 74}
CConvert (a,rectangular);
temp := COSH(2.0*a.y) - COS(2.0*a.x);
IF DeFuzz(temp) <> 0.0
THEN BEGIN
CSet (z,SIN(2.0*a.x)/temp,-SINH(2.0*a.y)/temp,rectangular)
END
ELSE z := Cinfinity
END {CCot};

{complex hyperbolic functions}
PROCEDURE CCosh (VAR z: Complex; a: Complex); {z = cosh(a)}
BEGIN {Abramowitz formula 4.5.50 on p. 84}
CConvert (a,rectangular);
CSet (z, COSH(a.x)*COS(a.y), SINH(a.x)*SIN(a.y), rectangular)
END {CCosh};

PROCEDURE CSinh (VAR z: Complex; a: Complex); {z = sinh(a)}
BEGIN {Abramowitz formula 4.5.49 on p.84}
CConvert (a,rectangular);
CSet (z, SINH(a.x)*COS(a.y), COSH(a.x)*SIN(a.y), rectangular)
END {CSinh};

PROCEDURE CTanh (VAR z: Complex; a: Complex); {z = tanh(a)}
VAR
temp: RealType;
BEGIN {Abramowitz formula 4.5.51 on p. 84}
CConvert (a,rectangular);
temp := COSH(2.0*a.x) + COS(2.0*a.y);
IF DeFuzz(temp) <> 0.0
THEN BEGIN
CSet (z,SINH(2.0*a.x)/temp,SIN(2.0*a.y)/temp,rectangular)
END
ELSE z := Cinfinity
END {CTanh};

PROCEDURE CSech (VAR z: Complex; a: Complex); {z = sech(a)}
VAR
temp: Complex;
BEGIN {Abramowitz formula 4.5.5 on p. 83}
CCosh (temp, a);
IF DeFuzz( Cabs(temp) ) <> 0.0
THEN CDiv (z, Cone,temp)
ELSE z := Cinfinity
END {CSec};

PROCEDURE CCsch (VAR z: Complex; a: Complex); {z = csch(a)}
VAR
temp: Complex;
BEGIN {Abramowitz formula 4.5.4 on p. 83}
CSinh (temp, a);
IF DeFuzz( Cabs(temp) ) <> 0.0
THEN CDiv (z, Cone,temp)
ELSE z := Cinfinity
END {CCsch};

PROCEDURE CCoth (VAR z: Complex; a: Complex); {z = coth(a)}
VAR
temp: RealType;
BEGIN {Abramowitz formula 4.5.52 on p. 84}
CConvert (a,rectangular);
temp := COSH(2.0*a.x) - COS(2.0*a.y);
IF DeFuzz(temp) <> 0.0
THEN BEGIN
CSet (z,SINH(2.0*a.x)/temp,-SIN(2.0*a.y)/temp,rectangular)
END
ELSE z := Cinfinity
END {CCoth};

{miscellaneous complex functions}
FUNCTION CAbs (z: Complex): RealType;
BEGIN
CASE z.form OF
rectangular: CAbs := SQRT( SQR(z.x) + SQR(z.y) );
polar: CAbs := z.r
END
END;

FUNCTION CAbsSqr (z: Complex): RealType; {CAbsSqr = |z|^2}
BEGIN
CASE z.form OF
rectangular: CAbsSqr := SQR(z.x) + SQR(z.y);
polar: CAbsSqr := SQR(z.r)
END
END {CAbsSqr};

PROCEDURE CIntPwr (VAR z: Complex; a: Complex; n: INTEGER); {z = a^n}
{CIntPwr directly applies DeMoivre's theorem to calculate
an integer power of a complex number. The formula holds
for both positive and negative values of 'n'.}
BEGIN
IF CAbsSqr(a) = 0.0
THEN
IF n = 0
THEN z := Cone {lim a^a = 1 as a -> 0}
ELSE z := Czero {0^n = 0, except for 0^0=1}
ELSE BEGIN
CConvert (a,polar);
z.form := polar;
z.r := Pwr(a.r,n);
z.theta := FixAngle(n*a.theta)
END
END {CIntPwr};

PROCEDURE CRealPwr (VAR z: Complex; a: Complex; x: RealType); {z = a^x}
BEGIN
IF CAbsSqr(a) = 0.0
THEN
IF Defuzz(x) = 0.0
THEN z := Cone {lim a^a = 1 as a -> 0}
ELSE z := Czero {0^n = 0, except for 0^0=1}
ELSE BEGIN
CConvert (a,polar);
z.form := polar;
z.r := Pwr(a.r,x);
z.theta := FixAngle(x*a.theta)
END
END {CRealPwr};

PROCEDURE CConjugate (VAR z: Complex; a: Complex); {z = a*}
BEGIN
z.form := a.form;
CASE a.form OF
polar:
BEGIN
z.r := a.r;
z.theta := FixAngle(-a.theta)
END;
rectangular:
BEGIN
z.x := a.x;
z.y := -a.y
END;
END;
END {CConjugate};

PROCEDURE CSqrt (VAR z: Complex; a: Complex); {z = SQRT(a)}
BEGIN
CRoot (z, a, 0,2) {return only one of the two values}
END {CSqrt};
{z = a^(1/n), n > 0}
PROCEDURE CRoot (VAR z: Complex; a: Complex; k,n: WORD);
{CRoot can calculate all 'n' roots of 'a' by varying 'k' from 0..n-1.}
{This is another application of DeMoivre's theorem. See CIntPwr.}
BEGIN
IF CAbs(a) = 0.0
THEN z := Czero {0^z = 0, except 0^0 is undefined}
ELSE BEGIN
CConvert (a,polar);
z.form := polar;
z.r := Pwr(a.r,1.0/n);
z.theta := FixAngle((a.theta + k*2.0*PI)/n)
END;
END;

PROCEDURE CI0 (VAR sum: Complex; z: Complex);
VAR
i : BYTE;
SizeSqr: RealType;
term : Complex;
zSQR25 : Complex;
BEGIN
CConvert (z,rectangular);
sum := Cone; {term 0}
Cmult(zSQR25, z,z);
zSQR25.x := 0.25 * zSQR25.x;
zSQR25.y := 0.25 * zSQR25.y; {瑉^2}
term := zSQR25;
CAdd (sum, sum,zSQR25); {term 1}
i := 1;
REPEAT
CMult(term,zSQR25,term);
INC (i);
term.x := term.x / SQR(i);
term.y := term.y / SQR(i);
CAdd (sum, sum,term); {sum := sum + term}
SizeSqr := SQR(term.x) + SQR(term.y)
UNTIL (i > MaxTerm) OR (SizeSqr < EpsilonSqr)
END;

PROCEDURE CJ0 (VAR sum: Complex; z: Complex);
VAR
addflag: BOOLEAN;
i : BYTE;
SizeSqr: RealType;
term : Complex;
zSQR25 : Complex;
BEGIN
CConvert (z,rectangular);
sum := Cone; {term 0}
Cmult(zSQR25, z,z);
zSQR25.x := 0.25 * zSQR25.x;
zSQR25.y := 0.25 * zSQR25.y; {瑉^2}
term := zSQR25;
CSub (sum, sum,zSQR25); {term 1}
addflag := FALSE;
i := 1;
REPEAT
CMult(term,zSQR25,term);
INC (i);
addflag := NOT addflag;
term.x := term.x / SQR(i);
term.y := term.y / SQR(i);
IF addflag
THEN CAdd (sum, sum,term) {sum := sum + term}
ELSE CSub (sum, sum,term); {sum := sum - term}
SizeSqr := SQR(term.x) + SQR(term.y)
UNTIL (i > MaxTerm) OR (SizeSqr < EpsilonSqr)
END {CJ0};

PROCEDURE CApproxLnGamma (VAR sum: Complex; z: Complex);
{This is the approximation used in the National Bureau of
Standards "Table of the Gamma Function for Complex Arguments,"
Applied Mathematics Series 34, 1954. The NBS table was created
using this approximation over the area 9 ?Re(z) ?10 and
0 ?Im(z) ?10. Other table values were computed using the
relationship ln ?z+1) = ln z + ln ?z).}
CONST
c: ARRAY[1..8] OF RealType
= (1/12, -1/360, 1/1260, -1/1680, 1/1188, -691/360360,
1/156, -3617/122400);
VAR
i : WORD;
powers: ARRAY[1..8] OF Complex;
temp1 : Complex;
temp2 : Complex;
BEGIN
CConvert (z,rectangular);
CLn (temp1,z); {ln(z}
CSet (temp2, z.x-0.5, z.y, rectangular); {z - 0.5}
CMult(sum, temp1,temp2); {(z - 0.5)*ln(z)}
CSub (sum, sum,z); {(z - 0.5)*ln(z) - z}
sum.x := sum.x + hLn2PI;

temp1 := Cone;
CDiv (powers[1], temp1, z); {z^-1}
CMult(temp2, powers[1],powers[1]); {z^-2}
FOR i := 2 TO 8 DO
CMult(powers[i], powers[i-1],temp2);
FOR i := 8 DOWNTO 1 DO BEGIN
CSet (temp1, c[i]*powers[i].x, c[i]*powers[i].y, rectangular);
CAdd (sum, sum,temp1);
END
END {CApproxLnGamma};

PROCEDURE CLnGamma (VAR z: Complex; a: Complex);
VAR
lna : Complex;
temp: Complex;
BEGIN
CConvert (a, rectangular);

IF (a.x <= 0.0) AND (DeFuzz(a.y) = 0.0)
THEN
IF DeFuzz(INT(a.x-1E-8) - a.x) = 0.0 {negative integer?}
THEN BEGIN
z := Cinfinity;
EXIT
END;

IF a.y < 0.0 {3rd or 4th quadrant?}
THEN BEGIN
CConjugate (a, a);
CLnGamma (z, a); {try again in 1st or 2nd quadrant}
CConjugate (z, z) {left this out! 1/3/91}
END
ELSE BEGIN
IF a.x < 9.0 {"left" of NBS table range}
THEN BEGIN
CLn (lna, a);
CSet (a, a.x+1.0, a.y, rectangular);
CLnGamma (temp,a);
CSub (z, temp,lna)
END
ELSE CApproxLnGamma (z,a) {NBS table range: 9 ?Re(z) ?10}
END
END {CLnGamma};

PROCEDURE CGamma (VAR z: Complex; a: Complex);
VAR
lnz: Complex;
BEGIN
CLnGamma (lnz,a);
IF lnz.x >= 75.0 {arbitrary cutoff for infinity}
THEN z := Cinfinity
ELSE
IF lnz.x < -200.0
THEN z := Czero {avoid underflow}
ELSE CExp (z, lnz);
END {CGamma};

{treat "fuzz" of real numbers}
PROCEDURE CDeFuzz (VAR z: Complex);
BEGIN
CASE z.form OF
rectangular:
BEGIN
z.x := DeFuzz(z.x);
z.y := DeFuzz(z.y);
END;
polar:
BEGIN
z.r := DeFuzz(z.r);
IF z.r = 0.0
THEN z.theta := 0.0 {canonical form when radius is zero}
ELSE z.theta := DeFuzz(z.theta)
END
END
END {CDeFuzz};

FUNCTION DeFuzz (x: RealType): RealType;
BEGIN
IF ABS(x) < fuzz
THEN Defuzz := 0
ELSE Defuzz := x
END {Defuzz};

PROCEDURE SetFuzz (value: RealType);
BEGIN
fuzz := value
END {SetFuzz};

{miscellaneous}
FUNCTION FixAngle (theta: RealType): RealType; {-PI < theta <= PI}
BEGIN
WHILE theta > PI DO
theta := theta - 2.0*PI;
WHILE theta <= -PI DO
theta := theta + 2.0*PI;
FixAngle := DeFuzz(theta)
END {FixAngle};

FUNCTION Pwr (x,y: RealType): RealType; {Pwr = x^y}
BEGIN
IF DeFuzz(x) = 0.0
THEN
IF DeFuzz(y) = 0.0
THEN Pwr := 1.0 {0^0 = 1 (i.e., lim x^x = 1 as x -> 0)}
ELSE Pwr := 0.0
ELSE Pwr := EXP( LN(x)*y )
END {Pwr};

FUNCTION Log10 (x: RealType): RealType;
BEGIN
Log10 := LN(x) / LN(10)
END {Log10};

FUNCTION LongMod (l1,l2: LongInt): LongInt;
BEGIN
LongMod := l1 - l2*(l1 DIV l2)
END {LongMod};

FUNCTION Cosh (x: RealType): RealType;
BEGIN
Cosh := 0.5*( EXP(x) + EXP(-x) )
END {Cosh};

FUNCTION Sinh (x: RealType): RealType;
BEGIN
Sinh := 0.5*( EXP(x) - EXP(-x) )
END {Sinh};

end.