/* gamma.c * * Gamma function * * * * SYNOPSIS: * * double x, y, gamma(); * extern int sgngam; * * y = gammafn( x ); * * * * DESCRIPTION: * * Returns gamma function of the argument. The result is * correctly signed, and the sign (+1 or -1) is also * returned in a global (extern) variable named sgngam. * This variable is also filled in by the logarithmic gamma * function lgam(). * * Arguments |x| <= 34 are reduced by recurrence and the function * approximated by a rational function of degree 6/7 in the * interval (2,3). Large arguments are handled by Stirling's * formula. Large negative arguments are made positive using * a reflection formula. * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * DEC -34, 34 10000 1.3e-16 2.5e-17 * IEEE -170,-33 20000 2.3e-15 3.3e-16 * IEEE -33, 33 20000 9.4e-16 2.2e-16 * IEEE 33, 171.6 20000 2.3e-15 3.2e-16 * * Error for arguments outside the test range will be larger * owing to error amplification by the exponential function. * */ /* lgam() * * Natural logarithm of gamma function * * * * SYNOPSIS: * * double x, y, lgam(); * extern int sgngam; * * y = lgam( x ); * * * * DESCRIPTION: * * Returns the base e (2.718...) logarithm of the absolute * value of the gamma function of the argument. * The sign (+1 or -1) of the gamma function is returned in a * global (extern) variable named sgngam. * * For arguments greater than 13, the logarithm of the gamma * function is approximated by the logarithmic version of * Stirling's formula using a polynomial approximation of * degree 4. Arguments between -33 and +33 are reduced by * recurrence to the interval [2,3] of a rational approximation. * The cosecant reflection formula is employed for arguments * less than -33. * * Arguments greater than MAXLGM return MAXNUM and an error * message. MAXLGM = 2.035093e36 for DEC * arithmetic or 2.556348e305 for IEEE arithmetic. * * * * ACCURACY: * * * arithmetic domain # trials peak rms * DEC 0, 3 7000 5.2e-17 1.3e-17 * DEC 2.718, 2.035e36 5000 3.9e-17 9.9e-18 * IEEE 0, 3 28000 5.4e-16 1.1e-16 * IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17 * The error criterion was relative when the function magnitude * was greater than one but absolute when it was less than one. * * The following test used the relative error criterion, though * at certain points the relative error could be much higher than * indicated. * IEEE -200, -4 10000 4.8e-16 1.3e-16 * */ /* gamma.c */ /* gamma function */ /* Cephes Math Library Release 2.2: July, 1992 Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */ #include "mconf.h" #include #ifdef UNK static double P[] = { 1.60119522476751861407E-4, 1.19135147006586384913E-3, 1.04213797561761569935E-2, 4.76367800457137231464E-2, 2.07448227648435975150E-1, 4.94214826801497100753E-1, 9.99999999999999996796E-1 }; static double Q[] = { -2.31581873324120129819E-5, 5.39605580493303397842E-4, -4.45641913851797240494E-3, 1.18139785222060435552E-2, 3.58236398605498653373E-2, -2.34591795718243348568E-1, 7.14304917030273074085E-2, 1.00000000000000000320E0 }; #define MAXGAM 171.624376956302725 static double LOGPI = 1.14472988584940017414; #endif #ifdef DEC static unsigned short P[] = { 0035047, 0162701, 0146301, 0005234, 0035634, 0023437, 0032065, 0176530, 0036452, 0137157, 0047330, 0122574, 0037103, 0017310, 0143041, 0017232, 0037524, 0066516, 0162563, 0164605, 0037775, 0004671, 0146237, 0014222, 0040200, 0000000, 0000000, 0000000 }; static unsigned short Q[] = { 0134302, 0041724, 0020006, 0116565, 0035415, 0072121, 0044251, 0025634, 0136222, 0003447, 0035205, 0121114, 0036501, 0107552, 0154335, 0104271, 0037022, 0135717, 0014776, 0171471, 0137560, 0034324, 0165024, 0037021, 0037222, 0045046, 0047151, 0161213, 0040200, 0000000, 0000000, 0000000 }; #define MAXGAM 34.84425627277176174 static unsigned short LPI[4] = { 0040222, 0103202, 0043475, 0006750, }; #define LOGPI *(double *)LPI #endif #ifdef IBMPC static unsigned short P[] = { 0x2153, 0x3998, 0xfcb8, 0x3f24, 0xbfab, 0xe686, 0x84e3, 0x3f53, 0x14b0, 0xe9db, 0x57cd, 0x3f85, 0x23d3, 0x18c4, 0x63d9, 0x3fa8, 0x7d31, 0xdcae, 0x8da9, 0x3fca, 0xe312, 0x3993, 0xa137, 0x3fdf, 0x0000, 0x0000, 0x0000, 0x3ff0 }; static unsigned short Q[] = { 0xd3af, 0x8400, 0x487a, 0xbef8, 0x2573, 0x2915, 0xae8a, 0x3f41, 0xb44a, 0xe750, 0x40e4, 0xbf72, 0xb117, 0x5b1b, 0x31ed, 0x3f88, 0xde67, 0xe33f, 0x5779, 0x3fa2, 0x87c2, 0x9d42, 0x071a, 0xbfce, 0x3c51, 0xc9cd, 0x4944, 0x3fb2, 0x0000, 0x0000, 0x0000, 0x3ff0 }; #define MAXGAM 171.624376956302725 static unsigned short LPI[4] = { 0xa1bd, 0x48e7, 0x50d0, 0x3ff2, }; #define LOGPI *(double *)LPI #endif #ifdef MIEEE static unsigned short P[] = { 0x3f24, 0xfcb8, 0x3998, 0x2153, 0x3f53, 0x84e3, 0xe686, 0xbfab, 0x3f85, 0x57cd, 0xe9db, 0x14b0, 0x3fa8, 0x63d9, 0x18c4, 0x23d3, 0x3fca, 0x8da9, 0xdcae, 0x7d31, 0x3fdf, 0xa137, 0x3993, 0xe312, 0x3ff0, 0x0000, 0x0000, 0x0000 }; static unsigned short Q[] = { 0xbef8, 0x487a, 0x8400, 0xd3af, 0x3f41, 0xae8a, 0x2915, 0x2573, 0xbf72, 0x40e4, 0xe750, 0xb44a, 0x3f88, 0x31ed, 0x5b1b, 0xb117, 0x3fa2, 0x5779, 0xe33f, 0xde67, 0xbfce, 0x071a, 0x9d42, 0x87c2, 0x3fb2, 0x4944, 0xc9cd, 0x3c51, 0x3ff0, 0x0000, 0x0000, 0x0000 }; #define MAXGAM 171.624376956302725 static unsigned short LPI[4] = { 0x3ff2, 0x50d0, 0x48e7, 0xa1bd, }; #define LOGPI *(double *)LPI #endif /* Stirling's formula for the gamma function */ #if UNK static double STIR[5] = { 7.87311395793093628397E-4, -2.29549961613378126380E-4, -2.68132617805781232825E-3, 3.47222221605458667310E-3, 8.33333333333482257126E-2, }; #define MAXSTIR 143.01608 static double SQTPI = 2.50662827463100050242E0; #endif #if DEC static unsigned short STIR[20] = { 0035516, 0061622, 0144553, 0112224, 0135160, 0131531, 0037460, 0165740, 0136057, 0134460, 0037242, 0077270, 0036143, 0107070, 0156306, 0027751, 0037252, 0125252, 0125252, 0146064, }; #define MAXSTIR 26.77 static unsigned short SQT[4] = { 0040440, 0066230, 0177661, 0034055, }; #define SQTPI *(double *)SQT #endif #if IBMPC static unsigned short STIR[20] = { 0x7293, 0x592d, 0xcc72, 0x3f49, 0x1d7c, 0x27e6, 0x166b, 0xbf2e, 0x4fd7, 0x07d4, 0xf726, 0xbf65, 0xc5fd, 0x1b98, 0x71c7, 0x3f6c, 0x5986, 0x5555, 0x5555, 0x3fb5, }; #define MAXSTIR 143.01608 static unsigned short SQT[4] = { 0x2706, 0x1ff6, 0x0d93, 0x4004, }; #define SQTPI *(double *)SQT #endif #if MIEEE static unsigned short STIR[20] = { 0x3f49, 0xcc72, 0x592d, 0x7293, 0xbf2e, 0x166b, 0x27e6, 0x1d7c, 0xbf65, 0xf726, 0x07d4, 0x4fd7, 0x3f6c, 0x71c7, 0x1b98, 0xc5fd, 0x3fb5, 0x5555, 0x5555, 0x5986, }; #define MAXSTIR 143.01608 static unsigned short SQT[4] = { 0x4004, 0x0d93, 0x1ff6, 0x2706, }; #define SQTPI *(double *)SQT #endif int sgngam = 0; extern int sgngam; extern double MAXLOG, MAXNUM; #ifndef PI extern double PI; #endif #ifdef INFINITIES extern double INFINITY; #endif /* Gamma function computed by Stirling's formula. * The polynomial STIR is valid for 33 <= x <= 172. */ static double stirf (double x) { double y, w, v; w = 1.0 / x; w = 1.0 + w * polevl (w, STIR, 4); y = exp (x); if (x > MAXSTIR) { /* Avoid overflow in pow() */ v = pow (x, 0.5 * x - 0.25); y = v * (v / y); } else { y = pow (x, x - 0.5) / y; } y = SQTPI * y * w; return (y); } double gammafn (double x) { double p, q, z; int i; sgngam = 1; if (g_isnan (x)) return (x); #ifdef INFINITIES if (x == G_INFINITY) return (x); if (x == -G_INFINITY) return (G_NAN); #endif if (!g_finite (x)) return (x); q = fabs (x); if (q > 33.0) { if (x < 0.0) { p = floor (q); if (p == q) { gamnan: mtherr ("gamma", MATHERR_DOMAIN); return (G_NAN); } i = (int) p; if ((i & 1) == 0) sgngam = -1; z = q - p; if (z > 0.5) { p += 1.0; z = q - p; } z = q * sin (PI * z); if (z == 0.0) { #ifdef INFINITIES return (sgngam * G_INFINITY); #else /*goverf:*/ mtherr ("gamma", MATHERR_OVERFLOW); return (sgngam * MAXNUM); #endif } z = fabs (z); z = PI / (z * stirf (q)); } else { z = stirf (x); } return (sgngam * z); } z = 1.0; while (x >= 3.0) { x -= 1.0; z *= x; } while (x < 0.0) { if (x > -1.E-9) goto small; z /= x; x += 1.0; } while (x < 2.0) { if (x < 1.e-9) goto small; z /= x; x += 1.0; } if (x == 2.0) return (z); x -= 2.0; p = polevl (x, P, 6); q = polevl (x, Q, 7); return (z * p / q); small: if (x == 0.0) { goto gamnan; } else return (z / ((1.0 + 0.5772156649015329 * x) * x)); } /* A[]: Stirling's formula expansion of log gamma * B[], C[]: log gamma function between 2 and 3 */ #ifdef UNK static double A[] = { 8.11614167470508450300E-4, -5.95061904284301438324E-4, 7.93650340457716943945E-4, -2.77777777730099687205E-3, 8.33333333333331927722E-2 }; static double B[] = { -1.37825152569120859100E3, -3.88016315134637840924E4, -3.31612992738871184744E5, -1.16237097492762307383E6, -1.72173700820839662146E6, -8.53555664245765465627E5 }; static double C[] = { /* 1.00000000000000000000E0, */ -3.51815701436523470549E2, -1.70642106651881159223E4, -2.20528590553854454839E5, -1.13933444367982507207E6, -2.53252307177582951285E6, -2.01889141433532773231E6 }; /* log( sqrt( 2*pi ) ) */ static double LS2PI = 0.91893853320467274178; #define MAXLGM 2.556348e305 #endif #ifdef DEC static unsigned short A[] = { 0035524, 0141201, 0034633, 0031405, 0135433, 0176755, 0126007, 0045030, 0035520, 0006371, 0003342, 0172730, 0136066, 0005540, 0132605, 0026407, 0037252, 0125252, 0125252, 0125132 }; static unsigned short B[] = { 0142654, 0044014, 0077633, 0035410, 0144027, 0110641, 0125335, 0144760, 0144641, 0165637, 0142204, 0047447, 0145215, 0162027, 0146246, 0155211, 0145322, 0026110, 0010317, 0110130, 0145120, 0061472, 0120300, 0025363 }; static unsigned short C[] = { /*0040200,0000000,0000000,0000000*/ 0142257, 0164150, 0163630, 0112622, 0143605, 0050153, 0156116, 0135272, 0144527, 0056045, 0145642, 0062332, 0145213, 0012063, 0106250, 0001025, 0145432, 0111254, 0044577, 0115142, 0145366, 0071133, 0050217, 0005122 }; /* log( sqrt( 2*pi ) ) */ static unsigned short LS2P[] = { 040153, 037616, 041445, 0172645, }; #define LS2PI *(double *)LS2P #define MAXLGM 2.035093e36 #endif #ifdef IBMPC static unsigned short A[] = { 0x6661, 0x2733, 0x9850, 0x3f4a, 0xe943, 0xb580, 0x7fbd, 0xbf43, 0x5ebb, 0x20dc, 0x019f, 0x3f4a, 0xa5a1, 0x16b0, 0xc16c, 0xbf66, 0x554b, 0x5555, 0x5555, 0x3fb5 }; static unsigned short B[] = { 0x6761, 0x8ff3, 0x8901, 0xc095, 0xb93e, 0x355b, 0xf234, 0xc0e2, 0x89e5, 0xf890, 0x3d73, 0xc114, 0xdb51, 0xf994, 0xbc82, 0xc131, 0xf20b, 0x0219, 0x4589, 0xc13a, 0x055e, 0x5418, 0x0c67, 0xc12a }; static unsigned short C[] = { /*0x0000,0x0000,0x0000,0x3ff0,*/ 0x12b2, 0x1cf3, 0xfd0d, 0xc075, 0xd757, 0x7b89, 0xaa0d, 0xc0d0, 0x4c9b, 0xb974, 0xeb84, 0xc10a, 0x0043, 0x7195, 0x6286, 0xc131, 0xf34c, 0x892f, 0x5255, 0xc143, 0xe14a, 0x6a11, 0xce4b, 0xc13e }; /* log( sqrt( 2*pi ) ) */ static unsigned short LS2P[] = { 0xbeb5, 0xc864, 0x67f1, 0x3fed }; #define LS2PI *(double *)LS2P #define MAXLGM 2.556348e305 #endif #ifdef MIEEE static unsigned short A[] = { 0x3f4a, 0x9850, 0x2733, 0x6661, 0xbf43, 0x7fbd, 0xb580, 0xe943, 0x3f4a, 0x019f, 0x20dc, 0x5ebb, 0xbf66, 0xc16c, 0x16b0, 0xa5a1, 0x3fb5, 0x5555, 0x5555, 0x554b }; static unsigned short B[] = { 0xc095, 0x8901, 0x8ff3, 0x6761, 0xc0e2, 0xf234, 0x355b, 0xb93e, 0xc114, 0x3d73, 0xf890, 0x89e5, 0xc131, 0xbc82, 0xf994, 0xdb51, 0xc13a, 0x4589, 0x0219, 0xf20b, 0xc12a, 0x0c67, 0x5418, 0x055e }; static unsigned short C[] = { 0xc075, 0xfd0d, 0x1cf3, 0x12b2, 0xc0d0, 0xaa0d, 0x7b89, 0xd757, 0xc10a, 0xeb84, 0xb974, 0x4c9b, 0xc131, 0x6286, 0x7195, 0x0043, 0xc143, 0x5255, 0x892f, 0xf34c, 0xc13e, 0xce4b, 0x6a11, 0xe14a }; /* log( sqrt( 2*pi ) ) */ static unsigned short LS2P[] = { 0x3fed, 0x67f1, 0xc864, 0xbeb5 }; #define LS2PI *(double *)LS2P #define MAXLGM 2.556348e305 #endif /* Logarithm of gamma function */ double lgam (double x) { double p, q, u, w, z; int i; sgngam = 1; if (g_isnan (x)) return (x); #ifdef INFINITIES if (!g_finite (x)) return (G_INFINITY); #endif if (x < -34.0) { q = -x; w = lgam (q); /* note this modifies sgngam! */ p = floor (q); if (p == q) { lgsing: #ifdef INFINITIES mtherr ("lgam", MATHERR_SING); return (G_INFINITY); #else goto loverf; #endif } i = (int) p; if ((i & 1) == 0) sgngam = -1; else sgngam = 1; z = q - p; if (z > 0.5) { p += 1.0; z = p - q; } z = q * sin (PI * z); if (z == 0.0) goto lgsing; /* z = log(PI) - log( z ) - w;*/ z = LOGPI - log (z) - w; return (z); } if (x < 13.0) { z = 1.0; p = 0.0; u = x; while (u >= 3.0) { p -= 1.0; u = x + p; z *= u; } while (u < 2.0) { if (u == 0.0) goto lgsing; z /= u; p += 1.0; u = x + p; } if (z < 0.0) { sgngam = -1; z = -z; } else sgngam = 1; if (u == 2.0) return (log (z)); p -= 2.0; x = x + p; p = x * polevl (x, B, 5) / p1evl (x, C, 6); return (log (z) + p); } if (x > MAXLGM) { #ifdef INFINITIES return (sgngam * G_INFINITY); #else loverf: mtherr ("lgam", MATHERR_OVERFLOW); return (sgngam * MAXNUM); #endif } q = (x - 0.5) * log (x) - x + LS2PI; if (x > 1.0e8) return (q); p = 1.0 / (x * x); if (x >= 1000.0) q += ((7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3) * p + 0.0833333333333333333333) / x; else q += polevl (p, A, 4) / x; return (q); }