/* gamma.c * * Gamma function * * * * SYNOPSIS: * * double x, y, gamma(); * extern int sgngam; * * y = gamma( 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.8: June, 2000 Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier */ #include <math.h> #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, PI; #ifdef ANSIPROT extern double pow ( double, double ); extern double log ( double ); extern double exp ( double ); extern double sin ( double ); extern double polevl ( double, void *, int ); extern double p1evl ( double, void *, int ); extern double floor ( double ); extern double fabs ( double ); extern int isnan ( double ); extern int isfinite ( double ); static double stirf ( double ); double lgam ( double ); #else double pow(), log(), exp(), sin(), polevl(), p1evl(), floor(), fabs(); int isnan(), isfinite(); static double stirf(); double lgam(); #endif #ifdef INFINITIES extern double INFINITY; #endif #ifdef NANS extern double NAN; #endif /* Gamma function computed by Stirling's formula. * The polynomial STIR is valid for 33 <= x <= 172. */ static double stirf(x) 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 gamma(x) double x; { double p, q, z; int i; sgngam = 1; #ifdef NANS if( isnan(x) ) return(x); #endif #ifdef INFINITIES #ifdef NANS if( x == INFINITY ) return(x); if( x == -INFINITY ) return(NAN); #else if( !isfinite(x) ) return(x); #endif #endif q = fabs(x); if( q > 33.0 ) { if( x < 0.0 ) { p = floor(q); if( p == q ) { #ifdef NANS gamnan: mtherr( "gamma", DOMAIN ); return (NAN); #else goto goverf; #endif } i = 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 * INFINITY); #else goverf: mtherr( "gamma", 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 ) { #ifdef INFINITIES #ifdef NANS goto gamnan; #else return( INFINITY ); #endif #else mtherr( "gamma", SING ); return( MAXNUM ); #endif } 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(x) double x; { double p, q, u, w, z; int i; sgngam = 1; #ifdef NANS if( isnan(x) ) return(x); #endif #ifdef INFINITIES if( !isfinite(x) ) return(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", SING ); return (INFINITY); #else goto loverf; #endif } i = 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 * INFINITY ); #else loverf: mtherr( "lgam", 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 ); }