diff options
Diffstat (limited to 'winsup/mingw/mingwex/math/tgammaf.c')
| -rw-r--r-- | winsup/mingw/mingwex/math/tgammaf.c | 265 |
1 files changed, 0 insertions, 265 deletions
diff --git a/winsup/mingw/mingwex/math/tgammaf.c b/winsup/mingw/mingwex/math/tgammaf.c deleted file mode 100644 index 07d294971..000000000 --- a/winsup/mingw/mingwex/math/tgammaf.c +++ /dev/null @@ -1,265 +0,0 @@ -/* gammaf.c - * - * Gamma function - * - * - * - * SYNOPSIS: - * - * float x, y, __tgammaf_r(); - * int* sgngamf; - * y = __tgammaf_r( x, sgngamf ); - * - * float x, y, tgammaf(); - * y = tgammaf( x); - * - * - * DESCRIPTION: - * - * Returns gamma function of the argument. The result is - * correctly signed. In the reentrant version the sign (+1 or -1) - * is returned in the variable referenced by sgngamf. - * - * 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 - * - * - */ - -/* -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> - */ - - -#ifndef __MINGW32__ -#include "mconf.h" -#else -#include "cephes_mconf.h" -#endif - -/* 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 const float STIR[] = { --2.705194986674176E-003, - 3.473255786154910E-003, - 8.333331788340907E-002, -}; -static const float MAXSTIR = 26.77; -static const float SQTPIF = 2.50662827463100050242; /* sqrt( 2 pi ) */ - -#ifndef __MINGW32__ - -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 - -#else /* __MINGW32__ */ -static float stirf(float); -#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 x ) -{ -float y, w, v; - -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 const 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 __tgammaf_r( float x, int* sgngamf) -{ -float p, q, z, nz; -int i, direction, negative; - -#ifdef NANS -if( isnan(x) ) - return(x); -#endif -#ifdef INFINITIES -#ifdef NANS -if( x == INFINITYF ) - return(x); -if( x == -INFINITYF ) - return(NANF); -#else -if( !isfinite(x) ) - return(x); -#endif -#endif - -*sgngamf = 1; -negative = 0; -nz = 0.0; -if( x < 0.0 ) - { - negative = 1; - q = -x; - p = floorf(q); - if( p == q ) - { -gsing: - _SET_ERRNO(EDOM); - mtherr( "tgammaf", SING ); -#ifdef INFINITIES - return (INFINITYF); -#else - return (MAXNUMF); -#endif - } - 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 ) - { - _SET_ERRNO(ERANGE); - mtherr( "tgamma", OVERFLOW ); -#ifdef INFINITIES - return( *sgngamf * INFINITYF); -#else - return( *sgngamf * MAXNUMF); -#endif - } - 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 ) - { - goto gsing; - } -else - { - p = z / ((1.0 + 0.5772156649015329 * x) * x); - goto gdone; - } -} - -/* This is the C99 version */ - -float tgammaf(float x) -{ - int local_sgngamf=0; - return (__tgammaf_r(x, &local_sgngamf)); -} |