/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:41 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dgamma.h" #include #include /* PARAMETER translations */ #define C1 0.9189385332046727417803297e0 #define HALF 0.5e0 #define ONE 1.0e0 #define PI 3.1415926535897932384626434e0 #define TWELVE 12.0e0 #define TWO 2.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ double /*FUNCTION*/ dgamma( double x) { LOGICAL32 parity; long int _l0, _l1, i, j, n; double const_, del, dgamma_v, f, fact, fp, res, sum, temp, x1, x2, xden, xnum, y, y1, ysq, z; static double eps, xgbig, xminin; static double xinf = 0.0e0; static double p[8]={-1.71618513886549492533811e0,2.47656508055759199108314e1, -3.79804256470945635097577e2,6.29331155312818442661052e2,8.66966202790413211295064e2, -3.14512729688483675254357e4,-3.61444134186911729807069e4,6.64561438202405440627855e4}; static double q[8]={-3.08402300119738975254353e1,3.15350626979604161529144e2, -1.01515636749021914166146e3,-3.10777167157231109440444e3,2.25381184209801510330112e4, 4.75584627752788110767815e3,-1.34659959864969306392456e5,-1.15132259675553483497211e5}; static double c[7]={-1.910444077728e-03,8.4171387781295e-04,-5.952379913043012e-04, 7.93650793500350248e-04,-2.777777777777681622553e-03,8.333333333333333331554247e-02, 5.7083835261e-03}; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const C = &c[0] - 1; double *const P = &p[0] - 1; double *const Q = &q[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2001-05-25 DGAMMA Krogh Minor change for making .f90 version. *>> 1996-03-30 DGAMMA Krogh Added external statement. *>> 1994-10-20 DGAMMA Krogh Changes to use M77CON *>> 1991-10-21 DGAMMA CLL Eliminated DGAM1 as a separate subroutine. *>> 1991-01-16 DGAMMA Lawson Replaced D2MACH/R2MACH with DGAM1. *>> 1985-08-02 DGAMMA Lawson Initial code. *--D replaces "?": ?GAMMA, ?ERM1, ?ERV1 * * ---------------------------------------------------------------------- * * THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A double precision * ARGUMENT X. PERMITS NEGATIVE AS WELL AS POSITIVE X. NOTE * THAT THE GAMMA FUNCTION HAS POLES AT ZERO AND AT NEGATIVE * ARGUMENTS. COMPUTATION IS BASED ON AN ALGORITHN OUTLINED IN * W.J.CODY, 'AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL * FUNCTIONS', LECTURE NOTES IN MATHEMATICS, 506, NUMERICAL ANALYSIS * DUNDEE, 1975, G. A. WATSON (ED.),SPRINGER VERLAG, BERLIN, 1976. * THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA * FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS * FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED. * THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM HART, ET. AL., * COMPUTER APPROXIMATIONS, WILEY AND SONS, NEW YORK, 1968. * LOWER ORDER APPROXIMATIONS CAN BE SUBSTITUTED FOR THESE ON * MACHINES WITH LESS PRECISE ARITHMETIC. * * Designed & programmed by W.J.CODY, Argonne National Lab.,1982. * Minor changes for the JPL library by C.L.LAWSON & S.CHAN,JPL,1983. * ************************************************************************ * *-- Begin mask code changes * EXPLANATION OF MACHINE-DEPENDENT CONSTANTS * * EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT * 1.0 + EPS .GT. 1.0 * XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER. * XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT * both XMININ and 1/XMININ are representable. * XGBIG - A value such that Gamma(XGBIG) = 0.875 * XINF. * (Computed and used in [D/S]GAMMA.) * XLBIG - A value such that LogGamma(XLBIG) = 0.875 * XINF. * (Computed and used in [D/S]LGAMA.) * * Values of XINF, XGBIG, and XLBIG for some machines: * * XINF XGBIG XLBIG Machines * * 2**127 = 0.170e39 34.81 0.180e37 Vax SP & DP; Unisys SP * 2**128 = 0.340e39 35.00 0.358e37 IEEE SP * 2**252 = 0.723e76 57.54 0.376e74 IBM30xx DP * 2**1023 = 0.899e308 171.46 0.112e306 Unisys DP * 2**1024 = 0.180e309 171.60 0.2216e306 IEEE DP * 2**1070 = 0.126e323 177.78 0.1501e320 CDC/7600 SP * 2**8191 = 0.550e2466 966.94 0.8464e2462 Cray SP & DP *-- End mask code changes * ************************************************************************ * * ERROR RETURNS * * THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR * WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED * TO BE FREE OF UNDERFLOW AND OVERFLOW. * * AUTHOR:W. J. CODY * APPLIED MATHMATICS DIVISION * ARGONNE NATIONAL LABORATORY * ARGONNE, IL 60439 * * LATEST MODIFICATION by Cody: MAY 18, 1982 * * ------------------------------------------------------------------ */ /* C1 = LOG base e of SQRT(2*PI) * */ /* ---------------------------------------------------------------------- * NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX * APPROXIMATION OVER (1,2). * ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- * COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF). * ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- * */ if (xinf == ZERO) { eps = DBL_EPSILON; xinf = DBL_MAX; if (DBL_MIN*DBL_MAX >= ONE) { xminin = DBL_MIN; } else { xminin = ONE/DBL_MAX; } /* Compute XGBIG * * XGBIG will satisfy Gamma(XGBIG) = 0.875 * XINF. * Use a Newton iteration and the following approximation for * the gamma function: * log(gamma(x)) ~ (x - .5)*log(x) - x + 0.5 * log(2 * PI) * */ temp = log( 0.875e0*xinf ); const_ = HALF*log( TWO*PI ) - temp; x1 = temp*0.34e0; for (j = 1; j <= 7; j++) { f = (x1 - HALF)*log( x1 ) - x1 + const_; fp = ((x1 - HALF)/x1) + log( x1 ) - ONE; del = -f/fp; x2 = x1 + del; if (fabs( del ) < 0.5e-5*x2) goto L_45; x1 = x2; } L_45: ; xgbig = x2; } parity = FALSE; fact = ONE; n = 0; y = x; if (y > ZERO) goto L_200; /* ---------------------------------------------------------------------- * ARGUMENT IS NEGATIVE OR ZERO * ---------------------------------------------------------------------- */ y = -x; j = (long)( y ); res = y - (double)( j ); if (res == ZERO) goto L_700; if (j != (j/2)*2) parity = TRUE; fact = -PI/sin( PI*res ); y += ONE; /* ---------------------------------------------------------------------- * ARGUMENT IS POSITIVE * ---------------------------------------------------------------------- */ L_200: if (y < eps) goto L_650; if (y >= TWELVE) goto L_300; y1 = y; if (y >= ONE) goto L_210; /* ---------------------------------------------------------------------- * 0.0 .LT. ARGUMENT .LT. 1.0 * ---------------------------------------------------------------------- */ z = y; y += ONE; goto L_250; /* ---------------------------------------------------------------------- * 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY * ---------------------------------------------------------------------- */ L_210: n = (long)( y ) - 1; y -= (double)( n ); z = y - ONE; /* ---------------------------------------------------------------------- * EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0 * ---------------------------------------------------------------------- */ L_250: xnum = ZERO; xden = ONE; for (i = 1; i <= 8; i++) { xnum = (xnum + P[i])*z; xden = xden*z + Q[i]; } res = (xnum/xden + HALF) + HALF; if (y == y1) goto L_900; if (y1 > y) goto L_280; /* ---------------------------------------------------------------------- * ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0 * ---------------------------------------------------------------------- */ res /= y1; goto L_900; /* ---------------------------------------------------------------------- * ADJUST RESULT FOR CASE 2.0 .LT. 12.0 * ---------------------------------------------------------------------- */ L_280: for (i = 1; i <= n; i++) { res *= y; y += ONE; } goto L_900; /* ---------------------------------------------------------------------- * EVALUATE FOR ARGUMENT .GE. 12.0 * ---------------------------------------------------------------------- */ L_300: if (y > xgbig) goto L_720; ysq = y*y; sum = C[7]; for (i = 1; i <= 6; i++) { sum = sum/ysq + C[i]; } sum = ((sum/y + C1) - y) + (y - HALF)*log( y ); res = exp( sum ); goto L_900; /* ---------------------------------------------------------------------- * ARGUMENT .LT. EPS * ---------------------------------------------------------------------- */ L_650: if (y < xminin) goto L_740; res = ONE/y; goto L_900; /* ---------------------------------------------------------------------- * RETURN FOR SINGULARITIES,EXTREME ARGUMENTS, ETC. * ---------------------------------------------------------------------- */ L_700: derm1( "DGAMMA", 1, 0, "POLE AT 0 AND NEG INTEGERS", "X", x, '.' ); goto L_780; L_720: derm1( "DGAMMA", 2, 0, "X SO LARGE VALUE WOULD OVERFLOW", "X", x, ',' ); derv1( "LIMIT", xgbig, '.' ); goto L_780; L_740: derm1( "DGAMMA", 3, 0, "X TOO NEAR TO A SINGULARITY.VALUE WOULD OVERFLOW." , "X", x, '.' ); L_780: dgamma_v = xinf; goto L_950; /* ---------------------------------------------------------------------- * FINAL ADJUSTMENTS AND RETURN * ---------------------------------------------------------------------- */ L_900: if (parity) res = -res; if (fact != ONE) res = fact/res; dgamma_v = res; L_950: return( dgamma_v ); } /* end of function */