' ******************************************************************
' Complex number library
' ------------------------------------------------------------------
' Based on ComplexMath Delphi library by E. F. Glynn
' http://www.efg2.com/Lab/Mathematics/Complex/index.html
' ******************************************************************

' ------------------------------------------------------------------
' Mathematical constants
' ------------------------------------------------------------------

#DEFINE MaxNum    1.797693134862315E+308  ' Max. floating point number: 2^1024
#DEFINE MaxLog  709.782712893384          ' Max. argument for Exp
#DEFINE MinLog -708.3964185322641         ' Min. argument for Exp
#DEFINE TwoPi     6.283185307179586       ' 2*Pi
#DEFINE PiDiv2    1.570796326794897       ' Pi/2
#DEFINE Ln2PiDiv2 0.9189385332046728      ' Ln(2*Pi)/2
#DEFINE Sqrt2     1.414213562373095#      ' Sqrt(2)

' ------------------------------------------------------------------
' Error codes for function evaluations
' ------------------------------------------------------------------

#DEFINE FOk         0  ' No error
#DEFINE FDomain    -1  ' Argument domain error
#DEFINE FSing      -2  ' Singularity
#DEFINE FOverflow  -3  ' Overflow range error
#DEFINE FUnderflow -4  ' Underflow range error

' ------------------------------------------------------------------
' Boolean constants
' ------------------------------------------------------------------

#DEFINE True -1
#DEFINE False 0

' ------------------------------------------------------------------
' Type definition
' ------------------------------------------------------------------

TYPE Complex
  X AS DOUBLE
  Y AS DOUBLE
END TYPE

' ------------------------------------------------------------------
' Global variables
' ------------------------------------------------------------------

COMMON SHARED AS INTEGER ErrCode 
' Error code from the latest function evaluation

DIM SHARED AS Complex C_zero   = (0, 0)
DIM SHARED AS Complex C_one    = (1, 0)
DIM SHARED AS Complex C_i      = (0, 1)     
DIM SHARED AS Complex C_MaxNum = (MaxNum, MaxNum)

' ------------------------------------------------------------------
' General functions
' ------------------------------------------------------------------

FUNCTION Cmplx(BYVAL X AS DOUBLE, BYVAL Y AS DOUBLE) AS Complex
' Returns the Complex number X + iY

  ErrCode = FOk
  RETURN TYPE<Complex> (X, Y)
END FUNCTION

FUNCTION Polar(BYVAL R AS DOUBLE, BYVAL Theta AS DOUBLE) AS Complex
' Returns the Complex number R * (cos(Theta) + i * sin(Theta))

  ErrCode = FOk
  
  IF R < 0 THEN
    ErrCode = FDomain
    RETURN C_zero
  ELSEIF R = 0 THEN
    RETURN C_zero
  END IF
 
  RETURN TYPE<Complex> (R * COS(Theta), R * SIN(Theta))
END FUNCTION

FUNCTION CReal(BYREF Z AS Complex) AS DOUBLE
' Real part of Z

  ErrCode = FOk
  RETURN Z.X
END FUNCTION

FUNCTION CImag(BYREF Z AS Complex) AS DOUBLE
' Imaginary part of Z

  ErrCode = FOk
  RETURN Z.Y
END FUNCTION

FUNCTION CSgn(BYREF Z AS Complex) AS INTEGER
' Complex sign

  ErrCode = FOk
  IF Z.X > 0 THEN
    RETURN 1
  ELSEIF Z.X < 0 THEN
    RETURN -1
  ELSEIF Z.Y > 0 THEN
    RETURN 1
  ELSEIF Z.Y < 0 THEN
    RETURN -1
  ELSE
    RETURN 0
  END IF
END FUNCTION

FUNCTION CAbs(BYREF Z AS Complex) AS DOUBLE
' Modulus of Z

  ErrCode = FOk
  
  IF Z.X = 0 THEN RETURN ABS(Z.Y)
  IF Z.Y = 0 THEN RETURN ABS(Z.X)

  DIM AS DOUBLE AbsX, AbsY, R

  AbsX = ABS(Z.X)
  AbsY = ABS(Z.Y)

  IF AbsX = AbsY THEN
    RETURN AbsX * Sqrt2
  ELSEIF AbsX > AbsY THEN
    R = AbsY / AbsX
    RETURN AbsX * SQR(1 + R * R)
  ELSE
    R = AbsX / AbsY
    RETURN AbsY * SQR(1 + R * R)
  END IF
