/* 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;