'****************************************************************** '* Purpose: This program computes the confluent * '* hypergeometric function M(a,b,x) using * '* subroutine CHGM * '* Input : a --- Parameter * '* b --- Parameter ( b <> 0,-1,-2,... ) * '* x --- Argument * '* Output: HG --- M(a,b,x) * '* Example: * '* a b x M(a,b,x) * '* ----------------------------------------- * '* 1.5 2.0 20.0 .1208527185D+09 * '* 4.5 2.0 20.0 .1103561117D+12 * '* -1.5 2.0 20.0 .1004836854D+05 * '* -4.5 2.0 20.0 -.3936045244D+03 * '* 1.5 2.0 50.0 .8231906643D+21 * '* 4.5 2.0 50.0 .9310512715D+25 * '* -1.5 2.0 50.0 .2998660728D+16 * '* -4.5 2.0 50.0 -.1807604439D+13 * '* -------------------------------------------------------------- * '* 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) * '****************************************************************** DEFDBL A-H, O-Z DEFINT I-N PRINT PRINT " a b x M(a,b,x) " PRINT " ----------------------------------------" A = 1.5: B = 2.0: X = 20.0 CALL CHGM(A, B, X, HG) PRINT " "; A; " "; B; " "; X; " "; HG A = 4.5 CALL CHGM(A, B, X, HG) PRINT " "; A; " "; B; " "; X; " "; HG A = -1.5 CALL CHGM(A, B, X, HG) PRINT " "; A; " "; B; " "; X; " "; HG A = -4.5 CALL CHGM(A, B, X, HG) PRINT " "; A; " "; B; " "; X; " "; HG A = 1.5: X = 50.0 CALL CHGM(A, B, X, HG) PRINT " "; A; " "; B; " "; X; " "; HG A = 4.5 CALL CHGM(A, B, X, HG) PRINT " "; A; " "; B; " "; X; " "; HG A = -1.5 CALL CHGM(A, B, X, HG) PRINT " "; A; " "; B; " "; X; " "; HG A = -4.5 CALL CHGM(A, B, X, HG) PRINT " "; A; " "; B; " "; X; " "; HG END SUB CHGM (A, B, X, HG) ' =================================================== ' Purpose: Compute confluent hypergeometric function ' M(a,b,x) ' Input : a --- Parameter ' b --- Parameter ( b <> 0,-1,-2,... ) ' x --- Argument ' Output: HG --- M(a,b,x) ' Routine called: GAMMA for computing â(x) ' =================================================== PI = 3.141592653589793# A0 = A A1 = A X0 = X HG = 0# IF B = 0# OR B = -ABS(INT(B)) THEN HG = 1D+15 ELSEIF A = 0.0 OR X = 0.0 THEN HG = 1.0 ELSEIF (A = -1.0) THEN HG = 1.0 - X / B ELSEIF (A = B) THEN HG = EXP(X) ELSEIF (A - B = 1.0) THEN HG = (1.0 + X / B) * EXP(X) ELSEIF (A = 1.0 AND B = 2.0) THEN HG = (EXP(X) - 1.0) / X ELSEIF (A = INT(A) AND A < 0.0) THEN M = INT(-A) R = 1.0 HG = 1.0 FOR K = 1 TO M R = R * (A + K - 1.0) / K / (B + K - 1.0) * X HG = HG + R NEXT K END IF IF (HG <> 0.0) GOTO 50 IF (X < 0.0) THEN A = B - A A0 = A X = ABS(X) END IF IF (A < 2.0) THEN NL = 0 IF (A >= 2.0) THEN NL = 1 LA = INT(A) A = A - LA - 1.0 END IF FOR N = 0 TO NL IF (A0 >= 2.0) THEN A = A + 1.0 IF (X <= 30.0 + ABS(B) OR A < 0.0) THEN HG = 1.0 RG = 1.0 FOR J = 1 TO 500 RG = RG * (A + J - 1.0) / (J * (B + J - 1.0)) * X HG = HG + RG IF (ABS(RG / HG) < 1E-15) THEN GOTO 25 NEXT J ELSE CALL GAMMA(A, TA) CALL GAMMA(B, TB) XG = B - A CALL GAMMA(XG, TBA) SUM1 = 1.0 SUM2 = 1.0 R1 = 1.0 R2 = 1.0 FOR I = 1 TO 8 R1 = -R1 * (A + I - 1.0) * (A - B + I) / (X * I) R2 = -R2 * (B - A + I - 1.0) * (A - I) / (X * I) SUM1 = SUM1 + R1 20 SUM2 = SUM2 + R2 NEXT I HG1 = TB / TBA * X ^ (-A) * COS(PI * A) * SUM1 HG2 = TB / TA * EXP(X) * X ^ (A - B) * SUM2 HG = HG1 + HG2 END IF 25 IF (N = 0) THEN Y0 = HG IF (N = 1) THEN Y1 = HG NEXT N IF (A0 >= 2.0) THEN FOR I = 1 TO LA - 1 HG = ((2.0 * A - B + X) * Y1 + (B - A) * Y0) / A Y0 = Y1 Y1 = HG 35 A = A + 1.0 NEXT I END IF IF (X0 < 0.0) THEN HG = HG * EXP(X0) A = A1 X = X0 50 END SUB SUB GAMMA (X, GA) ' ================================================== ' Purpose: Compute gamma function â(x) ' Input : x --- Argument of â(x) ' ( x is not equal to 0,-1,-2,úúú) ' Output: GA --- â(x) ' ================================================== DIM G(27) PI = 3.141592653589793# IF (X = INT(X)) THEN IF (X > 0!) THEN GA = 1! M1 = X - 1 FOR K = 2 TO M1 GA = GA * K NEXT K ELSE GA = 1E+15 END IF ELSE IF (ABS(X) > 1!) THEN Z = ABS(X) M = INT(Z) R = 1! FOR K = 1 TO M R = R * (Z - K) NEXT K Z = Z - M ELSE Z = X END IF G(1) = 1!: G(2) = .5772156649015329# G(3) = -.6558780715202538#: G(4) = -.0420026350340952# G(5) = .1665386113822915#: G(6) = -.0421977345555443# G(7) = -9.621971527877001D-03: G(8) = .007218943246663# G(9) = -.0011651675918591#: G(10) = -.0002152416741149# G(11) = .0001280502823882#: G(12) = -.0000201348547807# G(13) = -.0000012504934821#: G(14) = .000001133027232# G(15) = -.0000002056338417#: G(16) = .000000006116095# G(17) = .0000000050020075#: G(18) = -.0000000011812746# G(19) = .0000000001043427#: G(20) = 7.7823E-12 G(21) = -3.6968E-12: G(22) = 5.1E-13 G(23) = -2.06E-14: G(24) = -5.4E-15 G(25) = 1.4E-13: G(26) = 1E-15 GR = G(26) FOR K = 25 TO 1 STEP -1 GR = GR * Z + G(K) NEXT K GA = 1! / (GR * Z) IF (ABS(X) > 1!) THEN GA = GA * R IF (X < 0!) THEN GA = -PI / (X * GA * SIN(PI * X)) END IF END IF END SUB