END FUNCTION

FUNCTION CAbs2(BYREF Z AS Complex) AS DOUBLE
' Squared modulus of Z    
    
  ErrCode = FOk
  RETURN Z.X * Z.X + Z.Y * Z.Y
END FUNCTION

FUNCTION CArg(BYREF Z AS Complex) AS DOUBLE
' Argument of Z

  ErrCode = FOk
  RETURN ATAN2(Z.Y, Z.X)
END FUNCTION

FUNCTION CConj(BYREF Z AS Complex) AS Complex
' Complex conjugate
  
  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X, - Z.Y)    
END FUNCTION

FUNCTION CNeg(BYREF Z AS Complex) AS Complex
' Change sign
  
  ErrCode = FOk
  RETURN TYPE<Complex> (- Z.X, - Z.Y)    
END FUNCTION

FUNCTION CInv(BYREF Z AS Complex) AS Complex
' Inverse 1/Z
  
  IF Z.X = 0 AND Z.Y = 0 THEN
    ErrCode = FOverflow
    RETURN C_MaxNum
  END IF

  DIM AS DOUBLE Temp = 1 / (Z.X * Z.X + Z.Y * Z.Y)
  
  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X * Temp, - Z.Y * Temp)
END FUNCTION

FUNCTION CSqrt(BYREF Z AS Complex) AS Complex
' Complex square root of Z without destructive underflow or overflow 
'  cf. Numerical Recipes 

  IF Z.X = 0 AND Z.Y = 0 THEN RETURN C_zero
 
  DIM AS DOUBLE X, Y, U, V, W, R
  
  X = ABS(Z.X)
  Y = ABS(Z.Y)
  
  IF X >= Y THEN
    R = Y / X  
    W = SQR(X) * SQR(0.5 * (1 + SQR(1 + R * R)))
  ELSE 
    R = X / Y
    W = SQR(Y) * SQR(0.5 * (R + SQR(1 + R * R)))
  END IF
    
  IF Z.X >= 0.0 THEN
    U = W
    V = Z.Y / (2 * U)
  ELSE
    IF Z.Y >= 0 THEN V = W ELSE V = - W
    U = Z.Y / (2 * V)
  END IF
  
  ErrCode = FOk
  RETURN TYPE<Complex> (U, V)
END FUNCTION   

' ------------------------------------------------------------------
' Logarithm, Exponential, Roots
' ------------------------------------------------------------------

FUNCTION CLog(BYREF Z AS Complex) AS Complex
' Principal part of Complex logarithm

  DIM AS DOUBLE R, LnR

  ErrCode = FOk
  
  R = CAbs(Z)
  
  IF R < 0 THEN
    ErrCode = FDomain
    LnR = - MaxNum
  ELSEIF R = 0 THEN
    ErrCode = FSing
    LnR = - MaxNum
  ELSE
    LnR = LOG(R)
  END IF
  
  RETURN TYPE<Complex> (LnR, CArg(Z))
END FUNCTION  

FUNCTION CExp(BYREF Z AS Complex) AS Complex
' Complex exponential

  DIM AS DOUBLE ExpX

  ErrCode = FOk

  IF Z.X < MinLog THEN
    ErrCode = FUnderflow
    ExpX = 0
  ELSEIF Z.X > MaxLog THEN
    ErrCode = FOverflow
    ExpX = MaxNum
  ELSE
    ExpX = EXP(Z.X)
  END IF
  
  RETURN TYPE<Complex> (ExpX * COS(Z.Y), ExpX * SIN(Z.Y))
END FUNCTION  

FUNCTION CRoot(BYREF Z AS Complex, BYVAL K AS INTEGER, BYVAL N AS INTEGER) AS Complex
' CRoot can calculate all 'N' roots of 'Z' by varying 'K' from 0 to N-1
' This is an application of DeMoivre's theorem.

  IF N <= 0 OR K < 0 OR K >= N THEN
    ErrCode = FDomain
    RETURN C_zero  
  END IF

  ErrCode = FOk
  RETURN Polar(CAbs(Z) ^ (1 / N), (CArg(Z) + K * TwoPi) / N)
