Totally rework the math library, this time based on the MacOs X

math library (which is itself based on the math lib from FreeBSD).
 -Erik

1 files changed, 0 insertions, 423 deletions

diff --git a/libm/float/gammaf.c b/libm/float/gammaf.c
deleted file mode 100644
index e8c4694c4..000000000
--- a/libm/float/gammaf.c
+++ /dev/null
@@ -1,423 +0,0 @@
-/*							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 );
-}