{\$N+,E+} {Use math coprocessor, if any, emulate otherwise.}
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.
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}
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}

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;                  {CAbs = |z|}
BEGIN
CASE z.form OF
rectangular:  CAbs := SQRT( SQR(z.x) + SQR(z.y) );
polar:        CAbs := z.r
END
END {CAbs};

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 {CRoot};

{complex Bessel functions of order zero}
PROCEDURE CI0   (VAR sum:  Complex; z:  Complex);  {sum = I0(z)}
{I0(z) = ä ( (­z^2)^k / (k!)^2 ), k=0,1,2,...,ì}
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;       {­z^2}
term := zSQR25;
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 {CI0};

PROCEDURE CJ0   (VAR sum:  Complex; z:  Complex);  {sum = J0(z)}
{J0(z) = ä ( (-1)^k * (­z^2)^k / (k!)^2 ), k=0,1,2,...,ì}
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;       {­z^2}
term := zSQR25;
CSub (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);
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);
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};

BEGIN
SetFuzz (1.0E-12);
CSet ( Cone, 1.0, 0.0, rectangular);
CSet (Czero, 0.0, 0.0, rectangular);
CSet (Cinfinity, Infinity, Infinity, rectangular);
hLnPI := 0.5*(LN(PI));
hLn2PI := 0.5*(LN(2.0*PI));
ln2 := LN(2.0)
END.

{-------------  DEMO ----  CUT HERE ------ }
{\$N+,E+}
PROGRAM cdemo;

{This PROGRAM demonstrates the use of the ComplexOps UNIT.

(C) Copyright 1990, 1992, Earl F. Glynn, Overland Park, KS.  Compuserve 73257,3527.
non-commercial use.}

USES ComplexOps;

VAR
a      :  ARRAY[1..22] OF Complex;
csave  :  ARRAY[1..22] OF Complex;
k,m    :  WORD;
n      :  INTEGER;
x,y    :  RealType;
z,z1,z2:  Complex;

BEGIN

WRITELN ('Demo ComplexOPs PROCEDUREs and FUNCTIONs');
WRITELN;
WRITELN ('  Notes:  1.  CIS(w) = COS(w) + i*SIN(w), w = -PI..PI');
WRITELN ('          2.  z = x + i*y');
WRITELN;
WRITELN;

CSet (a[ 1],  0.0,  0.0, rectangular);
CSet (a[ 2],  0.5,  0.5, rectangular);
CSet (a[ 3], -0.5,  0.5, rectangular);
CSet (a[ 4], -0.5, -0.5, rectangular);
CSet (a[ 5],  0.5, -0.5, rectangular);
CSet (a[ 6],  1.0,  0.0, rectangular);
CSet (a[ 7],  1.0,  1.0, rectangular);
CSet (a[ 8],  0.0,  1.0, rectangular);
CSet (a[ 9], -1.0,  1.0, rectangular);
CSet (a[10], -1.0,  0.0, rectangular);
CSet (a[11], -1.0, -1.0, rectangular);
CSet (a[12],  0.0, -1.0, rectangular);
CSet (a[13],  1.0, -1.0, rectangular);
CSet (a[14],   5.,   0., rectangular);
CSet (a[15],   5.,   3., rectangular);
CSet (a[16],   0.,   3., rectangular);
CSet (a[17],  -5.,   3., rectangular);
CSet (a[18],  -5.,   0., rectangular);
CSet (a[19],  -5.,  -3., rectangular);
CSet (a[20],   0.,  -3., rectangular);
CSet (a[21],  -5.,  -3., rectangular);
CSet (a[22], -20.,  20., rectangular);

WRITELN ('Complex number definition/conversion/output:  CSet/CConvert/CStr');
WRITELN;
WRITELN ('   z rectangular':25,'z polar':28);
WRITELN ('     ---------------------------   ',
'-----------------------------');
FOR k := 1 TO 21 DO
WRITELN (k:3,'  ',CStr(a[k],12,8,rectangular),'  ',
CStr(a[k],12,8,polar));
WRITELN;
WRITELN;

WRITELN ('Complex arithmetic:  CAdd, CSub, CMult, CDiv');
WRITELN;

CSet (z1,  1, 1, rectangular);
WRITELN ('Let z1 = ':12,CStr(z1,8,3,rectangular):20,' or ',
CStr(z1,8,3,polar));
CSet (z2, SQRT(3), -1, rectangular);
WRITELN ('z2 = ':12,CStr(z2,8,3,rectangular):20,' or ',
CStr(z2,8,3,polar));
WRITELN;