END FUNCTION

' ------------------------------------------------------------------
' Operators
' ------------------------------------------------------------------

OPERATOR = (BYREF Z1 AS Complex, BYREF Z2 AS Complex) AS INTEGER
  ErrCode = FOk
  RETURN (Z1.X = Z2.X AND Z1.Y = Z2.Y) 
END OPERATOR

OPERATOR <> (BYREF Z1 AS Complex, BYREF Z2 AS Complex) AS INTEGER
  ErrCode = FOk
  RETURN (Z1.X <> Z2.X OR Z1.Y <> Z2.Y) 
END OPERATOR

OPERATOR + (BYREF Z1 AS Complex, BYREF Z2 AS Complex) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (Z1.X + Z2.X, Z1.Y + Z2.Y) 
END OPERATOR

OPERATOR + (BYVAL A AS DOUBLE, BYREF Z AS Complex) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X + A, Z.Y) 
END OPERATOR

OPERATOR + (BYREF Z AS Complex, BYVAL A AS DOUBLE) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X + A, Z.Y) 
END OPERATOR

OPERATOR - (BYREF Z1 AS Complex, BYREF Z2 AS Complex) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (Z1.X - Z2.X, Z1.Y - Z2.Y) 
END OPERATOR

OPERATOR - (BYVAL A AS DOUBLE, BYREF Z AS Complex) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (A - Z.X, - Z.Y) 
END OPERATOR

OPERATOR - (BYREF Z AS Complex, BYVAL A AS DOUBLE) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X - A, Z.Y) 
END OPERATOR

OPERATOR - (BYREF Z AS Complex) AS Complex
' Unary minus

  ErrCode = FOk
  RETURN TYPE<Complex> (- Z.X, - Z.Y) 
END OPERATOR

OPERATOR * (BYREF Z1 AS Complex, BYREF Z2 AS Complex) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (Z1.X * Z2.X - Z1.Y * Z2.Y, _
                        Z1.X * Z2.Y + Z1.Y * Z2.X) 
END OPERATOR

OPERATOR * (BYVAL A AS DOUBLE, BYREF Z AS Complex) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (A * Z.X, A * Z.Y) 
END OPERATOR

OPERATOR * (BYREF Z AS Complex, BYVAL A AS DOUBLE) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (A * Z.X, A * Z.Y)
END OPERATOR

OPERATOR / (BYREF Z1 AS Complex, BYREF Z2 AS Complex) AS Complex

  IF Z2.X = 0 AND Z2.Y = 0 THEN
    ErrCode = FOverflow
    RETURN C_MaxNum
  END IF

  DIM AS DOUBLE Temp = 1 / (Z2.X * Z2.X + Z2.Y * Z2.Y)
  
  ErrCode = FOk
  RETURN TYPE<Complex> ((Z1.X * Z2.X + Z1.Y * Z2.Y) * Temp, _
                        (Z1.Y * Z2.X - Z1.X * Z2.Y) * Temp)
END OPERATOR

OPERATOR / (BYVAL A AS DOUBLE, BYREF Z AS Complex) AS Complex

  IF Z.X = 0 AND Z.Y = 0 THEN
    ErrCode = FOverflow
    RETURN C_MaxNum
  END IF

  DIM AS DOUBLE Temp = A / (Z.X * Z.X + Z.Y * Z.Y)
  
  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X * Temp, - Z.Y * Temp)
END OPERATOR

OPERATOR / (BYREF Z AS Complex, BYVAL A AS DOUBLE) AS Complex

  IF A = 0 THEN
    ErrCode = FOverflow
    RETURN C_MaxNum
  END IF

  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X / A, Z.Y / A)
END OPERATOR

OPERATOR ^ (BYREF A AS Complex, BYVAL N AS INTEGER) AS Complex
' Computes A^N by repeated multiplications

  DIM AS INTEGER M
  DIM AS Complex B, Res
  
  ErrCode = FOk

  IF N < 0 THEN 
    M = - N  
    B = 1 / A
    IF ErrCode = FOverflow THEN RETURN C_MaxNum
  ELSE
    M = N
    B = A
  END IF  
  
  ' Legendre's algorithm for minimizing the number of multiplications
  
  Res = C_one
  DO WHILE M > 0
    IF (M AND 1) THEN Res = Res * B
    B = B * B
    M = M SHR 1
  LOOP
  
  RETURN Res
