diff options
author | Danny Smith <dannysmith@users.sourceforge.net> | 2002-11-27 06:41:25 +0300 |
---|---|---|
committer | Danny Smith <dannysmith@users.sourceforge.net> | 2002-11-27 06:41:25 +0300 |
commit | dc8597f966ca83f46897e7e269ede83566b98177 (patch) | |
tree | a6bf0cb29149d2cd72786d39eb444a73939a70b4 /winsup | |
parent | eb6d2e2f9a8be898e936b6ac07135d57a566014b (diff) |
* mingwex/math/lgamma.c: New file.
* mingwex/math/lgammaf.c: New file.
* mingwex/math/lgammal.c: New file.
* mingwex/math/tgamma.c: New file.
* mingwex/math/tgammaf.c: New file.
* mingwex/math/tgammal.c: New file.
* mingwex/math/cephes_mconf (polevlf): Add float version.
(p1evlf): Likewise.
Define _CEPHES_USE_ERRNO.
* mingwex/Makefile.in (MATH_DISTFILES): Add new files.
(MATH_OBJS): Add new objects.
* include/math.h (lgamma[fl]): Add prototypes.
(tgamma[fl]): Add prototypes.
Diffstat (limited to 'winsup')
-rw-r--r-- | winsup/mingw/include/math.h | 16 | ||||
-rw-r--r-- | winsup/mingw/mingwex/Makefile.in | 12 | ||||
-rw-r--r-- | winsup/mingw/mingwex/math/cephes_mconf.h | 98 | ||||
-rw-r--r-- | winsup/mingw/mingwex/math/lgamma.c | 359 | ||||
-rw-r--r-- | winsup/mingw/mingwex/math/lgammaf.c | 253 | ||||
-rw-r--r-- | winsup/mingw/mingwex/math/lgammal.c | 416 | ||||
-rw-r--r-- | winsup/mingw/mingwex/math/tgamma.c | 385 | ||||
-rw-r--r-- | winsup/mingw/mingwex/math/tgammaf.c | 265 | ||||
-rw-r--r-- | winsup/mingw/mingwex/math/tgammal.c | 501 |
9 files changed, 2300 insertions, 5 deletions
diff --git a/winsup/mingw/include/math.h b/winsup/mingw/include/math.h index fbfd712c0..93dcb5f09 100644 --- a/winsup/mingw/include/math.h +++ b/winsup/mingw/include/math.h @@ -501,7 +501,21 @@ extern long double powl (long double, long double); extern float sqrtf (float); extern long double sqrtl (long double); -/* 7.12.8 Error and gamma functions: TODO */ +/* TODO */ +/* 7.12.8.1 The erf functions */ +/* 7.12.8.2 The erfc functions */ + +/* 7.12.8.3 The lgamma functions */ + +extern double lgamma (double); +extern float lgammaf (float); +extern long double lgammal (long double); + +/* 77.12.8.4 The tgamma functions */ + +extern double tgamma (double); +extern float tgammaf (float); +extern long double tgammal (long double); /* 7.12.9.1 Double in C89 */ extern float ceilf (float); diff --git a/winsup/mingw/mingwex/Makefile.in b/winsup/mingw/mingwex/Makefile.in index 318613728..c38ada0d2 100644 --- a/winsup/mingw/mingwex/Makefile.in +++ b/winsup/mingw/mingwex/Makefile.in @@ -46,7 +46,8 @@ MATH_DISTFILES = \ fmodl.c fp_consts.c fp_consts.h fp_constsf.c fp_constsl.c \ fpclassify.c fpclassifyf.c fpclassifyl.c \ frexpf.c frexpl.S fucom.c hypotf.c hypotl.c ilogb.S ilogbf.S \ - ilogbl.S isnan.c isnanf.c isnanl.c ldexpf.c ldexpl.c llrint.c \ + ilogbl.S isnan.c isnanf.c isnanl.c ldexpf.c ldexpl.c \ + lgamma.c lgammaf.c lgammal.c llrint.c \ llrintf.c llrintl.c llround.c llroundf.c llroundl.c \ log10f.S log10l.S log1p.S log1pf.S log1pl.S log2.S log2f.S \ log2l.S logb.c logbf.c logbl.c logf.S logl.S lrint.c lrintf.c \ @@ -57,7 +58,8 @@ MATH_DISTFILES = \ remquof.S remquol.S rint.c rintf.c rintl.c round.c roundf.c \ roundl.c scalbn.S scalbnf.S scalbnl.S signbit.c signbitf.c \ signbitl.c sinf.S sinhf.c sinhl.c sinl.S sqrtf.c sqrtl.c \ - tanf.S tanhf.c tanhl.c tanl.S trunc.c truncf.c truncl.c + tanf.S tanhf.c tanhl.c tanl.S tgamma.c tgammaf.c tgammal.c \ + trunc.c truncf.c truncl.c CC = @CC@ # FIXME: Which is it, CC or CC_FOR_TARGET? @@ -115,7 +117,8 @@ MATH_OBJS = \ fmodl.o fp_consts.o fp_constsf.o fp_constsl.o \ fpclassify.o fpclassifyf.o fpclassifyl.o \ frexpf.o frexpl.o fucom.o hypotf.o hypotl.o ilogb.o ilogbf.o \ - ilogbl.o isnan.o isnanf.o isnanl.o ldexpf.o ldexpl.o llrint.o \ + ilogbl.o isnan.o isnanf.o isnanl.o ldexpf.o ldexpl.o \ + lgamma.o lgammaf.o lgammal.o llrint.o \ llrintf.o llrintl.o llround.o llroundf.o llroundl.o \ log10f.o log10l.o log1p.o log1pf.o log1pl.o log2.o log2f.o \ log2l.o logb.o logbf.o logbl.o logf.o logl.o lrint.o lrintf.o \ @@ -126,7 +129,8 @@ MATH_OBJS = \ remquof.o remquol.o rint.o rintf.o rintl.o round.o roundf.o \ roundl.o scalbn.o scalbnf.o scalbnl.o signbit.o signbitf.o \ signbitl.o sinf.o sinhf.o sinhl.o sinl.o sqrtf.o sqrtl.o \ - tanf.o tanhf.o tanhl.o tanl.o trunc.o truncf.o truncl.o + tanf.o tanhf.o tanhl.o tanl.o tgamma.o tgammaf.o tgammal.o \ + trunc.o truncf.o truncl.o FENV_OBJS = fesetround.o fegetround.o \ fegetenv.o fesetenv.o feupdateenv.o \ feclearexcept.o feholdexcept.o fegetexceptflag.o \ diff --git a/winsup/mingw/mingwex/math/cephes_mconf.h b/winsup/mingw/mingwex/math/cephes_mconf.h index 1dda63d53..85e0bdcf0 100644 --- a/winsup/mingw/mingwex/math/cephes_mconf.h +++ b/winsup/mingw/mingwex/math/cephes_mconf.h @@ -12,6 +12,8 @@ #define mtherr(fname, code) #define XPD 0, +#define _CEPHES_USE_ERRNO + #ifdef _CEPHES_USE_ERRNO #define _SET_ERRNO(x) errno = (x) #else @@ -275,3 +277,99 @@ while( --n ); return( y ); } +/* Float version */ + +/* polevlf.c + * p1evlf.c + * + * Evaluate polynomial + * + * + * + * SYNOPSIS: + * + * int N; + * float x, y, coef[N+1], polevlf[]; + * + * y = polevlf( x, coef, N ); + * + * + * + * DESCRIPTION: + * + * Evaluates polynomial of degree N: + * + * 2 N + * y = C + C x + C x +...+ C x + * 0 1 2 N + * + * Coefficients are stored in reverse order: + * + * coef[0] = C , ..., coef[N] = C . + * N 0 + * + * The function p1evl() assumes that coef[N] = 1.0 and is + * omitted from the array. Its calling arguments are + * otherwise the same as polevl(). + * + * + * SPEED: + * + * In the interest of speed, there are no checks for out + * of bounds arithmetic. This routine is used by most of + * the functions in the library. Depending on available + * equipment features, the user may wish to rewrite the + * program in microcode or assembly language. + * + */ + +/* +Cephes Math Library Release 2.1: December, 1988 +Copyright 1984, 1987, 1988 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +static __inline__ float polevlf(float x, const float* coef, int N ) +{ +float ans; +float *p; +int i; + +p = (float*)coef; +ans = *p++; + +/* +for( i=0; i<N; i++ ) + ans = ans * x + *p++; +*/ + +i = N; +do + ans = ans * x + *p++; +while( --i ); + +return( ans ); +} + +/* p1evl() */ +/* N + * Evaluate polynomial when coefficient of x is 1.0. + * Otherwise same as polevl. + */ + +static __inline__ float p1evlf( float x, const float *coef, int N ) +{ +float ans; +float *p; +int i; + +p = (float*)coef; +ans = x + *p++; +i = N-1; + +do + ans = ans * x + *p++; +while( --i ); + +return( ans ); +} diff --git a/winsup/mingw/mingwex/math/lgamma.c b/winsup/mingw/mingwex/math/lgamma.c new file mode 100644 index 000000000..f85094957 --- /dev/null +++ b/winsup/mingw/mingwex/math/lgamma.c @@ -0,0 +1,359 @@ +/* lgam() + * + * Natural logarithm of gamma function + * + * + * + * SYNOPSIS: + * + * double x, y, __lgamma_r(); + * int* sgngam; + * y = __lgamma_r( x, sgngam ); + * + * double x, y, lgamma(); + * y = lgamma( 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 the variable referenced by sgngam. + * + * For arguments greater than 13, the logarithm of the gamma + * function is approximated by the logarithmic version of + * Stirling's formula using a polynomial approximation of + * degree 4. Arguments between -33 and +33 are reduced by + * recurrence to the interval [2,3] of a rational approximation. + * The cosecant reflection formula is employed for arguments + * less than -33. + * + * Arguments greater than MAXLGM return MAXNUM and an error + * message. MAXLGM = 2.035093e36 for DEC + * arithmetic or 2.556348e305 for IEEE arithmetic. + * + * + * + * ACCURACY: + * + * + * arithmetic domain # trials peak rms + * DEC 0, 3 7000 5.2e-17 1.3e-17 + * DEC 2.718, 2.035e36 5000 3.9e-17 9.9e-18 + * IEEE 0, 3 28000 5.4e-16 1.1e-16 + * IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17 + * The error criterion was relative when the function magnitude + * was greater than one but absolute when it was less than one. + * + * The following test used the relative error criterion, though + * at certain points the relative error could be much higher than + * indicated. + * IEEE -200, -4 10000 4.8e-16 1.3e-16 + * + */ + +/* + * Cephes Math Library Release 2.8: June, 2000 + * Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier + */ + +/* + * 26-11-2002 Modified for mingw. + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + +#ifndef __MINGW32__ +#include "mconf.h" +#ifdef ANSIPROT +extern double pow ( double, double ); +extern double log ( double ); +extern double exp ( double ); +extern double sin ( double ); +extern double polevl ( double, void *, int ); +extern double p1evl ( double, void *, int ); +extern double floor ( double ); +extern double fabs ( double ); +extern int isnan ( double ); +extern int isfinite ( double ); +#else +double pow(), log(), exp(), sin(), polevl(), p1evl(), floor(), fabs(); +int isnan(), isfinite(); +#endif +#ifdef INFINITIES +extern double INFINITY; +#endif +#ifdef NANS +extern double NAN; +#endif +#else /* __MINGW32__ */ +#include "cephes_mconf.h" +#endif /* __MINGW32__ */ + + +/* A[]: Stirling's formula expansion of log gamma + * B[], C[]: log gamma function between 2 and 3 + */ +#ifdef UNK +static double A[] = { + 8.11614167470508450300E-4, +-5.95061904284301438324E-4, + 7.93650340457716943945E-4, +-2.77777777730099687205E-3, + 8.33333333333331927722E-2 +}; +static double B[] = { +-1.37825152569120859100E3, +-3.88016315134637840924E4, +-3.31612992738871184744E5, +-1.16237097492762307383E6, +-1.72173700820839662146E6, +-8.53555664245765465627E5 +}; +static double C[] = { +/* 1.00000000000000000000E0, */ +-3.51815701436523470549E2, +-1.70642106651881159223E4, +-2.20528590553854454839E5, +-1.13933444367982507207E6, +-2.53252307177582951285E6, +-2.01889141433532773231E6 +}; +/* log( sqrt( 2*pi ) ) */ +static double LS2PI = 0.91893853320467274178; +#define MAXLGM 2.556348e305 +static double LOGPI = 1.14472988584940017414; +#endif + +#ifdef DEC +static const unsigned short A[] = { +0035524,0141201,0034633,0031405, +0135433,0176755,0126007,0045030, +0035520,0006371,0003342,0172730, +0136066,0005540,0132605,0026407, +0037252,0125252,0125252,0125132 +}; +static const unsigned short B[] = { +0142654,0044014,0077633,0035410, +0144027,0110641,0125335,0144760, +0144641,0165637,0142204,0047447, +0145215,0162027,0146246,0155211, +0145322,0026110,0010317,0110130, +0145120,0061472,0120300,0025363 +}; +static const unsigned short C[] = { +/*0040200,0000000,0000000,0000000*/ +0142257,0164150,0163630,0112622, +0143605,0050153,0156116,0135272, +0144527,0056045,0145642,0062332, +0145213,0012063,0106250,0001025, +0145432,0111254,0044577,0115142, +0145366,0071133,0050217,0005122 +}; +/* log( sqrt( 2*pi ) ) */ +static const unsigned short LS2P[] = {040153,037616,041445,0172645,}; +#define LS2PI *(double *)LS2P +#define MAXLGM 2.035093e36 +static const unsigned short LPI[4] = { +0040222,0103202,0043475,0006750, +}; +#define LOGPI *(double *)LPI + +#endif + +#ifdef IBMPC +static const unsigned short A[] = { +0x6661,0x2733,0x9850,0x3f4a, +0xe943,0xb580,0x7fbd,0xbf43, +0x5ebb,0x20dc,0x019f,0x3f4a, +0xa5a1,0x16b0,0xc16c,0xbf66, +0x554b,0x5555,0x5555,0x3fb5 +}; +static const unsigned short B[] = { +0x6761,0x8ff3,0x8901,0xc095, +0xb93e,0x355b,0xf234,0xc0e2, +0x89e5,0xf890,0x3d73,0xc114, +0xdb51,0xf994,0xbc82,0xc131, +0xf20b,0x0219,0x4589,0xc13a, +0x055e,0x5418,0x0c67,0xc12a +}; +static const unsigned short C[] = { +/*0x0000,0x0000,0x0000,0x3ff0,*/ +0x12b2,0x1cf3,0xfd0d,0xc075, +0xd757,0x7b89,0xaa0d,0xc0d0, +0x4c9b,0xb974,0xeb84,0xc10a, +0x0043,0x7195,0x6286,0xc131, +0xf34c,0x892f,0x5255,0xc143, +0xe14a,0x6a11,0xce4b,0xc13e +}; +/* log( sqrt( 2*pi ) ) */ +static const unsigned short LS2P[] = { +0xbeb5,0xc864,0x67f1,0x3fed +}; +#define LS2PI *(double *)LS2P +#define MAXLGM 2.556348e305 +static const unsigned short LPI[4] = { +0xa1bd,0x48e7,0x50d0,0x3ff2, +}; +#define LOGPI *(double *)LPI +#endif + +#ifdef MIEEE +static const unsigned short A[] = { +0x3f4a,0x9850,0x2733,0x6661, +0xbf43,0x7fbd,0xb580,0xe943, +0x3f4a,0x019f,0x20dc,0x5ebb, +0xbf66,0xc16c,0x16b0,0xa5a1, +0x3fb5,0x5555,0x5555,0x554b +}; +static const unsigned short B[] = { +0xc095,0x8901,0x8ff3,0x6761, +0xc0e2,0xf234,0x355b,0xb93e, +0xc114,0x3d73,0xf890,0x89e5, +0xc131,0xbc82,0xf994,0xdb51, +0xc13a,0x4589,0x0219,0xf20b, +0xc12a,0x0c67,0x5418,0x055e +}; +static const unsigned short C[] = { +0xc075,0xfd0d,0x1cf3,0x12b2, +0xc0d0,0xaa0d,0x7b89,0xd757, +0xc10a,0xeb84,0xb974,0x4c9b, +0xc131,0x6286,0x7195,0x0043, +0xc143,0x5255,0x892f,0xf34c, +0xc13e,0xce4b,0x6a11,0xe14a +}; +/* log( sqrt( 2*pi ) ) */ +static const unsigned short LS2P[] = { +0x3fed,0x67f1,0xc864,0xbeb5 +}; +#define LS2PI *(double *)LS2P +#define MAXLGM 2.556348e305 +static unsigned short LPI[4] = { +0x3ff2,0x50d0,0x48e7,0xa1bd, +}; +#define LOGPI *(double *)LPI +#endif + + +/* Logarithm of gamma function */ +/* Reentrant version */ + +double __lgamma_r(double x, int* sgngam) +{ +double p, q, u, w, z; +int i; + +*sgngam = 1; +#ifdef NANS +if( isnan(x) ) + return(x); +#endif + +#ifdef INFINITIES +if( !isfinite(x) ) + return(INFINITY); +#endif + +if( x < -34.0 ) + { + q = -x; + w = __lgamma_r(q, sgngam); /* note this modifies sgngam! */ + p = floor(q); + if( p == q ) + { +lgsing: + _SET_ERRNO(EDOM); + mtherr( "lgam", SING ); +#ifdef INFINITIES + return (INFINITY); +#else + return (MAXNUM); +#endif + } + i = p; + if( (i & 1) == 0 ) + *sgngam = -1; + else + *sgngam = 1; + z = q - p; + if( z > 0.5 ) + { + p += 1.0; + z = p - q; + } + z = q * sin( PI * z ); + if( z == 0.0 ) + goto lgsing; +/* z = log(PI) - log( z ) - w;*/ + z = LOGPI - log( z ) - w; + return( z ); + } + +if( x < 13.0 ) + { + z = 1.0; + p = 0.0; + u = x; + while( u >= 3.0 ) + { + p -= 1.0; + u = x + p; + z *= u; + } + while( u < 2.0 ) + { + if( u == 0.0 ) + goto lgsing; + z /= u; + p += 1.0; + u = x + p; + } + if( z < 0.0 ) + { + *sgngam = -1; + z = -z; + } + else + *sgngam = 1; + if( u == 2.0 ) + return( log(z) ); + p -= 2.0; + x = x + p; + p = x * polevl( x, B, 5 ) / p1evl( x, C, 6); + return( log(z) + p ); + } + +if( x > MAXLGM ) + { + _SET_ERRNO(ERANGE); + mtherr( "lgamma", OVERFLOW ); +#ifdef INFINITIES + return( *sgngam * INFINITY ); +#else + return( *sgngam * MAXNUM ); +#endif + } + +q = ( x - 0.5 ) * log(x) - x + LS2PI; +if( x > 1.0e8 ) + return( q ); + +p = 1.0/(x*x); +if( x >= 1000.0 ) + q += (( 7.9365079365079365079365e-4 * p + - 2.7777777777777777777778e-3) *p + + 0.0833333333333333333333) / x; +else + q += polevl( p, A, 4 ) / x; +return( q ); +} + +/* This is the C99 version */ + +double lgamma(double x) +{ + int local_sgngam=0; + return (__lgamma_r(x, &local_sgngam)); +} diff --git a/winsup/mingw/mingwex/math/lgammaf.c b/winsup/mingw/mingwex/math/lgammaf.c new file mode 100644 index 000000000..20982f999 --- /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)); +} diff --git a/winsup/mingw/mingwex/math/lgammal.c b/winsup/mingw/mingwex/math/lgammal.c new file mode 100644 index 000000000..d2b306afd --- /dev/null +++ b/winsup/mingw/mingwex/math/lgammal.c @@ -0,0 +1,416 @@ +/* lgaml() + * + * Natural logarithm of gamma function + * + * + * + * SYNOPSIS: + * + * long double x, y, __lgammal_r(); + * int* sgngaml; + * y = __lgammal_r( x, sgngaml ); + * + * long double x, y, lgammal(); + * y = lgammal( 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 the variable referenced by sgngaml. + * + * For arguments greater than 33, the logarithm of the gamma + * function is approximated by the logarithmic version of + * Stirling's formula using a polynomial approximation of + * degree 4. Arguments between -33 and +33 are reduced by + * recurrence to the interval [2,3] of a rational approximation. + * The cosecant reflection formula is employed for arguments + * less than -33. + * + * Arguments greater than MAXLGML (10^4928) return MAXNUML. + * + * + * + * ACCURACY: + * + * + * arithmetic domain # trials peak rms + * IEEE -40, 40 100000 2.2e-19 4.6e-20 + * IEEE 10^-2000,10^+2000 20000 1.6e-19 3.3e-20 + * The error criterion was relative when the function magnitude + * was greater than one but absolute when it was less than one. + * + */ + +/* + * Copyright 1994 by Stephen L. Moshier + */ + +/* + * 26-11-2002 Modified for mingw. + * Danny Smith <dannysmith@users.sourceforge.net> + */ + +#ifndef __MINGW32__ +#include "mconf.h" +#ifdef ANSIPROT +extern long double fabsl ( long double ); +extern long double lgaml ( long double ); +extern long double logl ( long double ); +extern long double expl ( long double ); +extern long double gammal ( long double ); +extern long double sinl ( long double ); +extern long double floorl ( long double ); +extern long double powl ( long double, long double ); +extern long double polevll ( long double, void *, int ); +extern long double p1evll ( long double, void *, int ); +extern int isnanl ( long double ); +extern int isfinitel ( long double ); +#else +long double fabsl(), lgaml(), logl(), expl(), gammal(), sinl(); +long double floorl(), powl(), polevll(), p1evll(), isnanl(), isfinitel(); +#endif +#ifdef INFINITIES +extern long double INFINITYL; +#endif +#ifdef NANS +extern long double NANL; +#endif +#else /* __MINGW32__ */ +#include "cephes_mconf.h" +#endif /* __MINGW32__ */ + +#if UNK +static long double S[9] = { +-1.193945051381510095614E-3L, + 7.220599478036909672331E-3L, +-9.622023360406271645744E-3L, +-4.219773360705915470089E-2L, + 1.665386113720805206758E-1L, +-4.200263503403344054473E-2L, +-6.558780715202540684668E-1L, + 5.772156649015328608253E-1L, + 1.000000000000000000000E0L, +}; +#endif +#if IBMPC +static const unsigned short S[] = { +0xbaeb,0xd6d3,0x25e5,0x9c7e,0xbff5, XPD +0xfe9a,0xceb4,0xc74e,0xec9a,0x3ff7, XPD +0x9225,0xdfef,0xb0e9,0x9da5,0xbff8, XPD +0x10b0,0xec17,0x87dc,0xacd7,0xbffa, XPD +0x6b8d,0x7515,0x1905,0xaa89,0x3ffc, XPD +0xf183,0x126b,0xf47d,0xac0a,0xbffa, XPD +0x7bf6,0x57d1,0xa013,0xa7e7,0xbffe, XPD +0xc7a9,0x7db0,0x67e3,0x93c4,0x3ffe, XPD +0x0000,0x0000,0x0000,0x8000,0x3fff, XPD +}; +#endif +#if MIEEE +static long S[27] = { +0xbff50000,0x9c7e25e5,0xd6d3baeb, +0x3ff70000,0xec9ac74e,0xceb4fe9a, +0xbff80000,0x9da5b0e9,0xdfef9225, +0xbffa0000,0xacd787dc,0xec1710b0, +0x3ffc0000,0xaa891905,0x75156b8d, +0xbffa0000,0xac0af47d,0x126bf183, +0xbffe0000,0xa7e7a013,0x57d17bf6, +0x3ffe0000,0x93c467e3,0x7db0c7a9, +0x3fff0000,0x80000000,0x00000000, +}; +#endif + +#if UNK +static long double SN[9] = { + 1.133374167243894382010E-3L, + 7.220837261893170325704E-3L, + 9.621911155035976733706E-3L, +-4.219773343731191721664E-2L, +-1.665386113944413519335E-1L, +-4.200263503402112910504E-2L, + 6.558780715202536547116E-1L, + 5.772156649015328608727E-1L, +-1.000000000000000000000E0L, +}; +#endif +#if IBMPC +static const unsigned SN[] = { +0x5dd1,0x02de,0xb9f7,0x948d,0x3ff5, XPD +0x989b,0xdd68,0xc5f1,0xec9c,0x3ff7, XPD +0x2ca1,0x18f0,0x386f,0x9da5,0x3ff8, XPD +0x783f,0x41dd,0x87d1,0xacd7,0xbffa, XPD +0x7a5b,0xd76d,0x1905,0xaa89,0xbffc, XPD +0x7f64,0x1234,0xf47d,0xac0a,0xbffa, XPD +0x5e26,0x57d1,0xa013,0xa7e7,0x3ffe, XPD +0xc7aa,0x7db0,0x67e3,0x93c4,0x3ffe, XPD +0x0000,0x0000,0x0000,0x8000,0xbfff, XPD +}; +#endif +#if MIEEE +static long SN[27] = { +0x3ff50000,0x948db9f7,0x02de5dd1, +0x3ff70000,0xec9cc5f1,0xdd68989b, +0x3ff80000,0x9da5386f,0x18f02ca1, +0xbffa0000,0xacd787d1,0x41dd783f, +0xbffc0000,0xaa891905,0xd76d7a5b, +0xbffa0000,0xac0af47d,0x12347f64, +0x3ffe0000,0xa7e7a013,0x57d15e26, +0x3ffe0000,0x93c467e3,0x7db0c7aa, +0xbfff0000,0x80000000,0x00000000, +}; +#endif + + +/* A[]: Stirling's formula expansion of log gamma + * B[], C[]: log gamma function between 2 and 3 + */ + + +/* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x A(1/x^2) + * x >= 8 + * Peak relative error 1.51e-21 + * Relative spread of error peaks 5.67e-21 + */ +#if UNK +static long double A[7] = { + 4.885026142432270781165E-3L, +-1.880801938119376907179E-3L, + 8.412723297322498080632E-4L, +-5.952345851765688514613E-4L, + 7.936507795855070755671E-4L, +-2.777777777750349603440E-3L, + 8.333333333333331447505E-2L, +}; +#endif +#if IBMPC +static const unsigned short A[] = { +0xd984,0xcc08,0x91c2,0xa012,0x3ff7, XPD +0x3d91,0x0304,0x3da1,0xf685,0xbff5, XPD +0x3bdc,0xaad1,0xd492,0xdc88,0x3ff4, XPD +0x8b20,0x9fce,0x844e,0x9c09,0xbff4, XPD +0xf8f2,0x30e5,0x0092,0xd00d,0x3ff4, XPD +0x4d88,0x03a8,0x60b6,0xb60b,0xbff6, XPD +0x9fcc,0xaaaa,0xaaaa,0xaaaa,0x3ffb, XPD +}; +#endif +#if MIEEE +static long A[21] = { +0x3ff70000,0xa01291c2,0xcc08d984, +0xbff50000,0xf6853da1,0x03043d91, +0x3ff40000,0xdc88d492,0xaad13bdc, +0xbff40000,0x9c09844e,0x9fce8b20, +0x3ff40000,0xd00d0092,0x30e5f8f2, +0xbff60000,0xb60b60b6,0x03a84d88, +0x3ffb0000,0xaaaaaaaa,0xaaaa9fcc, +}; +#endif + +/* log gamma(x+2) = x B(x)/C(x) + * 0 <= x <= 1 + * Peak relative error 7.16e-22 + * Relative spread of error peaks 4.78e-20 + */ +#if UNK +static long double B[7] = { +-2.163690827643812857640E3L, +-8.723871522843511459790E4L, +-1.104326814691464261197E6L, +-6.111225012005214299996E6L, +-1.625568062543700591014E7L, +-2.003937418103815175475E7L, +-8.875666783650703802159E6L, +}; +static long double C[7] = { +/* 1.000000000000000000000E0L,*/ +-5.139481484435370143617E2L, +-3.403570840534304670537E4L, +-6.227441164066219501697E5L, +-4.814940379411882186630E6L, +-1.785433287045078156959E7L, +-3.138646407656182662088E7L, +-2.099336717757895876142E7L, +}; +#endif +#if IBMPC +static const unsigned short B[] = { +0x9557,0x4995,0x0da1,0x873b,0xc00a, XPD +0xfe44,0x9af8,0x5b8c,0xaa63,0xc00f, XPD +0x5aa8,0x7cf5,0x3684,0x86ce,0xc013, XPD +0x259a,0x258c,0xf206,0xba7f,0xc015, XPD +0xbe18,0x1ca3,0xc0a0,0xf80a,0xc016, XPD +0x168f,0x2c42,0x6717,0x98e3,0xc017, XPD +0x2051,0x9d55,0x92c8,0x876e,0xc016, XPD +}; +static const unsigned short C[] = { +/*0x0000,0x0000,0x0000,0x8000,0x3fff, XPD*/ +0xaa77,0xcf2f,0xae76,0x807c,0xc008, XPD +0xb280,0x0d74,0xb55a,0x84f3,0xc00e, XPD +0xa505,0xcd30,0x81dc,0x9809,0xc012, XPD +0x3369,0x4246,0xb8c2,0x92f0,0xc015, XPD +0x63cf,0x6aee,0xbe6f,0x8837,0xc017, XPD +0x26bb,0xccc7,0xb009,0xef75,0xc017, XPD +0x462b,0xbae8,0xab96,0xa02a,0xc017, XPD +}; +#endif +#if MIEEE +static long B[21] = { +0xc00a0000,0x873b0da1,0x49959557, +0xc00f0000,0xaa635b8c,0x9af8fe44, +0xc0130000,0x86ce3684,0x7cf55aa8, +0xc0150000,0xba7ff206,0x258c259a, +0xc0160000,0xf80ac0a0,0x1ca3be18, +0xc0170000,0x98e36717,0x2c42168f, +0xc0160000,0x876e92c8,0x9d552051, +}; +static long C[21] = { +/*0x3fff0000,0x80000000,0x00000000,*/ +0xc0080000,0x807cae76,0xcf2faa77, +0xc00e0000,0x84f3b55a,0x0d74b280, +0xc0120000,0x980981dc,0xcd30a505, +0xc0150000,0x92f0b8c2,0x42463369, +0xc0170000,0x8837be6f,0x6aee63cf, +0xc0170000,0xef75b009,0xccc726bb, +0xc0170000,0xa02aab96,0xbae8462b, +}; +#endif + +/* log( sqrt( 2*pi ) ) */ +static const long double LS2PI = 0.91893853320467274178L; +#define MAXLGM 1.04848146839019521116e+4928L + + +/* Logarithm of gamma function */ +/* Reentrant version */ + +long double __lgammal_r(long double x, int* sgngaml) +{ +long double p, q, w, z, f, nx; +int i; + +*sgngaml = 1; +#ifdef NANS +if( isnanl(x) ) + return(NANL); +#endif +#ifdef INFINITIES +if( !isfinitel(x) ) + return(INFINITYL); +#endif +if( x < -34.0L ) + { + q = -x; + w = __lgammal_r(q, sgngaml); /* note this modifies sgngam! */ + p = floorl(q); + if( p == q ) + { +lgsing: + _SET_ERRNO(EDOM); + mtherr( "lgammal", SING ); +#ifdef INFINITIES + return (INFINITYL); +#else + return (MAXNUML); +#endif + } + i = p; + if( (i & 1) == 0 ) + *sgngaml = -1; + else + *sgngaml = 1; + z = q - p; + if( z > 0.5L ) + { + p += 1.0L; + z = p - q; + } + z = q * sinl( PIL * z ); + if( z == 0.0L ) + goto lgsing; +/* z = LOGPI - logl( z ) - w; */ + z = logl( PIL/z ) - w; + return( z ); + } + +if( x < 13.0L ) + { + z = 1.0L; + nx = floorl( x + 0.5L ); + f = x - nx; + while( x >= 3.0L ) + { + nx -= 1.0L; + x = nx + f; + z *= x; + } + while( x < 2.0L ) + { + if( fabsl(x) <= 0.03125 ) + goto lsmall; + z /= nx + f; + nx += 1.0L; + x = nx + f; + } + if( z < 0.0L ) + { + *sgngaml = -1; + z = -z; + } + else + *sgngaml = 1; + if( x == 2.0L ) + return( logl(z) ); + x = (nx - 2.0L) + f; + p = x * polevll( x, B, 6 ) / p1evll( x, C, 7); + return( logl(z) + p ); + } + +if( x > MAXLGM ) + { + _SET_ERRNO(ERANGE); + mtherr( "lgammal", OVERFLOW ); +#ifdef INFINITIES + return( *sgngaml * INFINITYL ); +#else + return( *sgngaml * MAXNUML ); +#endif + } + +q = ( x - 0.5L ) * logl(x) - x + LS2PI; +if( x > 1.0e10L ) + return(q); +p = 1.0L/(x*x); +q += polevll( p, A, 6 ) / x; +return( q ); + + +lsmall: +if( x == 0.0L ) + goto lgsing; +if( x < 0.0L ) + { + x = -x; + q = z / (x * polevll( x, SN, 8 )); + } +else + q = z / (x * polevll( x, S, 8 )); +if( q < 0.0L ) + { + *sgngaml = -1; + q = -q; + } +else + *sgngaml = 1; +q = logl( q ); +return(q); +} + +/* This is the C99 version */ + +long double lgammal(long double x) +{ + int local_sgngaml=0; + return (__lgammal_r(x, &local_sgngaml)); +} diff --git a/winsup/mingw/mingwex/math/tgamma.c b/winsup/mingw/mingwex/math/tgamma.c new file mode 100644 index 000000000..c3912a890 --- /dev/null +++ b/winsup/mingw/mingwex/math/tgamma.c @@ -0,0 +1,385 @@ +/* gamma.c + * + * Gamma function + * + * + * + * SYNOPSIS: + * + * double x, y, __tgamma_r(); + * int* sgngam; + * y = __tgamma_r( x, sgngam ); + * + * double x, y, tgamma(); + * y = tgamma( 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 sgngam. + * + * Arguments |x| <= 34 are reduced by recurrence and the function + * approximated by a rational function of degree 6/7 in the + * interval (2,3). Large arguments are handled by Stirling's + * formula. Large negative arguments are made positive using + * a reflection formula. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -34, 34 10000 1.3e-16 2.5e-17 + * IEEE -170,-33 20000 2.3e-15 3.3e-16 + * IEEE -33, 33 20000 9.4e-16 2.2e-16 + * IEEE 33, 171.6 20000 2.3e-15 3.2e-16 + * + * Error for arguments outside the test range will be larger + * owing to error amplification by the exponential function. + * + */ + +/* +Cephes Math Library Release 2.8: June, 2000 +Copyright 1984, 1987, 1989, 1992, 2000 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 + +#ifdef UNK +static const double P[] = { + 1.60119522476751861407E-4, + 1.19135147006586384913E-3, + 1.04213797561761569935E-2, + 4.76367800457137231464E-2, + 2.07448227648435975150E-1, + 4.94214826801497100753E-1, + 9.99999999999999996796E-1 +}; +static const double Q[] = { +-2.31581873324120129819E-5, + 5.39605580493303397842E-4, +-4.45641913851797240494E-3, + 1.18139785222060435552E-2, + 3.58236398605498653373E-2, +-2.34591795718243348568E-1, + 7.14304917030273074085E-2, + 1.00000000000000000320E0 +}; +#define MAXGAM 171.624376956302725 +static const double LOGPI = 1.14472988584940017414; +#endif + +#ifdef DEC +static const unsigned short P[] = { +0035047,0162701,0146301,0005234, +0035634,0023437,0032065,0176530, +0036452,0137157,0047330,0122574, +0037103,0017310,0143041,0017232, +0037524,0066516,0162563,0164605, +0037775,0004671,0146237,0014222, +0040200,0000000,0000000,0000000 +}; +static const unsigned short Q[] = { +0134302,0041724,0020006,0116565, +0035415,0072121,0044251,0025634, +0136222,0003447,0035205,0121114, +0036501,0107552,0154335,0104271, +0037022,0135717,0014776,0171471, +0137560,0034324,0165024,0037021, +0037222,0045046,0047151,0161213, +0040200,0000000,0000000,0000000 +}; +#define MAXGAM 34.84425627277176174 +#endif + +#ifdef IBMPC +static const unsigned short P[] = { +0x2153,0x3998,0xfcb8,0x3f24, +0xbfab,0xe686,0x84e3,0x3f53, +0x14b0,0xe9db,0x57cd,0x3f85, +0x23d3,0x18c4,0x63d9,0x3fa8, +0x7d31,0xdcae,0x8da9,0x3fca, +0xe312,0x3993,0xa137,0x3fdf, +0x0000,0x0000,0x0000,0x3ff0 +}; +static const unsigned short Q[] = { +0xd3af,0x8400,0x487a,0xbef8, +0x2573,0x2915,0xae8a,0x3f41, +0xb44a,0xe750,0x40e4,0xbf72, +0xb117,0x5b1b,0x31ed,0x3f88, +0xde67,0xe33f,0x5779,0x3fa2, +0x87c2,0x9d42,0x071a,0xbfce, +0x3c51,0xc9cd,0x4944,0x3fb2, +0x0000,0x0000,0x0000,0x3ff0 +}; +#define MAXGAM 171.624376956302725 +#endif + +#ifdef MIEEE +static const unsigned short P[] = { +0x3f24,0xfcb8,0x3998,0x2153, +0x3f53,0x84e3,0xe686,0xbfab, +0x3f85,0x57cd,0xe9db,0x14b0, +0x3fa8,0x63d9,0x18c4,0x23d3, +0x3fca,0x8da9,0xdcae,0x7d31, +0x3fdf,0xa137,0x3993,0xe312, +0x3ff0,0x0000,0x0000,0x0000 +}; +static const unsigned short Q[] = { +0xbef8,0x487a,0x8400,0xd3af, +0x3f41,0xae8a,0x2915,0x2573, +0xbf72,0x40e4,0xe750,0xb44a, +0x3f88,0x31ed,0x5b1b,0xb117, +0x3fa2,0x5779,0xe33f,0xde67, +0xbfce,0x071a,0x9d42,0x87c2, +0x3fb2,0x4944,0xc9cd,0x3c51, +0x3ff0,0x0000,0x0000,0x0000 +}; +#define MAXGAM 171.624376956302725 +#endif + +/* Stirling's formula for the gamma function */ +#if UNK +static const double STIR[5] = { + 7.87311395793093628397E-4, +-2.29549961613378126380E-4, +-2.68132617805781232825E-3, + 3.47222221605458667310E-3, + 8.33333333333482257126E-2, +}; +#define MAXSTIR 143.01608 +static const double SQTPI = 2.50662827463100050242E0; +#endif +#if DEC +static const unsigned short STIR[20] = { +0035516,0061622,0144553,0112224, +0135160,0131531,0037460,0165740, +0136057,0134460,0037242,0077270, +0036143,0107070,0156306,0027751, +0037252,0125252,0125252,0146064, +}; +#define MAXSTIR 26.77 +static const unsigned short SQT[4] = { +0040440,0066230,0177661,0034055, +}; +#define SQTPI *(double *)SQT +#endif +#if IBMPC +static const unsigned short STIR[20] = { +0x7293,0x592d,0xcc72,0x3f49, +0x1d7c,0x27e6,0x166b,0xbf2e, +0x4fd7,0x07d4,0xf726,0xbf65, +0xc5fd,0x1b98,0x71c7,0x3f6c, +0x5986,0x5555,0x5555,0x3fb5, +}; +#define MAXSTIR 143.01608 +static const unsigned short SQT[4] = { +0x2706,0x1ff6,0x0d93,0x4004, +}; +#define SQTPI *(double *)SQT +#endif +#if MIEEE +static const unsigned short STIR[20] = { +0x3f49,0xcc72,0x592d,0x7293, +0xbf2e,0x166b,0x27e6,0x1d7c, +0xbf65,0xf726,0x07d4,0x4fd7, +0x3f6c,0x71c7,0x1b98,0xc5fd, +0x3fb5,0x5555,0x5555,0x5986, +}; +#define MAXSTIR 143.01608 +static const unsigned short SQT[4] = { +0x4004,0x0d93,0x1ff6,0x2706, +}; +#define SQTPI *(double *)SQT +#endif + +#ifndef __MINGW32__ +int sgngam = 0; +extern int sgngam; +extern double MAXLOG, MAXNUM, PI; +#ifdef ANSIPROT +extern double pow ( double, double ); +extern double log ( double ); +extern double exp ( double ); +extern double sin ( double ); +extern double polevl ( double, void *, int ); +extern double p1evl ( double, void *, int ); +extern double floor ( double ); +extern double fabs ( double ); +extern int isnan ( double ); +extern int isfinite ( double ); +static double stirf ( double ); +double lgam ( double ); +#else +double pow(), log(), exp(), sin(), polevl(), p1evl(), floor(), fabs(); +int isnan(), isfinite(); +static double stirf(); +double lgam(); +#endif +#ifdef INFINITIES +extern double INFINITY; +#endif +#ifdef NANS +extern double NAN; +#endif +#else /* __MINGW32__ */ +static double stirf ( double ); +#endif + +/* Gamma function computed by Stirling's formula. + * The polynomial STIR is valid for 33 <= x <= 172. + */ +static double stirf(x) +double x; +{ +double y, w, v; + +w = 1.0/x; +w = 1.0 + w * polevl( w, STIR, 4 ); +y = exp(x); +if( x > MAXSTIR ) + { /* Avoid overflow in pow() */ + v = pow( x, 0.5 * x - 0.25 ); + y = v * (v / y); + } +else + { + y = pow( x, x - 0.5 ) / y; + } +y = SQTPI * y * w; +return( y ); +} + + + +double __tgamma_r(double x, int* sgngam) +{ +double p, q, z; +int i; + +*sgngam = 1; +#ifdef NANS +if( isnan(x) ) + return(x); +#endif +#ifdef INFINITIES +#ifdef NANS +if( x == INFINITY ) + return(x); +if( x == -INFINITY ) + return(NAN); +#else +if( !isfinite(x) ) + return(x); +#endif +#endif +q = fabs(x); + +if( q > 33.0 ) + { + if( x < 0.0 ) + { + p = floor(q); + if( p == q ) + { +gsing: + _SET_ERRNO(EDOM); + mtherr( "tgamma", SING ); +#ifdef INFINITIES + return (INFINITY); +#else + return (MAXNUM); +#endif + } + i = p; + if( (i & 1) == 0 ) + *sgngam = -1; + z = q - p; + if( z > 0.5 ) + { + p += 1.0; + z = q - p; + } + z = q * sin( PI * z ); + if( z == 0.0 ) + { + _SET_ERRNO(ERANGE); + mtherr( "tgamma", OVERFLOW ); +#ifdef INFINITIES + return( *sgngam * INFINITY); +#else + return( *sgngam * MAXNUM); +#endif + } + z = fabs(z); + z = PI/(z * stirf(q) ); + } + else + { + z = stirf(x); + } + return( *sgngam * z ); + } + +z = 1.0; +while( x >= 3.0 ) + { + x -= 1.0; + z *= x; + } + +while( x < 0.0 ) + { + if( x > -1.E-9 ) + goto small; + z /= x; + x += 1.0; + } + +while( x < 2.0 ) + { + if( x < 1.e-9 ) + goto small; + z /= x; + x += 1.0; + } + +if( x == 2.0 ) + return(z); + +x -= 2.0; +p = polevl( x, P, 6 ); +q = polevl( x, Q, 7 ); +return( z * p / q ); + +small: +if( x == 0.0 ) + { + goto gsing; + } +else + return( z/((1.0 + 0.5772156649015329 * x) * x) ); +} + +/* This is the C99 version */ + +double tgamma(double x) +{ + int local_sgngam=0; + return (__tgamma_r(x, &local_sgngam)); +} diff --git a/winsup/mingw/mingwex/math/tgammaf.c b/winsup/mingw/mingwex/math/tgammaf.c new file mode 100644 index 000000000..38e7d75a1 --- /dev/null +++ b/winsup/mingw/mingwex/math/tgammaf.c @@ -0,0 +1,265 @@ +/* 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)); +} diff --git a/winsup/mingw/mingwex/math/tgammal.c b/winsup/mingw/mingwex/math/tgammal.c new file mode 100644 index 000000000..682a12e8e --- /dev/null +++ b/winsup/mingw/mingwex/math/tgammal.c @@ -0,0 +1,501 @@ +/* gammal.c + * + * Gamma function + * + * + * + * SYNOPSIS: + * + * long double x, y, __tgammal_r(); + * int* sgngaml; + * y = __tgammal_r( x, sgngaml ); + * + * long double x, y, tgammal(); + * y = tgammal( 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 |x| <= 13 are reduced by recurrence and the function + * approximated by a rational function of degree 7/8 in the + * interval (2,3). Large arguments are handled by Stirling's + * formula. Large negative arguments are made positive using + * a reflection formula. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -40,+40 10000 3.6e-19 7.9e-20 + * IEEE -1755,+1755 10000 4.8e-18 6.5e-19 + * + * Accuracy for large arguments is dominated by error in powl(). + * + */ + +/* +Copyright 1994 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 + +/* +gamma(x+2) = gamma(x+2) P(x)/Q(x) +0 <= x <= 1 +Relative error +n=7, d=8 +Peak error = 1.83e-20 +Relative error spread = 8.4e-23 +*/ + +#if UNK +static const long double P[8] = { + 4.212760487471622013093E-5L, + 4.542931960608009155600E-4L, + 4.092666828394035500949E-3L, + 2.385363243461108252554E-2L, + 1.113062816019361559013E-1L, + 3.629515436640239168939E-1L, + 8.378004301573126728826E-1L, + 1.000000000000000000009E0L, +}; +static const long double Q[9] = { +-1.397148517476170440917E-5L, + 2.346584059160635244282E-4L, +-1.237799246653152231188E-3L, +-7.955933682494738320586E-4L, + 2.773706565840072979165E-2L, +-4.633887671244534213831E-2L, +-2.243510905670329164562E-1L, + 4.150160950588455434583E-1L, + 9.999999999999999999908E-1L, +}; +#endif +#if IBMPC +static const short P[] = { +0x434a,0x3f22,0x2bda,0xb0b2,0x3ff0, XPD +0xf5aa,0xe82f,0x335b,0xee2e,0x3ff3, XPD +0xbe6c,0x3757,0xc717,0x861b,0x3ff7, XPD +0x7f43,0x5196,0xb166,0xc368,0x3ff9, XPD +0x9549,0x8eb5,0x8c3a,0xe3f4,0x3ffb, XPD +0x8d75,0x23af,0xc8e4,0xb9d4,0x3ffd, XPD +0x29cf,0x19b3,0x16c8,0xd67a,0x3ffe, XPD +0x0000,0x0000,0x0000,0x8000,0x3fff, XPD +}; +static const short Q[] = { +0x5473,0x2de8,0x1268,0xea67,0xbfee, XPD +0x334b,0xc2f0,0xa2dd,0xf60e,0x3ff2, XPD +0xbeed,0x1853,0xa691,0xa23d,0xbff5, XPD +0x296e,0x7cb1,0x5dfd,0xd08f,0xbff4, XPD +0x0417,0x7989,0xd7bc,0xe338,0x3ff9, XPD +0x3295,0x3698,0xd580,0xbdcd,0xbffa, XPD +0x75ef,0x3ab7,0x4ad3,0xe5bc,0xbffc, XPD +0xe458,0x2ec7,0xfd57,0xd47c,0x3ffd, XPD +0x0000,0x0000,0x0000,0x8000,0x3fff, XPD +}; +#endif +#if MIEEE +static const long P[24] = { +0x3ff00000,0xb0b22bda,0x3f22434a, +0x3ff30000,0xee2e335b,0xe82ff5aa, +0x3ff70000,0x861bc717,0x3757be6c, +0x3ff90000,0xc368b166,0x51967f43, +0x3ffb0000,0xe3f48c3a,0x8eb59549, +0x3ffd0000,0xb9d4c8e4,0x23af8d75, +0x3ffe0000,0xd67a16c8,0x19b329cf, +0x3fff0000,0x80000000,0x00000000, +}; +static const long Q[27] = { +0xbfee0000,0xea671268,0x2de85473, +0x3ff20000,0xf60ea2dd,0xc2f0334b, +0xbff50000,0xa23da691,0x1853beed, +0xbff40000,0xd08f5dfd,0x7cb1296e, +0x3ff90000,0xe338d7bc,0x79890417, +0xbffa0000,0xbdcdd580,0x36983295, +0xbffc0000,0xe5bc4ad3,0x3ab775ef, +0x3ffd0000,0xd47cfd57,0x2ec7e458, +0x3fff0000,0x80000000,0x00000000, +}; +#endif +/* +static const long double P[] = { +-3.01525602666895735709e0L, +-3.25157411956062339893e1L, +-2.92929976820724030353e2L, +-1.70730828800510297666e3L, +-7.96667499622741999770e3L, +-2.59780216007146401957e4L, +-5.99650230220855581642e4L, +-7.15743521530849602425e4L +}; +static const long double Q[] = { + 1.00000000000000000000e0L, +-1.67955233807178858919e1L, + 8.85946791747759881659e1L, + 5.69440799097468430177e1L, +-1.98526250512761318471e3L, + 3.31667508019495079814e3L, + 1.60577839621734713377e4L, +-2.97045081369399940529e4L, +-7.15743521530849602412e4L +}; +*/ +#define MAXGAML 1755.455L +/*static const long double LOGPI = 1.14472988584940017414L;*/ + +/* Stirling's formula for the gamma function +gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) +z(x) = x +13 <= x <= 1024 +Relative error +n=8, d=0 +Peak error = 9.44e-21 +Relative error spread = 8.8e-4 +*/ +#if UNK +static const long double STIR[9] = { + 7.147391378143610789273E-4L, +-2.363848809501759061727E-5L, +-5.950237554056330156018E-4L, + 6.989332260623193171870E-5L, + 7.840334842744753003862E-4L, +-2.294719747873185405699E-4L, +-2.681327161876304418288E-3L, + 3.472222222230075327854E-3L, + 8.333333333333331800504E-2L, +}; +#endif +#if IBMPC +static const short STIR[] = { +0x6ede,0x69f7,0x54e3,0xbb5d,0x3ff4, XPD +0xc395,0x0295,0x4443,0xc64b,0xbfef, XPD +0xba6f,0x7c59,0x5e47,0x9bfb,0xbff4, XPD +0x5704,0x1a39,0xb11d,0x9293,0x3ff1, XPD +0x30b7,0x1a21,0x98b2,0xcd87,0x3ff4, XPD +0xbef3,0x7023,0x6a08,0xf09e,0xbff2, XPD +0x3a1c,0x5ac8,0x3478,0xafb9,0xbff6, XPD +0xc3c9,0x906e,0x38e3,0xe38e,0x3ff6, XPD +0xa1d5,0xaaaa,0xaaaa,0xaaaa,0x3ffb, XPD +}; +#endif +#if MIEEE +static const long STIR[27] = { +0x3ff40000,0xbb5d54e3,0x69f76ede, +0xbfef0000,0xc64b4443,0x0295c395, +0xbff40000,0x9bfb5e47,0x7c59ba6f, +0x3ff10000,0x9293b11d,0x1a395704, +0x3ff40000,0xcd8798b2,0x1a2130b7, +0xbff20000,0xf09e6a08,0x7023bef3, +0xbff60000,0xafb93478,0x5ac83a1c, +0x3ff60000,0xe38e38e3,0x906ec3c9, +0x3ffb0000,0xaaaaaaaa,0xaaaaa1d5, +}; +#endif +#define MAXSTIR 1024.0L +static const long double SQTPI = 2.50662827463100050242E0L; + +/* 1/gamma(x) = z P(z) + * z(x) = 1/x + * 0 < x < 0.03125 + * Peak relative error 4.2e-23 + */ +#if UNK +static const long double S[9] = { +-1.193945051381510095614E-3L, + 7.220599478036909672331E-3L, +-9.622023360406271645744E-3L, +-4.219773360705915470089E-2L, + 1.665386113720805206758E-1L, +-4.200263503403344054473E-2L, +-6.558780715202540684668E-1L, + 5.772156649015328608253E-1L, + 1.000000000000000000000E0L, +}; +#endif +#if IBMPC +static const unsigned short S[] = { +0xbaeb,0xd6d3,0x25e5,0x9c7e,0xbff5, XPD +0xfe9a,0xceb4,0xc74e,0xec9a,0x3ff7, XPD +0x9225,0xdfef,0xb0e9,0x9da5,0xbff8, XPD +0x10b0,0xec17,0x87dc,0xacd7,0xbffa, XPD +0x6b8d,0x7515,0x1905,0xaa89,0x3ffc, XPD +0xf183,0x126b,0xf47d,0xac0a,0xbffa, XPD +0x7bf6,0x57d1,0xa013,0xa7e7,0xbffe, XPD +0xc7a9,0x7db0,0x67e3,0x93c4,0x3ffe, XPD +0x0000,0x0000,0x0000,0x8000,0x3fff, XPD +}; +#endif +#if MIEEE +static const long S[27] = { +0xbff50000,0x9c7e25e5,0xd6d3baeb, +0x3ff70000,0xec9ac74e,0xceb4fe9a, +0xbff80000,0x9da5b0e9,0xdfef9225, +0xbffa0000,0xacd787dc,0xec1710b0, +0x3ffc0000,0xaa891905,0x75156b8d, +0xbffa0000,0xac0af47d,0x126bf183, +0xbffe0000,0xa7e7a013,0x57d17bf6, +0x3ffe0000,0x93c467e3,0x7db0c7a9, +0x3fff0000,0x80000000,0x00000000, +}; +#endif +/* 1/gamma(-x) = z P(z) + * z(x) = 1/x + * 0 < x < 0.03125 + * Peak relative error 5.16e-23 + * Relative error spread = 2.5e-24 + */ +#if UNK +static const long double SN[9] = { + 1.133374167243894382010E-3L, + 7.220837261893170325704E-3L, + 9.621911155035976733706E-3L, +-4.219773343731191721664E-2L, +-1.665386113944413519335E-1L, +-4.200263503402112910504E-2L, + 6.558780715202536547116E-1L, + 5.772156649015328608727E-1L, +-1.000000000000000000000E0L, +}; +#endif +#if IBMPC +static const unsigned short SN[] = { +0x5dd1,0x02de,0xb9f7,0x948d,0x3ff5, XPD +0x989b,0xdd68,0xc5f1,0xec9c,0x3ff7, XPD +0x2ca1,0x18f0,0x386f,0x9da5,0x3ff8, XPD +0x783f,0x41dd,0x87d1,0xacd7,0xbffa, XPD +0x7a5b,0xd76d,0x1905,0xaa89,0xbffc, XPD +0x7f64,0x1234,0xf47d,0xac0a,0xbffa, XPD +0x5e26,0x57d1,0xa013,0xa7e7,0x3ffe, XPD +0xc7aa,0x7db0,0x67e3,0x93c4,0x3ffe, XPD +0x0000,0x0000,0x0000,0x8000,0xbfff, XPD +}; +#endif +#if MIEEE +static const long SN[27] = { +0x3ff50000,0x948db9f7,0x02de5dd1, +0x3ff70000,0xec9cc5f1,0xdd68989b, +0x3ff80000,0x9da5386f,0x18f02ca1, +0xbffa0000,0xacd787d1,0x41dd783f, +0xbffc0000,0xaa891905,0xd76d7a5b, +0xbffa0000,0xac0af47d,0x12347f64, +0x3ffe0000,0xa7e7a013,0x57d15e26, +0x3ffe0000,0x93c467e3,0x7db0c7aa, +0xbfff0000,0x80000000,0x00000000, +}; +#endif + +#ifndef __MINGW32__ +extern long double MAXLOGL, MAXNUML, PIL; +/* #define PIL 3.14159265358979323846L */ +/* #define MAXNUML 1.189731495357231765021263853E4932L */ + +#ifdef ANSIPROT +extern long double fabsl ( long double ); +extern long double lgaml ( long double ); +extern long double logl ( long double ); +extern long double expl ( long double ); +extern long double gammal ( long double ); +extern long double sinl ( long double ); +extern long double floorl ( long double ); +extern long double powl ( long double, long double ); +extern long double polevll ( long double, void *, int ); +extern long double p1evll ( long double, void *, int ); +extern int isnanl ( long double ); +extern int isfinitel ( long double ); +static long double stirf ( long double ); +#else +long double fabsl(), lgaml(), logl(), expl(), gammal(), sinl(); +long double floorl(), powl(), polevll(), p1evll(), isnanl(), isfinitel(); +static long double stirf(); +#endif +#ifdef INFINITIES +extern long double INFINITYL; +#endif +#ifdef NANS +extern long double NANL; +#endif + +#else /* __MINGW32__ */ +static long double stirf ( long double ); +#endif + + +/* Gamma function computed by Stirling's formula. */ + +static long double stirf(x) +long double x; +{ +long double y, w, v; + +w = 1.0L/x; +/* For large x, use rational coefficients from the analytical expansion. */ +if( x > 1024.0L ) + w = (((((6.97281375836585777429E-5L * w + + 7.84039221720066627474E-4L) * w + - 2.29472093621399176955E-4L) * w + - 2.68132716049382716049E-3L) * w + + 3.47222222222222222222E-3L) * w + + 8.33333333333333333333E-2L) * w + + 1.0L; +else + w = 1.0L + w * polevll( w, STIR, 8 ); +y = expl(x); +if( x > MAXSTIR ) + { /* Avoid overflow in pow() */ + v = powl( x, 0.5L * x - 0.25L ); + y = v * (v / y); + } +else + { + y = powl( x, x - 0.5L ) / y; + } +y = SQTPI * y * w; +return( y ); +} + + +long double __tgammal_r(long double x, int* sgngaml) +{ +long double p, q, z; +int i; + +*sgngaml = 1; +#ifdef NANS +if( isnanl(x) ) + return(NANL); +#endif +#ifdef INFINITIES +#ifdef NANS +if( x == INFINITYL ) + return(x); +if( x == -INFINITYL ) + return(NANL); +#else +if( !isfinite(x) ) + return(x); +#endif +#endif +q = fabsl(x); + +if( q > 13.0L ) + { + if( q > MAXGAML ) + goto goverf; + if( x < 0.0L ) + { + p = floorl(q); + if( p == q ) + { +gsing: + _SET_ERRNO(EDOM); + mtherr( "tgammal", SING ); +#ifdef INFINITIES + return (INFINITYL); +#else + return( *sgngaml * MAXNUML); +#endif + } + i = p; + if( (i & 1) == 0 ) + *sgngaml = -1; + z = q - p; + if( z > 0.5L ) + { + p += 1.0L; + z = q - p; + } + z = q * sinl( PIL * z ); + z = fabsl(z) * stirf(q); + if( z <= PIL/MAXNUML ) + { +goverf: + _SET_ERRNO(ERANGE); + mtherr( "tgammal", OVERFLOW ); +#ifdef INFINITIES + return( *sgngaml * INFINITYL); +#else + return( *sgngaml * MAXNUML); +#endif + } + z = PIL/z; + } + else + { + z = stirf(x); + } + return( *sgngaml * z ); + } + +z = 1.0L; +while( x >= 3.0L ) + { + x -= 1.0L; + z *= x; + } + +while( x < -0.03125L ) + { + z /= x; + x += 1.0L; + } + +if( x <= 0.03125L ) + goto small; + +while( x < 2.0L ) + { + z /= x; + x += 1.0L; + } + +if( x == 2.0L ) + return(z); + +x -= 2.0L; +p = polevll( x, P, 7 ); +q = polevll( x, Q, 8 ); +return( z * p / q ); + +small: +if( x == 0.0L ) + { + goto gsing; + } +else + { + if( x < 0.0L ) + { + x = -x; + q = z / (x * polevll( x, SN, 8 )); + } + else + q = z / (x * polevll( x, S, 8 )); + } +return q; +} + + +/* This is the C99 version. */ + +long double tgammal(long double x) +{ + int local_sgngaml=0; + return (__tgammal_r(x, &local_sgngaml)); +} + |