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

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

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

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

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

  % True -1
  % False 0

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

  type Complex double x,y

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

  sys ErrCode 
  ' Error code from the latest function evaluation

  Complex C_zero   = {0, 0}
  Complex C_one    = {1, 0}
  Complex C_i      = {0, 1}
  Complex C_MaxNum = {MaxNum, MaxNum}

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

  string operations="
  Cmplx
  Polar
  Cconj
  Cneg
  Cinv
  Csqrt
  Clog
  Cexp
  Croot
  Cpower
  CpolyI
  CpolyR
  CpolyC
  CfracR
  Csin
  Ccos
  Csincos
  Ctan
  CAsin
  CAcos
  CAtan
  Csinh
  Ccosh
  Csincosh
  Ctanh
  CAsinh
  CAcosh
  CAtanh
  CApproxLnGamma
  CLnGamma

  "

  ErrCode = FOk


  'OPERATORS
  ==========
  '
  macro Complex_op(a,b,c)
    '
    macro .move (a,b)
      #if typecodeof(b)=typecodeof(complex)
        a.x=b.x
        a.y=b.y
      #else
        a.x=b
        a.y=b
      #endif
    end macro
    '
    macro .neg (a,b)
      a.x=-b.x
      a.y=-b.y
    end macro
    '
    macro .+ (a,b,c)
      #if typecodeof(c)=typecodeof(complex)
        a.x=b.x+c.x
        a.y=b.y+c.y
      #else
        a.x=b.x+c
        a.y=b.y+c
      #endif
    end macro
    '
    macro .- (a,b,c)
      #if typecodeof(c)=typecodeof(complex)
        a.x=b.x-c.x
        a.y=b.y-c.y
      #else
        a.x=b.x-c
        a.y=b.y-c
      #endif
    end macro
    '
    macro .* (a,b,c)
      #if typecodeof(c)=typecodeof(complex)
        ErrCode = FOk
        a.x = b.x * c.x - b.y * c.y
        a.y = b.x * c.y + b.y * c.x
      #else
        a.x=b.x*c
        a.y=b.y*c
      #endif
    end macro
    '
    macro ./ (a,b,c)
      #if typecodeof(b)=typecodeof(complex)
        #if typecodeof(c)=typecodeof(complex)
          if (b.x = 0) and (b.y = 0)
            ErrCode = FOverflow
            a = C_MaxNum
          else
            double dvs = 1 / (c.x * c.x + c.y * c.y)
            ErrCode = FOk
            a.x=(b.x * c.x + b.y * c.y) * dvs
            a.y=(b.y * c.x - b.x * c.y) * dvs
          endif
        #else 'c is scalar
          if c = 0 then
            ErrCode = FOverflow
            a=C_MaxNum
          else
            ErrCode = FOk
            a.x=b.x/c
            a.y=b.y/c
          endif
        #endif
      #else 'b is scalar, assume c is complex
        if (c.x = 0) and (c.y = 0)
          ErrCode = FOverflow
          a=C_MaxNum
        else
          double Temp
          temp = b / ( (c.x * c.x) + (c.y * c.y) )
          ErrCode = FOk
          a={c.x * Temp, - c.y * Temp}
        endif
      #endif
    end macro
    '
    macro .= (a,b,c)
    if (b.x=c.x)and(b.y=c.y)
      a=1
    else
      a=0
    endif
    end macro
    '
    macro .<> (a,b,c)
    if (b.x=c.x)and(b.y=c.y)
      a=0
    else
      a=1
    endif
    end macro
    '
    macro .^ (a,b,c)
    #if typecodeof(C)<0x40 'integer. assume B is complex
      PowCCI(a,b,c)
    #elseif typecodeof(B)=typecodeof(Complex) 'assume C is scalar
      PowCCS(a,b,c)
    #elseif typecodeof(C)=typecodeof(Complex) 'assume B is scalar
      PowCSC(a,b,c)
    #endif
    end macro
  end macro 'complex_op


  function Cmplx(double x,y) as Complex
  =====================================
  ' Returns the Complex number X + iY
  ErrCode = FOk
  function={x,y}
  end function

  function Polar(double R,Theta) as Complex
  =========================================
  ' Returns the Complex number R * (cos(Theta) + i * sin(Theta))
  ErrCode = FOk
  if r < 0
    ErrCode = FDomain
    function=C_zero
  elseif r = 0
    function=C_zero
  else
    function={ R * cos(Theta), R * sin(Theta) }
  endif
  end function


  function CReal(Complex*z) as double
  ===================================
  ' Real part of Z
  ErrCode = FOk
  return z.x
  end function


  function CImag(Complex*z) as double
  ===================================
  ' Imaginary part of Z
  ErrCode = FOk
  return z.y
  end function

  function CSgn(Complex*z) as double
  ==================================
  ' Complex sign
  ErrCode = FOk
  if z.x > 0
    return 1
  elseif z.x < 0
    return -1
  elseif z.y > 0
    return 1
  elseif z.y < 0
    return -1
  else
    return 0
  endif
  end function

  function cabs(Complex*z) 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)
  '
  double AbsX, AbsY, R
  '
  AbsX = abs(z.x)
  AbsY = abs(z.y)

  if AbsX = AbsY
    return AbsX * Sqrt2
  elseif AbsX > AbsY
    R = AbsY / AbsX
    return AbsX * sqr(1 + R * R)
  else
    R = AbsX / AbsY
    return AbsY * sqr(1 + R * R)
  endif
  end function

  function CAbs2(Complex*z) as double
  ===================================
  ' Squared modulus of Z    
  ErrCode = FOk
  return z.x * z.x + z.y * z.y
  end function

  function CArg(Complex*z) as double
  ==================================
  ' Argument of Z
  ErrCode = FOk
  return atan(z.y, z.x)
  end function


  function CConj(Complex*z) as Complex
  ====================================
  ' Complex conjugate
  ErrCode = FOk
  function={z.x, -z.y}
  end function

  function CNeg(Complex*z) as Complex
  ===================================
  ' Change sign
  ErrCode = FOk
  function={- z.x, - z.y}
  end function


  function CInv(Complex*z) as Complex
  ===================================
  ' Inverse 1/Z
  '
  if z.x = 0 and z.y = 0
    ErrCode = FOverflow
    function=C_MaxNum
  else
    double Temp = 1 / (z.x * z.x + z.y * z.y)
    ErrCode = FOk
    function={z.x * Temp, - z.y * Temp}
  endif
  end function


  function CSqrt(Complex*z) as Complex
  ====================================
  ' Complex square root of Z without destructive underflow or overflow 
  '  cf. Numerical Recipes 
  '
  if z.x = 0 and z.y = 0
    function=C_zero
  else
    double X, Y, U, V, W, R
    X = abs(z.x)
    Y = abs(z.y)
    if X >= Y
      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)))
    endif
    if z.x >= 0.0
      U = W
      V = z.y / (2 * U)
    else
      if z.y >= 0 then V = W else V = - W
      U = z.y / (2 * V)
    endif
    ErrCode = FOk
    function={U, V}
  endif
  end function   

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

  function CLog(Complex*z) as Complex
  ===================================
  ' Principal part of Complex logarithm
  double R, LnR
  ErrCode = FOk
  '
  R = CAbs(z)
  if R < 0
    ErrCode = FDomain
    LnR = - MaxNum
  elseif R = 0
    ErrCode = FSing
    LnR = - MaxNum
  else
    LnR = log(R)
  endif
  function={LnR, CArg(Z)}
  end function  


  function CExp(Complex*z) as Complex
  ===================================
  ' Complex exponential
  double ExpX
  ErrCode = FOk
  if z.x < MinLog
    ErrCode = FUnderflow
    ExpX = 0
  elseif z.x > MaxLog
    ErrCode = FOverflow
    ExpX = MaxNum
  else
    ExpX = exp(z.x)
  endif
  function={ExpX * cos(z.y), ExpX * sin(z.y)}
  end function
  

  function CRoot(Complex*z, int K,N) 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)
    ErrCode = FDomain
    function=C_zero  
  else
    ErrCode = FOk
    function=Polar( pow(CAbs(z), 1/N) , (CArg(z) + K * TwoPi) / N)
  endif
  end function


  sub PowCCI (Complex *R,*A, int N)
  =================================