END OPERATOR

OPERATOR ^ (BYREF A AS Complex, BYREF B AS Complex) AS Complex

  ErrCode = FOk
  
  IF A.X = 0 AND A.Y = 0 THEN
    IF B.X = 0 AND B.Y = 0 THEN
      RETURN C_one   ' lim A^A = 1 when A -> 0
    ELSE
      RETURN C_zero  ' 0^B = 0
    END IF  
  END IF
  
  RETURN CExp(B * CLog(A))
END OPERATOR

OPERATOR ^ (BYREF A AS Complex, BYVAL B AS DOUBLE) AS Complex

  IF FRAC(B) = 0 THEN RETURN A ^ CINT(B)

  ErrCode = FOk
  
  IF A.X = 0 AND A.Y = 0 THEN
    IF B = 0 THEN
      RETURN C_one
    ELSEIF B > 0 THEN
      RETURN C_zero
    ELSE
      ErrCode = FSing
      RETURN C_MaxNum
    END IF
  END IF

  RETURN Polar(CAbs(A) ^ B, CArg(A) * B)
END OPERATOR

OPERATOR ^ (BYVAL A AS DOUBLE, BYREF B AS Complex) AS Complex
  
  IF A < 0 THEN
    ErrCode = FSing
    RETURN C_zero
  END IF
  
  ErrCode = FOk
  
  IF A = 0 THEN
    IF B.X = 0 AND B.Y = 0 THEN
      RETURN C_one
    ELSE
      RETURN C_zero
    END IF
  END IF
    
  RETURN CExp(B * LOG(A))
END OPERATOR

' ------------------------------------------------------------------
' Polynomials and rational fractions
' ------------------------------------------------------------------

DECLARE FUNCTION CPoly OVERLOAD (BYREF Z   AS Complex, _ 
                                 Coef()    AS Complex, _ 
                                 BYVAL Deg AS INTEGER) AS Complex

DECLARE FUNCTION CPoly (BYREF Z   AS Complex, _ 
                        Coef()    AS DOUBLE,  _ 
                        BYVAL Deg AS INTEGER) AS Complex
                       
DECLARE FUNCTION CFrac OVERLOAD (BYREF Z    AS Complex, _
                                 Coef()     AS Complex, _ 
                                 BYVAL Deg1 AS INTEGER, _ 
                                 BYVAL Deg2 AS INTEGER) AS Complex
                                 
DECLARE FUNCTION CFrac (BYREF Z    AS Complex, _
                        Coef()     AS DOUBLE, _ 
                        BYVAL Deg1 AS INTEGER, _ 
                        BYVAL Deg2 AS INTEGER) AS Complex
                        
' ------------------------------------------------------------------                        
                                 
FUNCTION CPoly(BYREF Z   AS Complex, _ 
               Coef()    AS Complex, _ 
               BYVAL Deg AS INTEGER) AS Complex

' Polynomial with complex coefficients
' P(Z) = Coef(0) + Coef(1) * Z + Coef(2) * Z^2 +... + Coef(Deg) * Z^Deg

  DIM AS INTEGER I
  DIM AS Complex P

  P = Coef(Deg)
  FOR I = Deg - 1 TO 0 STEP -1
    P = P * Z + Coef(I)
  NEXT I
  
  RETURN P
END FUNCTION

FUNCTION CPoly(BYREF Z   AS Complex, _ 
               Coef()    AS DOUBLE,  _ 
               BYVAL Deg AS INTEGER) AS Complex

' Polynomial with real coefficients

  DIM AS INTEGER I
  DIM AS Complex P = Cmplx(Coef(Deg), 0)
  
  FOR I = Deg - 1 TO 0 STEP -1
    P = P * Z + Coef(I)
  NEXT I
  
  RETURN P
END FUNCTION

FUNCTION CFrac(BYREF Z    AS Complex, _
               Coef()     AS Complex, _ 
               BYVAL Deg1 AS INTEGER, _ 
               BYVAL Deg2 AS INTEGER) AS Complex

