DECLARE SUB CPSI (X#, Y#, PSR#, PSI#) DECLARE SUB WCoeffs (X#, W#()) DECLARE SUB Coeffs (A#()) DECLARE FUNCTION XExp# (X#) DECLARE FUNCTION SinH# (X#) DECLARE FUNCTION CosH# (X#) DECLARE FUNCTION TanH# (X#) '********************************************************************* '* Calculate the Psi function for a complex argument * '* ----------------------------------------------------------------- * '* EXPLANATION: * '* * '* Purpose: This program computes the psi function psi(z) * '* for a complex argument using subroutine CPSI * '* Input : x --- Real part of z * '* y --- Imaginary part of z * '* Output: PSR --- Real part of psi(z) * '* PSI --- Imaginary part of psi(z) * '* Examples: * '* x y Re[psi(z)] Im[psi(z)] * '* ------------------------------------------- * '* 3.0 2.0 1.16459152 .67080728 * '* 3.0 -2.0 1.16459152 -.67080728 * '* -3.0 2.0 1.39536075 2.62465344 * '* -3.0 -2.0 1.39536075 -2.62465344 * '* * '* NOTE: Psi(z) = Gamma'(z) / Gamma(z) * '* ----------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Please enter x and y (z=x+iy): 3,2 * '* * '* x y Re[Psi(z)] Im[Psi(z)] * '* -------------------------------------------------------- * '* 3 2 1.164591515373977 .6708072826422302 * '* * '* ----------------------------------------------------------------- * '* 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 MCPSI; DEFDBL A-H, O-Z DEFINT I-N DIM aa AS DOUBLE DIM Z AS DOUBLE CLS PRINT INPUT " Please enter x and y (z=x+iy): ", X, Y PRINT PRINT " x y Re[Psi(z)] Im[Psi(z)]" PRINT " ---------------------------------------------------------" CALL CPSI(X, Y, PSR, PSI) PRINT " "; X; " "; Y; " "; PSR; " "; PSI PRINT END 'of main program SUB Coeffs (A()) A(1) = 1.648721270700128# A(2) = 1.284025416687742# A(3) = 1.133148453066826# A(4) = 1.064494458917859# A(5) = 1.031743407499103# A(6) = 1.015747708586686# A(7) = 1.007843097206448# A(8) = 1.003913889338348# A(9) = 1.001955033591003# END SUB ' hyperbolic cosine Function FUNCTION CosH (X) Y = XExp(X) Y = (Y + (1# / Y)) / 2# CosH = Y END FUNCTION SUB CPSI (X, Y, PSR, PSI) ' ============================================= ' Purpose: Compute the psi function for a ' complex argument ' Input : x --- Real part of z ' y --- Imaginary part of z ' Output: PSR --- Real part of psi(z) ' PSI --- Imaginary part of psi(z) ' ============================================= DIM A(8) 'Initialize Table A A(1) = -.08333333333333#: A(2) = 8.333333333333333D-03 A(3) = -3.968253968253968D-03: A(4) = 4.166666666666667D-03 A(5) = -7.575757575757576D-03: A(6) = 2.109279609279609D-02 A(7) = -8.333333333333333D-02: A(8) = .4432598039215686# PI = 4# * ATN(1#) IF Y = 0# AND X = INT(X) AND X <= 0# THEN PSR = 1D+200 PSI = 0# ELSE IF X < 0! THEN X1 = X Y1 = Y X = -X Y = -Y END IF X0 = X IF X < 8! THEN N = 8 - INT(X) X0 = X + 1# * N END IF IF X0 = 0# AND Y <> 0# THEN TH = .5 * PI IF X0 <> 0# THEN TH = ATN(Y / X0) Z2 = X0 * X0 + Y * Y Z0 = SQR(Z2) PSR = LOG(Z0) - .5 * X0 / Z2 PSI = TH + .5 * Y / Z2 FOR K = 1 TO 8 PSR = PSR + A(K) * Z2 ^ (-K) * COS(2# * K * TH) PSI = PSI - A(K) * Z2 ^ (-K) * SIN(2# * K * TH) NEXT K IF X < 8# THEN RR = 0# RI = 0# FOR K = 1 TO N RR = RR + (X0 - K) / ((X0 - K) ^ 2 + Y * Y) RI = RI + Y / ((X0 - K) ^ 2 + Y * Y) NEXT K PSR = PSR - RR PSI = PSI + RI END IF IF X1 < 0# THEN TN = SIN(PI * X) / COS(PI * X) TM = TanH(PI * Y) CT2 = TN * TN + TM * TM PSR = PSR + X / (X * X + Y * Y) + PI * (TN - TN * TM * TM) / CT2 PSI = PSI - Y / (X * X + Y * Y) - PI * TM * (1! + TN * TN) / CT2 X = X1 Y = Y1 END IF END IF END SUB 'CPSI '*---------------------------------------------* '* Hyperbolic sine Function * '* ------------------------------------------- * '* This Procedure uses the definition of the * '* hyperbolic sine and the modified cordic * '* exponential Function XExp(X) to approximate * '* SINH(X) over the entire range of real X. * '* ------------------------------------------- * FUNCTION SinH (X) ' Is X small enough to cause round off error? IF ABS(X) < .35 THEN GOTO 11 ' Calculate SINH(X) using exponential definition ' Get Y:=EXP(X) Y = XExp(X) Y = (Y - (1# / Y)) / 2# GOTO 20 11 'series approximation (for X small) Z = 1#: Y = 1# FOR I = 1 TO 8 Z = Z * X * X / ((2 * I) * (2 * I + 1)) Y = Y + Z NEXT I Y = X * Y 20 SinH = Y END FUNCTION ' hyperbolic tangent Function ' TANH(X]:=SINH(X)/COSH(X) FUNCTION TanH (X) V = SinH(X) Y = CosH(X) TanH = V / Y END FUNCTION SUB WCoeffs (X, W()) aa = .5 Z = X FOR I = 1 TO 9 W(I) = 0# IF Z > aa THEN W(I) = 1# Z = Z - W(I) * aa aa = aa / 2# NEXT I END SUB 'Modified cordic exponential subroutine FUNCTION XExp (X) ' This subroutine takes an input value and returns Y:=EXP(X) ' X may be any positive or negative real value. DIM A(9), W(9) E = EXP(1#) ' Get coefficients CALL Coeffs(A()) ' Reduce the range of X K = INT(X) X = X - 1# * K ' Determine the weighting coeffs, W(I) CALL WCoeffs(X, W()) ' Calculate products Y = 1# FOR I = 1 TO 9 IF W(I) > 0# THEN Y = Y * A(I) NEXT I ' Perform residual multiplication Y = Y * (1# + Z * (1# + Z / 2# * (1# + Z / 3# * (1# + Z / 4#)))) ' Account for factor EXP(K) IF K < 0 THEN E = 1# / E IF ABS(K) < 1 THEN GOTO 10 FOR I = 1 TO ABS(K) Y = Y * E NEXT I ' Restore X X = X + 1# * K 10 XExp = Y END FUNCTION