WRITELN ('z1 + z2 = ':12,CStr(z,8,3,rectangular):20,' or ',
CStr(z,8,3,polar));

CSub  (z,z1,z2);
WRITELN ('z1 - z2 = ':12,CStr(z,8,3,rectangular):20,' or ',
CStr(z,8,3,polar));

CMult (z,z1,z2);
WRITELN ('z1 * z2 = ':12,CStr(z,8,3,rectangular):20,' or ',
CStr(z,8,3,polar));

CDiv  (z,z1,z2);
WRITELN ('z1 / z2 = ':12,CStr(z,8,3,rectangular):20,' or ',
CStr(z,8,3,polar));
WRITELN;
WRITELN;

WRITELN ('Complex natural logarithm:  CLn = LN(z)');
WRITELN;
WRITELN ('  Notes:  1.  LN(z) is multivalued.');
WRITELN (' ':9,' 2.  Any multiple of 2*PI*i could be added to/',
'subtracted from LN(z).');
WRITELN (' ':9,' 3.  LN(1)=0; LN(-1)=PI*i; LN(+/-i)=+/-0.5*PI*i.');
WRITELN;
WRITELN ('LN(z)':35);
WRITELN ('z':11,'rectangular':27,'EXP( LN(z) ) = z':32);
WRITELN ('     ------------  ---------------------------  ',
'---------------------------');
FOR k := 1 TO 22 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
IF   CAbs(a[k]) = 0.0
THEN WRITELN ('undefined':18)
ELSE BEGIN
CLn (z,a[k]);
CExp (z1,z);
WRITELN (CStr(z,12,9,rectangular),'  ',CStr(z1,12,9,rectangular))
END
END;
WRITELN;
WRITELN;

WRITELN ('Complex exponential:  CExp = EXP(z)');
WRITELN;
WRITELN ('EXP(z)':35);
WRITELN ('z':11,'rectangular':27,'LN( EXP(z) ) = z':32);
WRITELN ('     ------------  ---------------------------  ',
'---------------------------');
FOR k := 1 TO 22 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
CExp (z,a[k]);
CLn (z1,z);
IF   CAbs(z) > 10.0
THEN m := 7
ELSE m := 9;
WRITELN (CStr(z,12,m,rectangular),'  ',CStr(z1,12,m,rectangular))
END;
WRITELN;
WRITELN;

WRITELN ('Complex power:  CPwr = z1^z2');
WRITELN;
WRITELN ('z^(-1+i)':36,'z^(-1+i)':29);
WRITELN ('z':11,'rectangular':27,'polar':26);
WRITELN ('     ------------  ---------------------------  ',
'-----------------------------');
CSet (z1, -1,1, rectangular);
FOR k := 1 TO 21 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
IF   CAbs(a[k]) = 0.0
THEN WRITELN ('undefined':18)
ELSE BEGIN
CPwr (z,a[k],z1);
WRITELN (CStr(z,12,9,rectangular),' ',CStr(z,12,9,polar))
END
END;
WRITELN;
WRITELN;

WRITELN ('Complex cosine:  CCos = COS(z)');
WRITELN;
WRITELN ('COS(z)':35,'COS(z)':29);
WRITELN ('z':11,'rectangular':27,'polar':26);
WRITELN ('     ------------  ---------------------------  ',
'-----------------------------');
FOR k := 1 TO 21 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
CCos (z,a[k]);
CIntPwr (csave[k], z,2);  {save COS^2}
IF   CAbs(z) > 10.0
THEN m := 7
ELSE m := 9;
WRITELN (CStr(z,12,m,rectangular),' ',CStr(z,12,m,polar))
END;
WRITELN;
WRITELN;

WRITELN ('Complex sine:  CSin = SIN(z)');
WRITELN;
WRITELN ('SIN(z)':35);
WRITELN ('z':11,'rectangular':27,'SIN^2(z)+COS^2(z)=1':32);
WRITELN ('     ------------  ---------------------------  ',
'---------------------------');
FOR k := 1 TO 21 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
CSin (z,a[k]);
CIntPwr (z1, z,2);      {SIN^2}
CAdd (z1, z1,csave[k]); {SIN^2 + COS^2}
IF   CAbs(z) > 10.0
THEN m := 7
ELSE m := 9;
WRITELN (CStr(z,12,m,rectangular),'  ',CStr(z1,12,9,rectangular))
END;
WRITELN;
WRITELN;

