uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 423 insertions, 0 deletions

diff --git a/libm/float/gammaf.c b/libm/float/gammaf.c
new file mode 100644
index 000000000..e8c4694c4
--- /dev/null
+++ b/libm/float/gammaf.c
@@ -0,0 +1,423 @@
+/*							gammaf.c
+ *
+ *	Gamma function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, gammaf();
+ * extern int sgngamf;
+ *
+ * y = gammaf( 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 sgngamf.
+ * This same variable is also filled in by the logarithmic
+ * gamma function lgam().
+ *
+ * Arguments between 0 and 10 are reduced by recurrence and the
+ * function is approximated by a polynomial function covering
+ * the interval (2,3).  Large arguments are handled by Stirling's
+ * formula. Negative arguments are made positive using
+ * a reflection formula.  
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE       0,-33      100,000     5.7e-7      1.0e-7
+ *    IEEE       -33,0      100,000     6.1e-7      1.2e-7
+ *
+ *
+ */
+/*							lgamf()
+ *
+ *	Natural logarithm of gamma function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, lgamf();
+ * extern int sgngamf;
+ *
+ * y = lgamf( 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 sgngamf.
+ *
+ * For arguments greater than 6.5, the logarithm of the gamma
+ * function is approximated by the logarithmic version of
+ * Stirling's formula.  Arguments between 0 and +6.5 are reduced by
+ * by recurrence to the interval [.75,1.25] or [1.5,2.5] of a rational
+ * approximation.  The cosecant reflection formula is employed for
+ * arguments less than zero.
+ *
+ * Arguments greater than MAXLGM = 2.035093e36 return MAXNUM and an
+ * error message.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *
+ *
+ * arithmetic      domain        # trials     peak         rms
+ *    IEEE        -100,+100       500,000    7.4e-7       6.8e-8
+ * The error criterion was relative when the function magnitude
+ * was greater than one but absolute when it was less than one.
+ * The routine has low relative error for positive arguments.
+ *
+ * The following test used the relative error criterion.
+ *    IEEE    -2, +3              100000     4.0e-7      5.6e-8
+ *
+ */
+
+/*							gamma.c	*/
+/*	gamma function	*/
+
+/*
+Cephes Math Library Release 2.7:  July, 1998
+Copyright 1984, 1987, 1989, 1992, 1998 by Stephen L. Moshier
+*/
+
+
+#include <math.h>
+
+/* define MAXGAM 34.84425627277176174 */
+
+/* Stirling's formula for the gamma function
+ * gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) ( 1 + 1/x P(1/x) )
+ * .028 < 1/x < .1
+ * relative error < 1.9e-11
+ */
+static float STIR[] = {
+-2.705194986674176E-003,
+ 3.473255786154910E-003,
+ 8.333331788340907E-002,
+};
+static float MAXSTIR = 26.77;
+static float SQTPIF = 2.50662827463100050242; /* sqrt( 2 pi ) */
+
+int sgngamf = 0;
+extern int sgngamf;
+extern float MAXLOGF, MAXNUMF, PIF;
+
+#ifdef ANSIC
+float expf(float);
+float logf(float);
+float powf( float, float );
+float sinf(float);
+float gammaf(float);
+float floorf(float);
+static float stirf(float);
+float polevlf( float, float *, int );
+float p1evlf( float, float *, int );
+#else
+float expf(), logf(), powf(), sinf(), floorf();
+float polevlf(), p1evlf();
+static float stirf();
+#endif
+
+/* Gamma function computed by Stirling's formula,
+ * sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
+ * The polynomial STIR is valid for 33 <= x <= 172.
+ */
+static float stirf( float xx )
+{
+float x, y, w, v;
+
+x = xx;
+w = 1.0/x;
+w = 1.0 + w * polevlf( w, STIR, 2 );
+y = expf( -x );
+if( x > MAXSTIR )
+	{ /* Avoid overflow in pow() */
+	v = powf( x, 0.5 * x - 0.25 );
+	y *= v;
+	y *= v;
+	}
+else
+	{
+	y = powf( x, x - 0.5 ) * y;
+	}
+y = SQTPIF * y * w;
+return( y );
+}
+
+
+/* gamma(x+2), 0 < x < 1 */
+static float P[] = {
+ 1.536830450601906E-003,
+ 5.397581592950993E-003,
+ 4.130370201859976E-003,
+ 7.232307985516519E-002,
+ 8.203960091619193E-002,
+ 4.117857447645796E-001,
+ 4.227867745131584E-001,
+ 9.999999822945073E-001,
+};
+
+float gammaf( float xx )
+{
+float p, q, x, z, nz;
+int i, direction, negative;
+
+x = xx;
+sgngamf = 1;
+negative = 0;
+nz = 0.0;
+if( x < 0.0 )
+	{
+	negative = 1;
+	q = -x;
+	p = floorf(q);
+	if( p == q )
+		goto goverf;
+	i = p;
+	if( (i & 1) == 0 )
+		sgngamf = -1;
+	nz = q - p;
+	if( nz > 0.5 )
+		{
+		p += 1.0;
+		nz = q - p;
+		}
+	nz = q * sinf( PIF * nz );
+	if( nz == 0.0 )
+		{
+goverf:
+		mtherr( "gamma", OVERFLOW );
+		return( sgngamf * MAXNUMF);
+		}
+	if( nz < 0 )
+		nz = -nz;
+	x = q;
+	}
+if( x >= 10.0 )
+	{
+	z = stirf(x);
+	}
+if( x < 2.0 )
+	direction = 1;
+else
+	direction = 0;
+z = 1.0;
+while( x >= 3.0 )
+	{
+	x -= 1.0;
+	z *= x;
+	}
+/*
+while( x < 0.0 )
+	{
+	if( x > -1.E-4 )
+		goto small;
+	z *=x;
+	x += 1.0;
+	}
+*/
+while( x < 2.0 )
+	{
+	if( x < 1.e-4 )
+		goto small;
+	z *=x;
+	x += 1.0;
+	}
+
+if( direction )
+	z = 1.0/z;
+
+if( x == 2.0 )
+	return(z);
+
+x -= 2.0;
+p = z * polevlf( x, P, 7 );
+
+gdone:
+
+if( negative )
+	{
+	p = sgngamf * PIF/(nz * p );
+	}
+return(p);
+
+small:
+if( x == 0.0 )
+	{
+	mtherr( "gamma", SING );
+	return( MAXNUMF );
+	}
+else
+	{
+	p = z / ((1.0 + 0.5772156649015329 * x) * x);
+	goto gdone;
+	}
+}
+
+
+
+
+/* log gamma(x+2), -.5 < x < .5 */
+static float B[] = {
+ 6.055172732649237E-004,
+-1.311620815545743E-003,
+ 2.863437556468661E-003,
+-7.366775108654962E-003,
+ 2.058355474821512E-002,
+-6.735323259371034E-002,
+ 3.224669577325661E-001,
+ 4.227843421859038E-001
+};
+
+/* log gamma(x+1), -.25 < x < .25 */
+static float C[] = {
+ 1.369488127325832E-001,
+-1.590086327657347E-001,
+ 1.692415923504637E-001,
+-2.067882815621965E-001,
+ 2.705806208275915E-001,
+-4.006931650563372E-001,
+ 8.224670749082976E-001,
+-5.772156501719101E-001
+};
+
+/* log( sqrt( 2*pi ) ) */
+static float LS2PI  =  0.91893853320467274178;
+#define MAXLGM 2.035093e36
+static float PIINV =  0.318309886183790671538;
+
+/* Logarithm of gamma function */
+
+
+float lgamf( float xx )
+{
+float p, q, w, z, x;
+float nx, tx;
+int i, direction;
+
+sgngamf = 1;
+
+x = xx;
+if( x < 0.0 )
+	{
+	q = -x;
+	w = lgamf(q); /* note this modifies sgngam! */
+	p = floorf(q);
+	if( p == q )
+		goto loverf;
+	i = p;
+	if( (i & 1) == 0 )
+		sgngamf = -1;
+	else
+		sgngamf = 1;
+	z = q - p;
+	if( z > 0.5 )
+		{
+		p += 1.0;
+		z = p - q;
+		}
+	z = q * sinf( PIF * z );
+	if( z == 0.0 )
+		goto loverf;
+	z = -logf( PIINV*z ) - w;
+	return( z );
+	}
+
+if( x < 6.5 )
+	{
+	direction = 0;
+	z = 1.0;
+	tx = x;
+	nx = 0.0;
+	if( x >= 1.5 )
+		{
+		while( tx > 2.5 )
+			{
+			nx -= 1.0;
+			tx = x + nx;
+			z *=tx;
+			}
+		x += nx - 2.0;
+iv1r5:
+		p = x * polevlf( x, B, 7 );
+		goto cont;
+		}
+	if( x >= 1.25 )
+		{
+		z *= x;
+		x -= 1.0; /* x + 1 - 2 */
+		direction = 1;
+		goto iv1r5;
+		}
+	if( x >= 0.75 )
+		{
+		x -= 1.0;
+		p = x * polevlf( x, C, 7 );
+		q = 0.0;
+		goto contz;
+		}
+	while( tx < 1.5 )
+		{
+		if( tx == 0.0 )
+			goto loverf;
+		z *=tx;
+		nx += 1.0;
+		tx = x + nx;
+		}
+	direction = 1;
+	x += nx - 2.0;
+	p = x * polevlf( x, B, 7 );
+
+cont:
+	if( z < 0.0 )
+		{
+		sgngamf = -1;
+		z = -z;
+		}
+	else
+		{
+		sgngamf = 1;
+		}
+	q = logf(z);
+	if( direction )
+		q = -q;
+contz:
+	return( p + q );
+	}
+
+if( x > MAXLGM )
+	{
+loverf:
+	mtherr( "lgamf", OVERFLOW );
+	return( sgngamf * MAXNUMF );
+	}
+
+/* Note, though an asymptotic formula could be used for x >= 3,
+ * there is cancellation error in the following if x < 6.5.  */
+q = LS2PI - x;
+q += ( x - 0.5 ) * logf(x);
+
+if( x <= 1.0e4 )
+	{
+	z = 1.0/x;
+	p = z * z;
+	q += ((    6.789774945028216E-004 * p
+		 - 2.769887652139868E-003 ) * p
+		+  8.333316229807355E-002 ) * z;
+	}
+return( q );
+}