diff options
Diffstat (limited to 'libm/float/gammaf.c')
| -rw-r--r-- | libm/float/gammaf.c | 423 |
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 ); +} |