WRITELN ('Complex tangent:  CTan = TAN(z)');
WRITELN;
WRITELN ('TAN(z)':35,'TAN(z)':29);
WRITELN ('z':11,'rectangular':27,'polar':26);
WRITELN ('     ------------  ---------------------------  ',
'-----------------------------');
FOR k := 1 TO 21 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
CTan (z,a[k]);
IF   CAbs(z) > 10.0
THEN m := 7
ELSE m := 9;
WRITELN (CStr(z,12,m,rectangular),' ',CStr(z,12,m,polar))
END;
WRITELN;
WRITELN;

WRITELN ('Complex hyperbolic cosine:  CCosh = COSH(z)');
WRITELN;
WRITELN ('COSH(z)':36,'COSH(z)':29);
WRITELN ('z':11,'rectangular':27,'polar':26);
WRITELN ('     ------------  ---------------------------  ',
'-----------------------------');
FOR k := 1 TO 21 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
CCosh (z,a[k]);
CIntPwr (csave[k], z,2);  {save COSH^2}
IF   CAbs(z) > 10.0
THEN m := 7
ELSE m := 9;
WRITELN (CStr(z,12,m,rectangular),' ',CStr(z,12,m,polar))
END;
WRITELN;
WRITELN;

WRITELN ('Complex hyperbolic sine:  CSinh = SINH(z)');
WRITELN;
WRITELN ('SINH(z)':36);
WRITELN ('z':11,'rectangular':27,'COSH^2(z)-SINH^2(z)=1':34);
WRITELN ('     ------------  ---------------------------  ',
'---------------------------');
FOR k := 1 TO 21 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
CSinh (z,a[k]);
CIntPwr (z1, z,2);      {SINH^2}
CSub (z1, csave[k],z1); {COSH^2 - SINH^2}
IF   CAbs(z) > 10.0
THEN m := 7
ELSE m := 9;
WRITELN (CStr(z,12,m,rectangular),'  ',CStr(z1,12,9,rectangular))
END;
WRITELN;
WRITELN;

WRITELN ('Complex hyperbolic tangent:  CTanh = TANH(z)');
WRITELN;
WRITELN ('TANH(z)':36,'TANH(z)':29);
WRITELN ('z':11,'rectangular':27,'polar':26);
WRITELN ('     ------------  ---------------------------  ',
'-----------------------------');
FOR k := 1 TO 21 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
CTanh (z,a[k]);
IF   CAbs(z) > 10.0
THEN m := 4
ELSE m := 9;
WRITELN (CStr(z,12,m,rectangular),' ',CStr(z,12,m,polar))
END;
WRITELN;
WRITELN;

WRITELN ('Absolute value of complex number:  CAbs = ABS(z)');
WRITELN;
WRITELN ('z':11,'ABS(z)':17);
WRITELN ('     ------------  ------------');
FOR k := 1 TO 21 DO BEGIN
WRITELN (k:3,' ',CStr(a[k],5,1,rectangular),'  ',CAbs(a[k]):12:9)
END;
WRITELN;

WRITELN ('Complex integer power:  CIntPwr = z^n  ',
'(using DeMoivre''s Theorem)');
WRITELN;
WRITELN ('z^3':34,'z^3':29);
WRITELN ('z':11,'rectangular':27,'polar':26);
WRITELN ('     ------------  ---------------------------  ',
'-----------------------------');
FOR k := 1 TO 21 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
IF   CAbs(a[k]) = 0.0
THEN WRITELN ('undefined':18)
ELSE BEGIN
CIntPwr (z,a[k],3);
IF   CAbs(z) > 10.0
THEN m := 7
ELSE m := 9;
WRITELN (CStr(z,12,m,rectangular),' ',CStr(z,12,m,polar))
END
END;
WRITELN;
WRITELN;

WRITELN ('Complex conjugate:  CConjugate = z*');
WRITELN;
WRITELN ('z*':35,'z*':29);
WRITELN ('z':11,'rectangular':28,'polar':26);
WRITELN ('     ------------  ---------------------------  ',
'-----------------------------');
FOR k := 1 TO 21 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
CConjugate (z,a[k]);
WRITELN (CStr(z,12,8,rectangular),' ',CStr(z,12,8,polar))
END;
WRITELN;
WRITELN;

