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