/* 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 <gnan.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;
#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);
}
syntax highlighted by Code2HTML, v. 0.9.1