/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:59 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "slgama.h" #include #include float /*FUNCTION*/ slgama( float x) { long int _l0, i, j; float corr, del, lomega, res, slgama_v, t1, xden, xm1, xm2, xm4, xnum, y, y1, y2, ysq; static float eps, frtbig, xlbig; static float xinf = 0.e0; static float one = 1.0e0; static float half = 0.5e0; static float twelve = 12.0e0; static float zero = 0.0e0; static float four = 4.e0; static float five = 5.e0; static float thrhal = 1.5e0; static float sqrtpi = 0.9189385332046727417803297e0; static float p65 = 0.65e0; static float t5 = 0.76551212348464539649e0; static float aln4 = 1.3862943611198906188e0; static float d2 = 4.227843350984671393993777e-1; static float p2[8]={4.974607845568932035012064e0,5.424138599891070494101986e2, 1.550693864978364947665077e4,1.847932904445632425417223e5,1.088204769468828767498470e6, 3.338152967987029735917223e6,5.106661678927352456275255e6,3.074109054850539556250927e6}; static float q2[8]={1.830328399370592604055942e2,7.765049321445005871323047e3, 1.331903827966074194402448e5,1.136705821321969608938755e6,5.267964117437946917577538e6, 1.346701454311101692290052e7,1.782736530353274213975932e7,9.533095591844353613395747e6}; static float d4 = 1.791759469228055000094023e0; static float p4[8]={1.474502166059939948905062e4,2.426813369486704502836312e6, 1.214755574045093227939592e8,2.663432449630976949898078e9,2.940378956634553899906876e10, 1.702665737765398868392998e11,4.926125793377430887588120e11,5.606251856223951465078242e11}; static float q4[8]={2.690530175870899333379843e3,6.393885654300092398984238e5, 4.135599930241388052042842e7,1.120872109616147941376570e9,1.488613728678813811542398e10, 1.016803586272438228077304e11,3.417476345507377132798597e11,4.463158187419713286462081e11}; static float 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 */ float *const C = &c[0] - 1; float *const P2 = &p2[0] - 1; float *const P4 = &p4[0] - 1; float *const Q2 = &q2[0] - 1; float *const Q4 = &q4[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. *>> 1996-03-30 SLGAMA Krogh Added external statement. *>> 1995-11-16 SLGAMA Krogh SFTRAN => Fortran 77 *>> 1995-10-24 SLGAMA Krogh Removed blanks in numbers for C conversion. *>> 1994-10-19 SLGAMA Krogh Changes to use M77CON *>> 1992-05-19 SLGAMA CLL Corrected computation of FRTBIG. *>> 1991-10-21 SLGAMA CLL Eliminated DGAM1 as a separate subroutine. *>> 1991-01-16 SLGAMA Lawson Using DGAM1 in place of D2MACH/R2MACH *>> 1986-03-18 SLGAMA Lawson Initial code. *--S replaces "?": ?LGAMA, ?GAMMA, ?RAT1, ?ERM1, ?ERV1 * * Designed and programmed by W.J.CODY, Argonne National Lab, 1982. * Minor changes for use in the JPL MATH77 library by C.L.LAWSON & * S.CHAN, JPL, 1983. * 1992-05-19 CLL. Noted that FRTBIG was being computed using XLBIG * before XLBIG was computed. Thus branch on Y .GT. FRTBIG was * unreliable and was likely to cause wrong formulas to be used for * X between 12 and the correct value of FRTBIG. Corrected this. * ---------------------------------------------------------------------- * This routine calculated the LOG(GAMMA) function for a double * precision argument X. Computation is based on an algorithm * outlined in references 1 and 2. The program uses rational * functions that approximate LOG(GAMMA) to least 18 * significant decimal digits. Approximations for X .LT. 12.0 * are unpublished. Lower order approximations can be * substituted on machines with less precise arithmetic. * * Explanaton of machine-dependent constants * * XINF The largest machine representable floating-point * number. * * EPS The smallest positive floating-point number such * that 1.0 + EPS .GT. 1.0 * * 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.) * * FRTBIG The fourth root of XLBIG. * *-- Begin mask code changes * 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 Messages * 1) X .LE. 0., setting result large. * 2) X too large., setting result large. * * ------------------------------------------------------------------ */ /* T5 = LN(SQRT(4*PI)) - HALF */ /* ALN4 = LN(4) */ /* ---------------------------------------------------------------------- * NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX * APPROXIMATION OVER (1.5,4.0). * ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- * NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX * APPROXIMATION OVER (4.0,12.0) * ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- * COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12,INF). * ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */ if (xinf == zero) { xinf = FLT_MAX; eps = FLT_EPSILON; /* Compute XLBIG * * XLBIG must satisfy LogGamma(XLBIG) = 0.875 * XINF. * Solve this equation by successive substitution using the * approximation: log(gamma(x)) ~ x*(log(x) - 1). * Scale this equation by substituting y * 0.875 * XINF for x, * and solve the scaled equation for y. y will range from * 0.0121 when XINF = 2**127 to 0.000176 when XINF = 2**8191. * */ lomega = logf( 0.875e0*xinf ); y1 = 0.01e0; for (j = 1; j <= 7; j++) { y2 = 1.0e0/(lomega + logf( y1 ) - 1.0e0); del = y2 - y1; if (fabsf( del ) < 0.5e-5*y2) goto L_20; y1 = y2; } L_20: ; xlbig = 0.875e0*xinf*y2; frtbig = sqrtf( sqrtf( xlbig ) ); } y = x; if (y <= zero) { serm1( "SLGAMA", 2, 0, "X .LE. 0.,SETTING RESULT LARGE", "X" , x, '.' ); goto L_700; } if (y > xlbig) { serm1( "SLGAMA", 2, 0, "X TOO LARGE,SETTING RESULT LARGE", "X", x, ',' ); serv1( "LIMIT FOR X", xlbig, '.' ); goto L_700; } if (y > twelve) goto L_400; if (y > five) { res = logf( sgamma( y ) ); goto L_900; } if (y > four) goto L_300; if (y > thrhal) goto L_200; if (y >= p65) { xm1 = (y - half) - half; res = srat1( xm1 ); } else if (y > half) { /* *** Here for .5 < Y < .65 we use the formula * LGAM(Y) = LGAM(2Y) - LGAM(Y+ 1/2) - Y*LOG(4) + LOG(SQRT(4*PI)) * */ xm1 = y - half; t1 = -srat1( xm1 ); xm1 += xm1; res = (((t1 + srat1( xm1 )) - y*aln4) + half) + t5; } else if (y > eps) { xm1 = y; res = srat1( xm1 ) - logf( y ); } else { res = -logf( y ); } goto L_900; /* ---------------------------------------------------------------------- * 1.5 .LT. X .LE. 4.0 * ---------------------------------------------------------------------- */ L_200: ; xm2 = (y - one) - one; xden = one; xnum = zero; for (i = 1; i <= 8; i++) { xnum = xnum*xm2 + P2[i]; xden = xden*xm2 + Q2[i]; } res = xm2*(d2 + xm2*(xnum/xden)); goto L_900; /* --------------------------------------------------------------------- * 4.0 .LT. X .LE. 5.0 * --------------------------------------------------------------------- */ L_300: ; xm4 = y - four; xden = -one; xnum = zero; for (i = 1; i <= 8; i++) { xnum = xnum*xm4 + P4[i]; xden = xden*xm4 + Q4[i]; } res = d4 + xm4*(xnum/xden); goto L_900; /* --------------------------------------------------------------------- * EVALUATE FOR ARGUMENT .GE. 12.0 * --------------------------------------------------------------------- */ L_400: ; res = zero; if (y > frtbig) goto L_460; res = C[7]; ysq = y*y; for (i = 1; i <= 6; i++) { res = res/ysq + C[i]; } L_460: ; res /= y; corr = logf( y ); res += sqrtpi - half*corr; res += y*(corr - one); goto L_900; /* --------------------------------------------------------------------- * RETURN FOR BAD ARGUMENTS * --------------------------------------------------------------------- */ L_700: ; res = xinf; /* --------------------------------------------------------------------- * FINAL ADJUSTMENTS AND RETURN * --------------------------------------------------------------------- */ L_900: ; slgama_v = res; return( slgama_v ); } /* end of function */ /* ================================================================== */ float /*FUNCTION*/ srat1( float xm1) { long int i; float srat1_v, xden, xnum; static float d1 = -5.772156649015328605195174e-1; static float p1[8]={4.945235359296727046734888e0,2.018112620856775083915565e2, 2.290838373831346393026739e3,1.131967205903380828685045e4,2.855724635671635335736389e4, 3.848496228443793359990269e4,2.637748787624195437963534e4,7.225813979700288197698961e3}; static float q1[8]={6.748212550303777196073036e1,1.113332393857199323513008e3, 7.738757056935398733233834e3,2.763987074403340708898585e4,5.499310206226157329794414e4, 6.161122180066002127833352e4,3.635127591501940507276287e4,8.785536302431013170870835e3}; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const P1 = &p1[0] - 1; float *const Q1 = &q1[0] - 1; /* end of OFFSET VECTORS */ /* Evaluate the rational function RAT1(XM1) which * with XM1 = X - 1 approximates LOG(GAMMA(X)) * for .5 .LE. X .LE. 1.5 . This rational function * has poor error amplification properties for * .5 .le. X .le. 0.65 so we use it for * 0.65 .le. X .le. 1.5. * * D1, P1(), and Q1() are coefficients for a rational minimax * approximation over (0.5, 1.5). * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ xden = 1.0e0; xnum = 0.0e0; for (i = 1; i <= 8; i++) { xnum = xnum*xm1 + P1[i]; xden = xden*xm1 + Q1[i]; } srat1_v = xm1*(d1 + xm1*(xnum/xden)); return( srat1_v ); } /* end of function */