summaryrefslogtreecommitdiff
path: root/libm/double/gamma.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/gamma.c')
-rw-r--r--libm/double/gamma.c685
1 files changed, 685 insertions, 0 deletions
diff --git a/libm/double/gamma.c b/libm/double/gamma.c
new file mode 100644
index 000000000..341b4e915
--- /dev/null
+++ b/libm/double/gamma.c
@@ -0,0 +1,685 @@
+/* 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 );
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-11 14:08:53 +0000