diff options
Diffstat (limited to 'winsup/mingw/mingwex/math/lgammaf.c')
| -rw-r--r-- | winsup/mingw/mingwex/math/lgammaf.c | 253 |
1 files changed, 0 insertions, 253 deletions
diff --git a/winsup/mingw/mingwex/math/lgammaf.c b/winsup/mingw/mingwex/math/lgammaf.c deleted file mode 100644 index 20982f999..000000000 --- a/winsup/mingw/mingwex/math/lgammaf.c +++ /dev/null @@ -1,253 +0,0 @@ -/* 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)); -} |