' Rational fraction with complex coefficients
'
'               Coef(0) + Coef(1) * Z + ... + Coef(Deg1) * Z^Deg1
' F(Z) = --------------------------------------------------------------------
'        Coef(Deg1 + 1) + Coef(Deg1+2) * Z + ... + Coef(Deg1+Deg2+1) * Z^Deg2


  DIM AS INTEGER I
  DIM AS Complex P = Coef(Deg1)
  DIM AS Complex Q = Coef(Deg1 + Deg2 + 1)

  FOR I = Deg1 - 1 TO 0 STEP -1
    P = P * Z + Coef(I)
  NEXT I

  FOR I = Deg1 + Deg2 TO Deg1 + 1 STEP -1
    Q = Q * Z + Coef(I)
  NEXT I
  
  RETURN P / Q
END FUNCTION 

FUNCTION CFrac(BYREF Z    AS Complex, _
               Coef()     AS DOUBLE,  _ 
               BYVAL Deg1 AS INTEGER, _ 
               BYVAL Deg2 AS INTEGER) AS Complex

' Rational fraction with real coefficients

  DIM AS INTEGER I
  DIM AS Complex P = Cmplx(Coef(Deg1), 0)
  DIM AS Complex Q = Cmplx(Coef(Deg1 + Deg2 + 1), 0)

  FOR I = Deg1 - 1 TO 0 STEP -1
    P = P * Z + Coef(I)
  NEXT I

  FOR I = Deg1 + Deg2 TO Deg1 + 1 STEP -1
    Q = Q * Z + Coef(I)
  NEXT I
  
  RETURN P / Q
END FUNCTION 

' ------------------------------------------------------------------
' Trigonometric functions
' ------------------------------------------------------------------

SUB SinhCosh(BYVAL X AS DOUBLE, BYREF SinhX AS DOUBLE, BYREF CoshX AS DOUBLE)
' Compute Sinh(X) and Cosh(X) with only one exponential     
  
  IF X < MinLog OR X > MaxLog THEN
    ErrCode = FOverflow
    CoshX = MaxNum
    SinhX = SGN(X) * CoshX
    EXIT SUB
  END IF
  
  DIM AS DOUBLE ExpX, ExpMinusX
  
  ErrCode = FOk
  ExpX = EXP(X)
  ExpMinusX = 1 / ExpX
  SinhX = 0.5 * (ExpX - ExpMinusX)
  CoshX = 0.5 * (ExpX + ExpMinusX)
END SUB  

FUNCTION CSin(BYREF Z AS Complex) AS Complex
  
  DIM AS DOUBLE SinhY, CoshY
  
  SinhCosh Z.Y, SinhY, CoshY
  RETURN TYPE<Complex> (SIN(Z.X) * CoshY, COS(Z.X) * SinhY)
END FUNCTION

FUNCTION CCos(BYREF Z AS Complex) AS Complex
  
  DIM AS DOUBLE SinhY, CoshY
  
  SinhCosh Z.Y, SinhY, CoshY
  RETURN TYPE<Complex> (COS(Z.X) * CoshY, - SIN(Z.X) * SinhY)
END FUNCTION

SUB CSinCos(BYREF Z AS Complex, BYREF SinZ AS Complex, BYREF CosZ AS Complex)
    
  DIM AS DOUBLE SinX, CosX, SinhY, CoshY
  
  SinX = SIN(Z.X)
  CosX = COS(Z.X)
  
  SinhCosh Z.Y, SinhY, CoshY
  
  SinZ = TYPE<Complex> (SinX * CoshY, CosX * SinhY)
  CosZ = TYPE<Complex> (CosX * CoshY, - SinX * SinhY)
END SUB    
    
FUNCTION CTan(BYREF Z AS Complex) AS Complex

  DIM AS DOUBLE X2, Y2, SinX2, CosX2, SinhY2, CoshY2, Temp
  
  X2 = 2 * Z.X
  Y2 = 2 * Z.Y
  
  SinX2 = SIN(X2)
  CosX2 = COS(X2)
  
  SinhCosh Y2, SinhY2, CoshY2
  
  IF ErrCode = FOk THEN Temp = CosX2 + CoshY2 ELSE Temp = CoshY2
  
  IF Temp <> 0 THEN
    RETURN TYPE<Complex> (SinX2 / Temp, SinhY2 / Temp)
  ELSE
    ' Z = Pi/2 + k*Pi 
    ErrCode = FSing  
    RETURN TYPE<Complex> (MaxNum, 0)
  END IF