' Computes A^N by repeated multiplications

  int M
  complex B, Res
  
  ErrCode = FOk

  if N < 0
    M = - N  
   B = 1 / A
    IF ErrCode = FOverflow then R=C_MaxNum
  else
    M = N
    B = A
  endif
  
  ' Legendre's algorithm for minimizing the number of multiplications
  
  Res = C_one
  while M > 0
    if M and 1 then Res = Res * B
    B = B * B
    M = M >>1
  wend
  
  R=Res
  end sub


  sub PowCCS(complex*A,*B,double C)
  =================================
  if frac(B) = 0
    PowCCI(A,B,C)
    exit sub
  endif
  ErrCode = FOk
  if (B.X = 0) and (B.Y = 0)
    if C = 0
      A=C_one
      exit sub
    elseif C > 0
      A=C_zero
      exit sub
    endif
    ErrCode = FSing
    A=C_MaxNum
    exit sub
  endif
  a=Polar(CAbs(B) ^ C, CArg(B) * C)
  end sub


  sub PowCSC(complex*A,double B,complex*C)
  ========================================
  if B < 0 THEN
    ErrCode = FSing
    A=C_zero
    exit sub
  endif
  ErrCode = FOk
  if B = 0
    if C.X = 0 AND C.Y = 0
      A=C_one
      exit sub
    else
      A=C_zero
      exit sub
    endif
  endif
  A=CExp(C * log(B))
  end sub


  function CPower(Complex*V, int N) as Complex
  ============================================
  ' Computes A^N by repeated multiplications
  int M
  Complex B,R
  ErrCode = FOk
  if N < 0 
    M = - N  
    B = 1 / V
    if ErrCode = FOverflow
      function=C_MaxNum
      exit function
    endif
  else
    M = N
    B = V
  endif  
  '
  ' Legendre's algorithm for minimizing the number of multiplications
  '
  R = C_one
  while M > 0
    if (M and 1) then R = R * B
    B = B * B
    M = M >> 1
  wend
  function=R
  end function


  function CPower (Complex *V,*B) as Complex
  ==========================================
  ErrCode = FOk
  if (V.x = 0) and (V.y = 0)
    if (B.x = 0) and (B.y = 0)
      function=C_one   ' lim V^V = 1 when V -> 0
      exit function
    else
      function=C_zero  ' 0^B = 0
    endif  
  endif
  function=Cexp(B * CLog(V))
  end function

  ' ------------------------------------------------------------------
  ' 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
                        
  ' ------------------------------------------------------------------                        


  !* Coef(int Deg) as Complex
                                 
  function CPolyI(Complex*Z, sys pCoef, int Deg) as Complex
  =========================================================
  ' Polynomial with complex coefficients
  ' P(Z) = Coef(0) + Coef(1) * Z + Coef(2) * Z^2 +... + Coef(Deg) * Z^Deg
  int I
  Complex P
  @Coef=pCoef
  P = Coef(Deg)
  for I = Deg - 1 to 0 step -1
    P = P * Z + cmplx(Coef(I),0)
  next
  function=P
  end function

  !* Coef(int Deg) as double

  function CPolyR(Complex*Z, sys pCoef, int Deg) as Complex
  =========================================================
  ' Polynomial with real coefficients

  int I
  Complex P
  @Coef=pCoef
  P.x = Coef(Deg)
  for I = Deg - 1 to 0 step -1
    P = P * Z + cmplx(Coef(I),0)
  next
  function=P
  end function

  !* Coef(int Deg) as Complex

  function CPolyC(Complex*Z, sys pCoef, int Deg1,Deg2) 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
  int I
  @coef=pCoef
  Complex P = Coef(Deg1)
  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
  '
  function=P / Q
  end function 

  !* Coef(int Deg) as double

  function CFracR(Complex*Z, sys pCoef, Deg1,Deg2) as Complex
  ===========================================================
  ' Rational fraction with real coefficients
  '
  sys I
  @Coef=pCoef
  Complex P = Coef(Deg1)
  Complex Q = Cmplx(Coef(Deg1 + Deg2 + 1), 0)
  '
  for I = Deg1 - 1 to 0 step -1
    P = P * Z + Cmplx(Coef(I),0)
  next I
  '
  for I = Deg1 + Deg2 to Deg1 + 1 step -1
    Q = Q * Z + Cmplx(Coef(I),0)
  next I
  '
  function=P / Q
  end function 

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



  sub SinhCosh(double X, *SinhX, *CoshX)
  ======================================
  ' 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
  endif
  '
  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(Complex*z) as Complex
  ===================================
  double SinhY, CoshY
  SinhCosh z.y, SinhY, CoshY
  function={sin(z.x) * CoshY, cos(z.x) * SinhY}
  end function

  function CCos(Complex*z) as Complex
  ===================================
  double SinhY, CoshY
  SinhCosh(z.y, SinhY, CoshY)
  function={cos(z.x) * CoshY, - sin(z.x) * SinhY}
  end function

  sub CSinCos(Complex *z, *SinZ, *CosZ)
  =====================================    
  double SinX, CosX, SinhY, CoshY
  SinX = SIN(z.x)
  CosX = COS(z.x)
  SinhCosh(z.y, SinhY, CoshY)
  SinZ = {SinX * CoshY, CosX * SinhY}
  CosZ = {CosX * CoshY, - SinX * SinhY}
  end sub    
    
  function CTan(Complex*z) as Complex
  ===================================
  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
    function={SinX2 / Temp, SinhY2 / Temp}
  else
    ' Z = Pi/2 + k*Pi 
    ErrCode = FSing  
    function={MaxNum, 0} : return
  endif
  end function

  function CASin(Complex*z) as Complex
  ====================================
  double  Rp, Rm, S, T, X2, XX, YY
  Complex Z1
  '
  Z1 = {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
  function={asin(Rp - Rm), CSgn(Z1) * log(T + sqr(T * T - 1))}
  end function

  function CACos(Complex*z) as Complex
  ====================================
  ' ArcCos(Z) = Pi/2 - ArcSin(Z)
  Complex A = CASin(Z)
  ErrCode = FOk
  function={PiDiv2 - A.X, - A.Y}
  end function

  function CATan(Complex*z) as Complex
  ====================================
  if (z.x = 0) and (abs(z.y) = 1) THEN  ' Z = +/- i 
    ErrCode = FSing
    function={0, SGN(z.y) * MaxNum}
    exit function
  endif

  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
  
  function={0.5 * (atan(z.x, - Ym1) - atan(- z.x, Yp1)),
    0.25 * log((XX + Yp1 * Yp1) / (XX + Ym1 * Ym1))}
  end function

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

  function CSinh(Complex*z) as Complex
  ====================================
  double SinhX, CoshX, SinY, CosY
  SinY = SIN(z.y)
  CosY = COS(z.y)
  SinhCosh(z.x, SinhX, CoshX)
  ErrCode = FOk
  function={SinhX * CosY, CoshX * SinY}
  end function


  function CCosh(Complex*z) as Complex
  ====================================
  double SinhX, CoshX, SinY, CosY
  SinY = sin(z.y)
  CosY = cos(z.y)
  SinhCosh( z.x,SinhX, CoshX)
  ErrCode = FOk
  function={CoshX * CosY, SinhX * SinY}
  end function

  sub CSinhCosh(Complex*z, *SinhZ, *CoshZ)
  ========================================
  double SinY, CosY, SinhX, CoshX
  SinY = sin(z.y)
  CosY = cos(z.y)
  SinhCosh(z.x, SinhX, CoshX)
  SinhZ = {SinhX * CosY, CoshX * SinY}
  CoshZ = {CoshX * CosY, SinhX * SinY}
  end sub    

  function CTanh(Complex*z) as Complex
  ====================================
  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  ' Z = i * (Pi/2 + k*Pi) 
    ErrCode = FSing
    function={0, MaxNum}
  else
    ErrCode = FOk   
    function={SinhX2 / Temp, SinY2 / Temp}
  endif
  end function

  function CASinh(Complex*z) as Complex
  =====================================
  ' ArcSinh(Z) = -i*ArcSin(i*Z)
  '
  Complex iZ = {- z.y, z.x}
  Complex b = CASin(iZ)
  ErrCode = FOk
  function={b.y, - b.x}
  end function

  function CACosh(Complex*z) as Complex
  =====================================
  ' ArcCosh(Z) = CSgn(Y + i(1-X)) * i * ArcCos(Z) where Z = X+iY
  '
  Complex b = {z.y, 1 - z.x}  ' Y + i(1-X)
  double  S = CSgn(b)
  b = CACos(Z)
  ErrCode = FOk
  function={- S * b.y, S * b.x}
  end function

  function CATanh(Complex*z) as Complex
  =====================================
  ' ArcTanh(Z) = -i*ArcTan(i*Z) }  
  '
  if (abs(z.x) = 1) and (z.y = 0)  ' Z = +/- 1
    ErrCode = FSing
    function={sgn(z.x) * MaxNum, 0}
  else
    Complex iZ = {- z.y, z.x}
    Complex b = CATan(iZ)
    ErrCode = FOk
    function={b.y, - b.x}
  endif
  end function

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


  function CApproxLnGamma(Complex*Z) 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) 
  '
  indexbase 1
  double C[] = {
  8.33333333333333E-02, -2.77777777777778E-03, _
  7.93650793650794E-04, -5.95238095238095E-04, _
  8.41750841750842E-04, -1.91752691752692E-03, _
  6.41025641025641E-03, -2.95506535947712E-02}
  
  Complex Powers[8], Temp, Sum, tmp
  int I
  Powers[1] = 1 / Z
  Temp = Powers[1] * Powers[1]
  for I = 2 to 8
    Powers[I] = Powers[I - 1] * Temp
  next  
  Sum = (Z - 0.5 ) * CLog(Z) - Z + Ln2PiDiv2
  for I = 8 to 1 step -1
    Sum = Sum + C[i] * Powers[i]
  next I
  function=Sum
  end function

  function CLnGamma(Complex*z) as Complex
  =======================================
  ErrCode = FOk
  
  ' Negative integer
  if z.x <= 0 and frac(z.x) = 0 and z.y = 0
    ErrCode = FSing
    function=C_MaxNum
    exit function
  endif
  '
  ' 3rd or 4th quadrant
  if z.y < 0
    function=CConj(CLnGamma(CConj(z)))
    exit function
  endif
  ' 
  ' "left" of NBS table range
  if z.x < 9 THEN
    function=CLnGamma(z+1) - CLog(z)
    exit function
  endif
  '
  ' NBS table range: 9 < Re(z) < 10
  function=CApproxLnGamma(z)
  end function    
