uClibc now has a math library. muahahahaha!

 -Erik

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 );
+}