END FUNCTION

FUNCTION CASin(BYREF Z AS Complex) AS Complex
    
  DIM AS DOUBLE  Rp, Rm, S, T, X2, XX, YY
  DIM AS Complex Z1
  
  Z1 = TYPE<Complex> (Z.Y, - Z.X)
  
  X2 = 2 * Z.X
  XX = Z.X * Z.X
  YY = Z.Y * Z.Y
  S  = XX + YY + 1
  Rp = 0.5 * SQR(S + X2)
  Rm = 0.5 * SQR(S - X2)
  T  = Rp + Rm
   
  ErrCode = FOk
  
  RETURN TYPE<Complex> (ASIN(Rp - Rm), _
                        CSgn(Z1) * LOG(T + SQR(T * T - 1))) 
END FUNCTION

FUNCTION CACos(BYREF Z AS Complex) AS Complex
' ArcCos(Z) = Pi/2 - ArcSin(Z)

  DIM AS Complex A = CASin(Z)
  
  ErrCode = FOk
  RETURN TYPE<Complex> (PiDiv2 - A.X, - A.Y)
END FUNCTION

FUNCTION CATan(BYREF Z AS Complex) AS Complex
    
  IF (Z.X = 0) AND (ABS(Z.Y) = 1) THEN  ' Z = +/- i 
    ErrCode = FSing
    RETURN TYPE<Complex> (0, SGN(Z.Y) * MaxNum)
  END IF

  DIM AS DOUBLE XX, YY, Yp1, Ym1
  
  XX  = Z.X * Z.X
  YY  = Z.Y * Z.Y
  Yp1 = Z.Y + 1
  Ym1 = Z.Y - 1

  ErrCode = FOk
  
  RETURN TYPE<Complex> (0.5 * (ATAN2(Z.X, - Ym1) - ATAN2(- Z.X, Yp1)), _
                        0.25 * LOG((XX + Yp1 * Yp1) / (XX + Ym1 * Ym1)))  
END FUNCTION

' ------------------------------------------------------------------
' Hyperbolic functions
' ------------------------------------------------------------------

FUNCTION CSinh(BYREF Z AS Complex) AS Complex
  
  DIM AS DOUBLE SinhX, CoshX, SinY, CosY
  
  SinY = SIN(Z.Y)
  CosY = COS(Z.Y)
  
  SinhCosh Z.X, SinhX, CoshX
  
  ErrCode = FOk
  RETURN TYPE<Complex> (SinhX * CosY, CoshX * SinY) 
END FUNCTION

FUNCTION CCosh(BYREF Z AS Complex) AS Complex
    
  DIM AS DOUBLE SinhX, CoshX, SinY, CosY
  
  SinY = SIN(Z.Y)
  CosY = COS(Z.Y)
  
  SinhCosh Z.X, SinhX, CoshX
  
  ErrCode = FOk
  RETURN TYPE<Complex> (CoshX * CosY, SinhX * SinY)  
END FUNCTION

SUB CSinhCosh(BYREF Z AS Complex, BYREF SinhZ AS Complex, BYREF CoshZ AS Complex)
    
  DIM AS DOUBLE SinY, CosY, SinhX, CoshX
  
  SinY = SIN(Z.Y)
  CosY = COS(Z.Y)
  
  SinhCosh Z.X, SinhX, CoshX
  
  SinhZ = TYPE<Complex> (SinhX * CosY, CoshX * SinY)
  CoshZ = TYPE<Complex> (CoshX * CosY, SinhX * SinY)
END SUB    

FUNCTION CTanh(BYREF Z AS Complex) AS Complex

  DIM AS DOUBLE X2, Y2, SinY2, CosY2, SinhX2, CoshX2, Temp
  
  X2 = 2.0 * Z.X
  Y2 = 2.0 * Z.Y
  
  SinY2 = SIN(Y2)
  CosY2 = COS(Y2)
  
  SinhCosh X2, SinhX2, CoshX2
  
  IF ErrCode = FOk THEN Temp = CoshX2 + CosY2 ELSE Temp = CoshX2

  IF Temp = 0 THEN  ' Z = i * (Pi/2 + k*Pi) 
    ErrCode = FSing
    RETURN TYPE<Complex> (0, MaxNum)
  END IF  

  ErrCode = FOk   
  RETURN TYPE<Complex> (SinhX2 / Temp, SinY2 / Temp)