WRITELN ('Complex square root:  CSqrt = SQRT(z)');
WRITELN;
WRITELN ('SQRT(z)':36,'SQRT(z)':28);
WRITELN ('z':11,'root 1':25,'root 2':28);
WRITELN ('     ------------  ---------------------------  ',
'---------------------------');
FOR k := 1 TO 21 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
CSqrt (z,a[k]);       {same as CRoot (z,a[k],0,2)}
CRoot (z1,a[k],1,2);
WRITELN (CStr(z,12,9,rectangular),'  ',CStr(z1,12,9,rectangular))
END;
WRITELN;
WRITELN;

WRITELN ('The three cube roots of -1+i:  (-1+i)^(1/3)');
WRITELN ('(See Schaum''s Outline Series "Complex Variables", 1964, ',
'p. 18, problem 29.)');
WRITELN;
WRITELN ('z^(1/3)':35,'z^(1/3)':29);
WRITELN ('z':11,'rectangular':27,'polar':26);
WRITELN ('     ------------  ---------------------------  ',
'-----------------------------');
CSet (z1, -1,1, rectangular);
FOR k := 0 TO 2 DO BEGIN
WRITE (k:3,' ',CStr(z1,5,1,rectangular),'  ');
CRoot (z,z1,k,3);
WRITELN (CStr(z,12,9,rectangular),' ',CStr(z,12,9,polar))
END;
WRITELN;
WRITELN;

WRITELN ('Complex Bessel function:  CI0 = I0(z)');
WRITELN;
WRITELN ('I0(z)':36,'I0(z)':29);
WRITELN ('z':11,'rectangular':27,'polar':26);
WRITELN ('     ------------  ---------------------------  ',
'-----------------------------');
FOR k := 1 TO 21 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
CI0 (z,a[k]);
IF   CAbs(z) > 10.0
THEN m := 7
ELSE m := 9;
WRITELN (CStr(z,12,m,rectangular),' ',CStr(z,12,m,polar))
END;
WRITELN;
WRITELN;

WRITELN ('Complex Bessel function:  CJ0 = J0(z)');
WRITELN;
WRITELN ('J0(z)':36,'J0(z)':29);
WRITELN ('z':11,'rectangular':27,'polar':26);
WRITELN ('     ------------  ---------------------------  ',
'-----------------------------');
FOR k := 1 TO 21 DO BEGIN
WRITE (k:3,' ',CStr(a[k],5,1,rectangular),'  ');
CJ0 (z,a[k]);
IF   CAbs(z) > 10.0
THEN m := 7
ELSE m := 9;
WRITELN (CStr(z,12,m,rectangular),' ',CStr(z,12,m,polar))
END;
WRITELN;
WRITELN;

WRITELN ('Removing "Fuzz" from real numbers for zero test:');
WRITELN;  {Note:  CStr calls CConvert that calls CDefuzz}
CSet (z, -3.21E-14,7.65E-14, rectangular);
WRITELN ('  Before:  ',z.x:18:15,' +',z.y:18:15,'i');
CDeFuzz (z);
WRITELN ('   After:  ',CStr(z,18,15,rectangular));
WRITELN;
CSet (z, -3.21E-14,PI, polar);
WRITELN ('  Before:  ',z.r:18:15,'*CIS(',z.theta:18:15,')');
CDeFuzz (z);
WRITELN ('   After:  ',CStr(z,18,15,polar));
WRITELN;
WRITELN;

WRITELN ('Miscellaneous:  FixAngle -- keep angle in interval (-PI..PI)');
WRITELN;

WRITELN ('    -------- --------');
FOR n := -4 TO 8 DO BEGIN
x := n*PI/2.0;
y := FixAngle(x);
WRITELN (n:3,' ',x:8:5,' ',y:8:5)
END;
WRITELN;
WRITELN;

WRITELN ('Real power function:  Pwr = x^y');
WRITELN;
WRITELN ('        x        y         x^y');
WRITELN ('    -------- -------- ------------');
WRITELN (' ':4,2.1:8:5,' ',-2.5:8:5,Pwr(2.1,-2.5):12:9);
WRITELN (' ':4,2.1:8:5,' ', 2.5:8:5,Pwr(2.1, 2.5):12:9);
WRITELN (' ':4,1.4:8:5,' ', 8.9:8:5,Pwr(1.2, 8.9):12:9);
WRITELN (' ':4,0.0:8:5,' ', 2.0:8:5,Pwr(0.0, 2.0):12:9);
WRITELN (' ':4,4.2:8:5,' ', 0.0:8:5,Pwr(4.2, 0.0):12:9);
WRITELN;

END {cdemo}.