/* Copyright (c) 1995 by R. A. Vowels, from "Introduction to PL/I, Algorithms, and */ /* Structured Programming". Permission is given to reproduce and to use these procedures */ /* as part of a program, and to include them as part of a larger work to be sold for profit. */ /* However, the user is not permitted to sell the procedures separately. Provided always */ /* that these procedures and this copyright notice are reproduced in full. */ /* This function procedure implements the Gamma function. See C. Lanczos, "A Precision */ /* Approximation of the Gamma Function", SIAM Journal on Numerical Analysis, Series B, */ /* Vol. 1, 1964, pp. 86-96. Accuracy: when z is an integer, 1 <= z <= 31, is 18 */ /* significant digits, otherwise better than 1 part in 2*10**(-10). */ /* For speed, the function uses table lookup for Gamma(1) through Gamma(31), and Lanczos' */ /* formula for other values. */ GAMMA: PROCEDURE (Z) OPTIONS (REORDER) RECURSIVE RETURNS ( FLOAT DECIMAL (15)); DECLARE Z FLOAT DECIMAL (15); %INCLUDE MACFACT; DECLARE Y FLOAT DECIMAL (15); %DECLARE N FIXED; %N = 30; DECLARE FACTORIALS (0:N) FLOAT DECIMAL (15) STATIC INITIAL ( FACT (N) ); DECLARE G FLOAT DECIMAL (15) STATIC INITIAL (5); DECLARE PI FLOAT DECIMAL (15) STATIC INITIAL (3.14159265358979324E0); IF Z = 0 THEN SIGNAL ERROR; IF Z < 1 THEN RETURN (GAMMA(Z+1) / Z); IF (Z >= 1) & (Z <= N+1) THEN DO; IF Z = TRUNC(Z) THEN RETURN (FACTORIALS (Z-1)); END; %DEACTIVATE N; Y = SQRT (2*PI) *(Z + G - 0.5)**(Z - 0.5) / EXP (Z + G - 0.5) * (1.000000000178E0 + 76.180091729406E0 / Z - 86.505320327112E0/(Z+1) + 24.014098222230E0 / (Z + 2) - 1.231739516140E0 / (Z + 3) + 0.001208580030E0 / (Z + 4) - 0.000005363820E0 / (Z + 5) ); RETURN ( Y ); END GAMMA; /* This function procedure implements the LOGGAMMA function, which forms the natural */ /* logarithm of the Gamma function. See C. Lanczos, "A Precision Approximation of the */ /* Gamma Function", SIAM Journal on Numerical Analysis, Series B, Vol. 1, 1964, pp. 86-96. */ /* Accuracy: when z is an integer, 1 <= z <= 31, is 18 significant digits, otherwise better */ /* than 1 part in 2*10**(-10). */ /* For speed, the function uses table lookup for LogGamma(1) through LogGamma(31), */ /* and Lanczos' formula for other values. */ LOGGAMMA: PROCEDURE (Z) OPTIONS (REORDER) RETURNS (FLOAT DECIMAL (15)); DECLARE Z FLOAT DECIMAL (15); DECLARE Y FLOAT DECIMAL (15); %DECLARE N FIXED; %N = 30; DECLARE FACTORIALS (0:N) FLOAT DECIMAL (15) STATIC INITIAL ( FACT (N) ); DECLARE G FLOAT DECIMAL (15) STATIC INITIAL (5); DECLARE PI FLOAT DECIMAL (15) STATIC INITIAL (3.14159265358979324E0); IF Z = 0 THEN SIGNAL ERROR; IF Z < 1 THEN RETURN (LOG (GAMMA(Z+1) / Z)); IF (Z >= 1) & (Z <= N+1) THEN DO; IF Z = TRUNC(Z) THEN RETURN (LOG (FACTORIALS (Z-1)) ); END; %DEACTIVATE N; Y = (Z - 0.5)*LOG (Z + G - 0.5) - (Z + G - 0.5) + LOG (SQRT (2*PI)* (1.000000000178E0 + 76.180091729406E0 / Z - 86.505320327112E0/(Z+1) + 24.014098222230E0 / (Z + 2) - 1.231739516140 / (Z + 3) + 0.001208580030E0 / (Z + 4) - 0.000005363820E0 / (Z + 5) ) ); RETURN ( Y ); END LOGGAMMA;