END FUNCTION

FUNCTION CASinh(BYREF Z AS Complex) AS Complex
' ArcSinh(Z) = -i*ArcSin(i*Z)

  DIM AS Complex iZ = (- Z.Y, Z.X)
  DIM AS Complex A = CASin(iZ)
  
  ErrCode = FOk
  RETURN TYPE<Complex> (A.Y, - A.X)
END FUNCTION

FUNCTION CACosh(BYREF Z AS Complex) AS Complex
' ArcCosh(Z) = CSgn(Y + i(1-X)) * i * ArcCos(Z) where Z = X+iY

  DIM AS Complex A = (Z.Y, 1 - Z.X)  ' Y + i(1-X)
  DIM AS DOUBLE  S = CSgn(A)
  
  A = CACos(Z)
  
  ErrCode = FOk
  RETURN TYPE<Complex> (- S * A.Y, S * A.X)
END FUNCTION

FUNCTION CATanh(BYREF Z AS Complex) AS Complex
' ArcTanh(Z) = -i*ArcTan(i*Z) }  

  IF ABS(Z.X) = 1 AND Z.Y = 0 THEN  ' Z = +/- 1
    ErrCode = FSing
    RETURN TYPE<complex> (SGN(Z.X) * MaxNum, 0)
  END IF

  DIM AS Complex iZ = (- Z.Y, Z.X)
  DIM AS Complex A = CATan(iZ)
  
  ErrCode = FOk
  RETURN TYPE<complex> (A.Y, - A.X)
END FUNCTION

' ------------------------------------------------------------------
' Gamma function
' ------------------------------------------------------------------

FUNCTION CApproxLnGamma(BYREF Z AS Complex) AS 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) 

  DIM AS DOUBLE C(1 TO 8) = _
    {8.33333333333333E-02, -2.77777777777778E-03, _
     7.93650793650794E-04, -5.95238095238095E-04, _
     8.41750841750842E-04, -1.91752691752692E-03, _
     6.41025641025641E-03, -2.95506535947712E-02}
  
  DIM AS Complex Powers(1 TO 8), Temp, Sum
  
  DIM AS INTEGER I
  
  Powers(1) = 1 / Z
  
  Temp = Powers(1) * Powers(1)
  
  FOR I = 2 TO 8
    Powers(I) = Powers(I - 1) * Temp
  NEXT I  
  
  Sum = (Z - 0.5) * CLog(Z) - Z + Ln2PiDiv2 
  
  FOR I = 8 TO 1 STEP -1
    Sum = Sum + C(i) * Powers(i)
  NEXT I
  
  RETURN Sum
END FUNCTION

FUNCTION CLnGamma(BYREF Z AS Complex) AS Complex
  
  ErrCode = FOk
  
  ' Negative integer
  IF Z.X <= 0 AND FRAC(Z.X) = 0 AND Z.Y = 0 THEN  
    ErrCode = FSing
    RETURN C_MaxNum
  END IF
  
  ' 3rd or 4th quadrant
  IF Z.Y < 0 THEN                     
    RETURN CConj(CLnGamma(CConj(Z))) 
  END IF  

  ' "left" of NBS table range
  IF Z.X < 9 THEN  
    RETURN CLnGamma(Z + 1) - CLog(Z)
  END IF

  ' NBS table range: 9 < Re(z) < 10 
  RETURN CApproxLnGamma(Z)  
END FUNCTION    

' ------------------------------------------------------------------
' Printing functions
' ------------------------------------------------------------------

SUB CPrint(BYREF Z AS Complex, BYREF Mask AS STRING = "####.####")
' Prints a complex in rectangular form

  PRINT USING Mask; Z.X;
  IF Z.Y >= 0 THEN PRINT " + "; ELSE PRINT " - ";
  PRINT USING Mask; ABS(Z.Y); 
  PRINT " * i";
END SUB

SUB CPrintPolar(BYREF Z AS Complex, BYREF Mask AS STRING = "####.####")
' Prints a complex in polar form: Z = |Z| * exp(Arg(Z) * i)

  PRINT USING Mask; CAbs(Z); 
  PRINT " * exp(";
  PRINT USING Mask; CArg(Z); 
  PRINT " * i)"
END SUB


