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