diff options
Diffstat (limited to 'winsup/mingw/mingwex/math/lgammaf.c')
-rw-r--r-- | winsup/mingw/mingwex/math/lgammaf.c | 253 |
1 files changed, 253 insertions, 0 deletions
diff --git a/winsup/mingw/mingwex/math/lgammaf.c b/winsup/mingw/mingwex/math/lgammaf.c new file mode 100644 index 00000000000..20982f999fc --- /dev/null +++ b/winsup/mingw/mingwex/math/lgammaf.c @@ -0,0 +1,253 @@ +/* lgamf() + * + * Natural logarithm of gamma function + * + * + * + * SYNOPSIS: + * + * float x, y, __lgammaf_r(); + * int* sgngamf; + * y = __lgammaf_r( x, sgngamf ); + * + * float x, y, lgammaf(); + * y = lgammaf( x); + * + * + * + * DESCRIPTION: + * + * Returns the base e (2.718...) logarithm of the absolute + * value of the gamma function of the argument. In the reentrant + * version the sign (+1 or -1) of the gamma function is returned in + * variable referenced by 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 + * + */ + + +/* + Cephes Math Library Release 2.7: July, 1998 + Copyright 1984, 1987, 1989, 1992, 1998 by Stephen L. Moshier +*/ + +/* + 26-11-2002 Modified for mingw. + Danny Smith <dannysmith@users.sourceforge.net> +*/ + + +/* log gamma(x+2), -.5 < x < .5 */ +static const 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 const 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 const float LS2PI = 0.91893853320467274178; +#define MAXLGM 2.035093e36 +static const float PIINV = 0.318309886183790671538; + +#ifndef __MINGW32__ +#include "mconf.h" +float floorf(float); +float polevlf( float, float *, int ); +float p1evlf( float, float *, int ); +#else +#include "cephes_mconf.h" +#endif + +/* Reentrant version */ +/* Logarithm of gamma function */ + +float __lgammaf_r( float x, int* sgngamf ) +{ +float p, q, w, z; +float nx, tx; +int i, direction; + +*sgngamf = 1; +#ifdef NANS +if( isnan(x) ) + return(x); +#endif + +#ifdef INFINITIES +if( !isfinite(x) ) + return(x); +#endif + + +if( x < 0.0 ) + { + q = -x; + w = __lgammaf_r(q, sgngamf); /* note this modifies sgngam! */ + p = floorf(q); + if( p == q ) + { +lgsing: + _SET_ERRNO(EDOM); + mtherr( "lgamf", SING ); +#ifdef INFINITIES + return (INFINITYF); +#else + return( *sgngamf * MAXNUMF ); +#endif + } + 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 lgsing; + 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 lgsing; + 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 ) + { + _SET_ERRNO(ERANGE); + mtherr( "lgamf", OVERFLOW ); +#ifdef INFINITIES + return( *sgngamf * INFINITYF ); +#else + return( *sgngamf * MAXNUMF ); +#endif + + } + +/* 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 ); +} + +/* This is the C99 version */ + +float lgammaf(float x) +{ + int local_sgngamf=0; + return (__lgammaf_r(x, &local_sgngamf)); +} |