'**************************************************************
'*       Calculate Function Gamma with a complex argument     *
'* ---------------------------------------------------------- *
'* EXPLANATION:                                               *
'* Purpose: This program computes the gamma function G(z)     *
'*          or Ln[G(z)] for a complex argument using          *
'*          subroutine CGAMA                                  *
'* Input :  x  --- Real part of z                             *
'*          y  --- Imaginary part of z                        *
'*          KF --- Function code                              *
'*          KF=0 for Ln[G(z)]                                 *
'*          KF=1 for G(z)                                     *
'* Output:  GR --- Real part of Ln[G(z)] or G(z)              *
'*          GI --- Imaginary part of Ln[G(z)] or G(z)         *
'* Examples:                                                  *
'*    x         y           Re[G(z)]           Im[G(z)]       *
'*  --------------------------------------------------------  *
'*   2.50      5.00     .2267360319E-01    -.1172284404E-01   *
'*   5.00     10.00     .1327696517E-01     .3639011746E-02   *
'*   2.50     -5.00     .2267360319E-01     .1172284404E-01   *
'*   5.00    -10.00     .1327696517E-01    -.3639011746E-02   *
'*                                                            *
'*    x         y          Re[LnG(z)]         Im[LnG(z)]      *
'*  --------------------------------------------------------  *
'*   2.50      5.00    -.3668103262E+01     .5806009801E+01   *
'*   5.00     10.00    -.4285507444E+01     .1911707090E+02   *
'*   2.50     -5.00    -.3668103262E+01    -.5806009801E+01   *
'*   5.00    -10.00    -.4285507444E+01    -.1911707090E+02   *
'* ---------------------------------------------------------- *
'* SAMPLE RUNS:                                               *
'*                                                            *
'* Please enter KF, x and y: 1,2.5,5                          *
'*                                                            *
'*    x    y          Re[G(z)]             Im[G(z)]           *
'*  --------------------------------------------------------  *
'*   2.5   5    2.267360318980015D-02  -1.172284404171513D-02 *
'*                                                            *
'* Please enter KF, x and y: 0,2.5,5                          *
'*                                                            *
'*    x    y         Re[LnG(z)]           Im[LnG(z)]          *
'*  --------------------------------------------------------  *
'*   2.5   5   -3.66810326235451        5.806009800636287     *
'*                                                            *
'* ---------------------------------------------------------- *
'* REFERENCE:                                                 *
'*    "Fortran Routines for Computation of Special Functions, *
'*     jin.ece.uiuc.edu/routines/routines.html".              *
'*                                                            *
'*                        Basic Release By J-P Moreau, Paris. *
'*                                (www.jpmoreau.fr)           *
'**************************************************************
'PROGRAM MCGAMA
DEFDBL A-H, O-Z
DEFINT I-N

    CLS
    PRINT
    INPUT "  Please enter KF, x and y: ", KF, X, Y
    PRINT
    IF KF = 1 THEN
      PRINT "    x     y           Re[G(z)]           Im[G(z)]"
    ELSE
      PRINT "    x     y          Re[LnG(z)]          Im[LnG(z)]"
    END IF
    PRINT "   ----------------------------------------------------------"

    GOSUB 1000 'call CGAMA(X,Y,KF,GR,GI)

    PRINT "  "; X; "  "; Y; "  "; GR; "  "; GI
    PRINT

END 'of main program

    'auxiliary functions
500 'Function SINH(xx)
    expx = EXP(xx)
    SINH = .5 * (expx - 1# / expx)
    RETURN

600 'Function COSH(xx)
      expx = EXP(xx)
      COSH = .5 * (expx + 1# / expx)
    RETURN


1000 'Subroutine CGAMA(X,Y,KF,GR,GI)
' ===========================================================
'       Purpose: Compute the gamma function G(z) or Ln[G(z)]
'                for a complex argument
'       Input :  x  --- Real part of z
'                y  --- Imaginary part of z
'                KF --- Function code
'                       KF=0 for Ln[G(z)]
'                       KF=1 for G(z)
'       Output:  GR --- Real part of Ln[G(z)] or G(z)
'                GI --- Imaginary part of Ln[G(z)] or G(z)
' ===========================================================
    DIM A(10)
'   G0,GR1,GI1,SR,SI,T,TH,TH1,TH2,X1,Y1,Z1,Z2: Double
'   J,K,NA: Integer

    PI = 4# * ATN(1#)

    'Initialize table A
    A(1) = 8.333333333333333D-02: A(2) = -2.777777777777778D-03
    A(3) = 7.936507936507937D-04: A(4) = -5.952380952380952D-04
    A(5) = 8.417508417508418D-04: A(6) = -1.917526917526918D-03
    A(7) = 6.41025641025641D-03: A(8) = -2.955065359477124D-02
    A(9) = .1796443723688307#: A(10) = -1.3924322169059#

    IF Y = 0! AND X = INT(X) AND X <= 0! THEN
       GR = 1E+30  'arbitrary big number
       GI = 0#
       RETURN
    ELSEIF X < 0! THEN
       X1 = X
       Y1 = Y
       X = -X
       Y = -Y
    END IF
    X0 = X
    IF X <= 7# THEN
      NA = INT(7# - X)
      X0 = X + NA
    END IF
    Z1 = SQR(X0 * X0 + Y * Y)
    TH = ATN(Y / X0)
    GR = (X0 - .5) * LOG(Z1) - TH * Y - X0 + .5 * LOG(2! * PI)
    GI = TH * (X0 - .5) + Y * LOG(Z1) - Y
    FOR K = 1 TO 10
      T = Z1 ^ (1 - 2 * K)
      GR = GR + A(K) * T * COS((2# * K - 1#) * TH)
      GI = GI - A(K) * T * SIN((2# * K - 1#) * TH)
    NEXT K
    IF X <= 7# THEN
      GR1 = 0#
      GI1 = 0#
      FOR J = 0 TO NA - 1
        GR1 = GR1 + .5 * LOG((X + J) ^ 2 + Y * Y)
        GI1 = GI1 + ATN(Y / (X + J))
      NEXT J
      GR = GR - GR1
      GI = GI - GI1
    END IF
    IF X1 < 0# THEN
      Z1 = SQR(X * X + Y * Y)
      TH1 = ATN(Y / X)
      xx = PI * Y: GOSUB 500: GOSUB 600
      SR = -SIN(PI * X) * COSH
      SI = -COS(PI * X) * SINH
      Z2 = SQR(SR * SR + SI * SI)
      TH2 = ATN(SI / SR)
      IF SR < 0# THEN TH2 = PI + TH2
      GR = LOG(PI / (Z1 * Z2)) - GR
      GI = -TH1 - TH2 - GI
      X = X1
      Y = Y1
    END IF
    IF KF = 1 THEN
      G0 = EXP(GR)
      GR = G0 * COS(GI)
      GI = G0 * SIN(GI)
    END IF
    RETURN

'end of file mcgama.bas