#include #include "mconf.h" double p1evl(double, double [], int); double polevl(double, double [], int); #ifndef ANSIPROT double lgam(), exp(), log(), fabs(), igam(), igamc(); #endif int merror = 0; /* Notice: the order of appearance of the following * messages is bound to the error codes defined * in mconf.h. */ static char *ermsg[7] = { "unknown", /* error code 0 */ "domain", /* error code 1 */ "singularity", /* et seq. */ "overflow", "underflow", "total loss of precision", "partial loss of precision" }; /* const.c * * Globally declared constants * * * * SYNOPSIS: * * extern double nameofconstant; * * * * * DESCRIPTION: * * This file contains a number of mathematical constants and * also some needed size parameters of the computer arithmetic. * The values are supplied as arrays of hexadecimal integers * for IEEE arithmetic; arrays of octal constants for DEC * arithmetic; and in a normal decimal scientific notation for * other machines. The particular notation used is determined * by a symbol (DEC, IBMPC, or UNK) defined in the include file * mconf.h. * * The default size parameters are as follows. * * For DEC and UNK modes: * MACHEP = 1.38777878078144567553E-17 2**-56 * MAXLOG = 8.8029691931113054295988E1 log(2**127) * MINLOG = -8.872283911167299960540E1 log(2**-128) * MAXNUM = 1.701411834604692317316873e38 2**127 * * For IEEE arithmetic (IBMPC): * MACHEP = 1.11022302462515654042E-16 2**-53 * MAXLOG = 7.09782712893383996843E2 log(2**1024) * MINLOG = -7.08396418532264106224E2 log(2**-1022) * MAXNUM = 1.7976931348623158E308 2**1024 * * The global symbols for mathematical constants are * PI = 3.14159265358979323846 pi * PIO2 = 1.57079632679489661923 pi/2 * PIO4 = 7.85398163397448309616E-1 pi/4 * SQRT2 = 1.41421356237309504880 sqrt(2) * SQRTH = 7.07106781186547524401E-1 sqrt(2)/2 * LOG2E = 1.4426950408889634073599 1/log(2) * SQ2OPI = 7.9788456080286535587989E-1 sqrt( 2/pi ) * LOGE2 = 6.93147180559945309417E-1 log(2) * LOGSQ2 = 3.46573590279972654709E-1 log(2)/2 * THPIO4 = 2.35619449019234492885 3*pi/4 * TWOOPI = 6.36619772367581343075535E-1 2/pi * * These lists are subject to change. */ /* const.c */ /* Cephes Math Library Release 2.3: March, 1995 Copyright 1984, 1995 by Stephen L. Moshier */ #ifdef UNK #if 1 double MACHEP = 1.11022302462515654042E-16; /* 2**-53 */ #else double MACHEP = 1.38777878078144567553E-17; /* 2**-56 */ #endif double UFLOWTHRESH = 2.22507385850720138309E-308; /* 2**-1022 */ #ifdef DENORMAL double MAXLOG = 7.09782712893383996732E2; /* log(MAXNUM) */ /* double MINLOG = -7.44440071921381262314E2; */ /* log(2**-1074) */ double MINLOG = -7.451332191019412076235E2; /* log(2**-1075) */ #else double MAXLOG = 7.08396418532264106224E2; /* log 2**1022 */ double MINLOG = -7.08396418532264106224E2; /* log 2**-1022 */ #endif double MAXNUM = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */ double PI = 3.14159265358979323846; /* pi */ double PIO2 = 1.57079632679489661923; /* pi/2 */ double PIO4 = 7.85398163397448309616E-1; /* pi/4 */ double SQRT2 = 1.41421356237309504880; /* sqrt(2) */ double SQRTH = 7.07106781186547524401E-1; /* sqrt(2)/2 */ double LOG2E = 1.4426950408889634073599; /* 1/log(2) */ double SQ2OPI = 7.9788456080286535587989E-1; /* sqrt( 2/pi ) */ double LOGE2 = 6.93147180559945309417E-1; /* log(2) */ double LOGSQ2 = 3.46573590279972654709E-1; /* log(2)/2 */ double THPIO4 = 2.35619449019234492885; /* 3*pi/4 */ double TWOOPI = 6.36619772367581343075535E-1; /* 2/pi */ #ifdef INFINITIES double INFINITY = 1.0/0.0; /* 99e999; */ #else double INFINITY = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */ #endif #ifdef NANS double NAN = 1.0/0.0 - 1.0/0.0; #else double NAN = 0.0; #endif #ifdef MINUSZERO double NEGZERO = -0.0; #else double NEGZERO = 0.0; #endif #endif #ifdef IBMPC /* 2**-53 = 1.11022302462515654042E-16 */ unsigned short MACHEP[4] = {0x0000,0x0000,0x0000,0x3ca0}; unsigned short UFLOWTHRESH[4] = {0x0000,0x0000,0x0000,0x0010}; #ifdef DENORMAL /* log(MAXNUM) = 7.09782712893383996732224E2 */ unsigned short MAXLOG[4] = {0x39ef,0xfefa,0x2e42,0x4086}; /* log(2**-1074) = - -7.44440071921381262314E2 */ /*unsigned short MINLOG[4] = {0x71c3,0x446d,0x4385,0xc087};*/ unsigned short MINLOG[4] = {0x3052,0xd52d,0x4910,0xc087}; #else /* log(2**1022) = 7.08396418532264106224E2 */ unsigned short MAXLOG[4] = {0xbcd2,0xdd7a,0x232b,0x4086}; /* log(2**-1022) = - 7.08396418532264106224E2 */ unsigned short MINLOG[4] = {0xbcd2,0xdd7a,0x232b,0xc086}; #endif /* 2**1024*(1-MACHEP) = 1.7976931348623158E308 */ unsigned short MAXNUM[4] = {0xffff,0xffff,0xffff,0x7fef}; unsigned short PI[4] = {0x2d18,0x5444,0x21fb,0x4009}; unsigned short PIO2[4] = {0x2d18,0x5444,0x21fb,0x3ff9}; unsigned short PIO4[4] = {0x2d18,0x5444,0x21fb,0x3fe9}; unsigned short SQRT2[4] = {0x3bcd,0x667f,0xa09e,0x3ff6}; unsigned short SQRTH[4] = {0x3bcd,0x667f,0xa09e,0x3fe6}; unsigned short LOG2E[4] = {0x82fe,0x652b,0x1547,0x3ff7}; unsigned short SQ2OPI[4] = {0x3651,0x33d4,0x8845,0x3fe9}; unsigned short LOGE2[4] = {0x39ef,0xfefa,0x2e42,0x3fe6}; unsigned short LOGSQ2[4] = {0x39ef,0xfefa,0x2e42,0x3fd6}; unsigned short THPIO4[4] = {0x21d2,0x7f33,0xd97c,0x4002}; unsigned short TWOOPI[4] = {0xc883,0x6dc9,0x5f30,0x3fe4}; #ifdef INFINITIES unsigned short INFINITY[4] = {0x0000,0x0000,0x0000,0x7ff0}; #else unsigned short INFINITY[4] = {0xffff,0xffff,0xffff,0x7fef}; #endif #ifdef NANS unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x7ffc}; #else unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x0000}; #endif #ifdef MINUSZERO unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x8000}; #else unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x0000}; #endif #endif #ifdef MIEEE /* 2**-53 = 1.11022302462515654042E-16 */ unsigned short MACHEP[4] = {0x3ca0,0x0000,0x0000,0x0000}; unsigned short UFLOWTHRESH[4] = {0x0010,0x0000,0x0000,0x0000}; #ifdef DENORMAL /* log(2**1024) = 7.09782712893383996843E2 */ unsigned short MAXLOG[4] = {0x4086,0x2e42,0xfefa,0x39ef}; /* log(2**-1074) = - -7.44440071921381262314E2 */ /* unsigned short MINLOG[4] = {0xc087,0x4385,0x446d,0x71c3}; */ unsigned short MINLOG[4] = {0xc087,0x4910,0xd52d,0x3052}; #else /* log(2**1022) = 7.08396418532264106224E2 */ unsigned short MAXLOG[4] = {0x4086,0x232b,0xdd7a,0xbcd2}; /* log(2**-1022) = - 7.08396418532264106224E2 */ unsigned short MINLOG[4] = {0xc086,0x232b,0xdd7a,0xbcd2}; #endif /* 2**1024*(1-MACHEP) = 1.7976931348623158E308 */ unsigned short MAXNUM[4] = {0x7fef,0xffff,0xffff,0xffff}; unsigned short PI[4] = {0x4009,0x21fb,0x5444,0x2d18}; unsigned short PIO2[4] = {0x3ff9,0x21fb,0x5444,0x2d18}; unsigned short PIO4[4] = {0x3fe9,0x21fb,0x5444,0x2d18}; unsigned short SQRT2[4] = {0x3ff6,0xa09e,0x667f,0x3bcd}; unsigned short SQRTH[4] = {0x3fe6,0xa09e,0x667f,0x3bcd}; unsigned short LOG2E[4] = {0x3ff7,0x1547,0x652b,0x82fe}; unsigned short SQ2OPI[4] = {0x3fe9,0x8845,0x33d4,0x3651}; unsigned short LOGE2[4] = {0x3fe6,0x2e42,0xfefa,0x39ef}; unsigned short LOGSQ2[4] = {0x3fd6,0x2e42,0xfefa,0x39ef}; unsigned short THPIO4[4] = {0x4002,0xd97c,0x7f33,0x21d2}; unsigned short TWOOPI[4] = {0x3fe4,0x5f30,0x6dc9,0xc883}; #ifdef INFINITIES unsigned short INFINITY[4] = {0x7ff0,0x0000,0x0000,0x0000}; #else unsigned short INFINITY[4] = {0x7fef,0xffff,0xffff,0xffff}; #endif #ifdef NANS unsigned short NAN[4] = {0x7ff8,0x0000,0x0000,0x0000}; #else unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x0000}; #endif #ifdef MINUSZERO unsigned short NEGZERO[4] = {0x8000,0x0000,0x0000,0x0000}; #else unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x0000}; #endif #endif #ifdef DEC /* 2**-56 = 1.38777878078144567553E-17 */ unsigned short MACHEP[4] = {0022200,0000000,0000000,0000000}; unsigned short UFLOWTHRESH[4] = {0x0080,0x0000,0x0000,0x0000}; /* log 2**127 = 88.029691931113054295988 */ unsigned short MAXLOG[4] = {041660,007463,0143742,025733,}; /* log 2**-128 = -88.72283911167299960540 */ unsigned short MINLOG[4] = {0141661,071027,0173721,0147572,}; /* 2**127 = 1.701411834604692317316873e38 */ unsigned short MAXNUM[4] = {077777,0177777,0177777,0177777,}; unsigned short PI[4] = {040511,007732,0121041,064302,}; unsigned short PIO2[4] = {040311,007732,0121041,064302,}; unsigned short PIO4[4] = {040111,007732,0121041,064302,}; unsigned short SQRT2[4] = {040265,002363,031771,0157145,}; unsigned short SQRTH[4] = {040065,002363,031771,0157144,}; unsigned short LOG2E[4] = {040270,0125073,024534,013761,}; unsigned short SQ2OPI[4] = {040114,041051,0117241,0131204,}; unsigned short LOGE2[4] = {040061,071027,0173721,0147572,}; unsigned short LOGSQ2[4] = {037661,071027,0173721,0147572,}; unsigned short THPIO4[4] = {040426,0145743,0174631,007222,}; unsigned short TWOOPI[4] = {040042,0174603,067116,042025,}; /* Approximate infinity by MAXNUM. */ unsigned short INFINITY[4] = {077777,0177777,0177777,0177777,}; unsigned short NAN[4] = {0000000,0000000,0000000,0000000}; #ifdef MINUSZERO unsigned short NEGZERO[4] = {0000000,0000000,0000000,0100000}; #else unsigned short NEGZERO[4] = {0000000,0000000,0000000,0000000}; #endif #endif #ifndef UNK extern unsigned short MACHEP[]; extern unsigned short UFLOWTHRESH[]; extern unsigned short MAXLOG[]; extern unsigned short UNDLOG[]; extern unsigned short MINLOG[]; extern unsigned short MAXNUM[]; extern unsigned short PI[]; extern unsigned short PIO2[]; extern unsigned short PIO4[]; extern unsigned short SQRT2[]; extern unsigned short SQRTH[]; extern unsigned short LOG2E[]; extern unsigned short SQ2OPI[]; extern unsigned short LOGE2[]; extern unsigned short LOGSQ2[]; extern unsigned short THPIO4[]; extern unsigned short TWOOPI[]; extern unsigned short INFINITY[]; extern unsigned short NAN[]; extern unsigned short NEGZERO[]; #endif extern double MACHEP, MAXLOG; static double big = 4.503599627370496e15; static double biginv = 2.22044604925031308085e-16; #ifdef UNK static double P[] = { 1.60119522476751861407E-4, 1.19135147006586384913E-3, 1.04213797561761569935E-2, 4.76367800457137231464E-2, 2.07448227648435975150E-1, 4.94214826801497100753E-1, 9.99999999999999996796E-1 }; static 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 double LOGPI = 1.14472988584940017414; #endif #ifdef DEC static 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 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 static unsigned short LPI[4] = { 0040222,0103202,0043475,0006750, }; #define LOGPI *(double *)LPI #endif #ifdef IBMPC static 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 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 static unsigned short LPI[4] = { 0xa1bd,0x48e7,0x50d0,0x3ff2, }; #define LOGPI *(double *)LPI #endif #ifdef MIEEE static 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 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 static unsigned short LPI[4] = { 0x3ff2,0x50d0,0x48e7,0xa1bd, }; #define LOGPI *(double *)LPI #endif /* Stirling's formula for the gamma function */ #if UNK static double STIR[5] = { 7.87311395793093628397E-4, -2.29549961613378126380E-4, -2.68132617805781232825E-3, 3.47222221605458667310E-3, 8.33333333333482257126E-2, }; #define MAXSTIR 143.01608 static double SQTPI = 2.50662827463100050242E0; #endif #if DEC static 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 unsigned short SQT[4] = { 0040440,0066230,0177661,0034055, }; #define SQTPI *(double *)SQT #endif #if IBMPC static 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 unsigned short SQT[4] = { 0x2706,0x1ff6,0x0d93,0x4004, }; #define SQTPI *(double *)SQT #endif #if MIEEE static 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 unsigned short SQT[4] = { 0x4004,0x0d93,0x1ff6,0x2706, }; #define SQTPI *(double *)SQT #endif int sgngam = 0; extern int sgngam; extern double MAXLOG, MAXNUM, PI; #ifndef ANSIPROT 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 /* igam.c * * Incomplete gamma integral * * * * SYNOPSIS: * * double a, x, y, igam(); * * y = igam( a, x ); * * DESCRIPTION: * * The function is defined by * * x * - * 1 | | -t a-1 * igam(a,x) = ----- | e t dt. * - | | * | (a) - * 0 * * * In this implementation both arguments must be positive. * The integral is evaluated by either a power series or * continued fraction expansion, depending on the relative * values of a and x. * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,30 200000 3.6e-14 2.9e-15 * IEEE 0,100 300000 9.9e-14 1.5e-14 */ /* igamc() * * Complemented incomplete gamma integral * * * * SYNOPSIS: * * double a, x, y, igamc(); * * y = igamc( a, x ); * * DESCRIPTION: * * The function is defined by * * * igamc(a,x) = 1 - igam(a,x) * * inf. * - * 1 | | -t a-1 * = ----- | e t dt. * - | | * | (a) - * x * * * In this implementation both arguments must be positive. * The integral is evaluated by either a power series or * continued fraction expansion, depending on the relative * values of a and x. * * ACCURACY: * * Tested at random a, x. * a x Relative error: * arithmetic domain domain # trials peak rms * IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15 * IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15 */ /* Cephes Math Library Release 2.0: April, 1987 Copyright 1985, 1987 by Stephen L. Moshier Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */ double igamc( a, x ) double a, x; { double ans, ax, c, yc, r, t, y, z; double pk, pkm1, pkm2, qk, qkm1, qkm2; if( (x <= 0) || ( a <= 0) ) return( 1.0 ); if( (x < 1.0) || (x < a) ) return( 1.e0 - igam(a,x) ); ax = a * log(x) - x - lgam(a); if( ax < -MAXLOG ) { mtherr( "igamc", UNDERFLOW ); return( 0.0 ); } ax = exp(ax); /* continued fraction */ y = 1.0 - a; z = x + y + 1.0; c = 0.0; pkm2 = 1.0; qkm2 = x; pkm1 = x + 1.0; qkm1 = z * x; ans = pkm1/qkm1; do { c += 1.0; y += 1.0; z += 2.0; yc = y * c; pk = pkm1 * z - pkm2 * yc; qk = qkm1 * z - qkm2 * yc; if( qk != 0 ) { r = pk/qk; t = fabs( (ans - r)/r ); ans = r; } else t = 1.0; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; if( fabs(pk) > big ) { pkm2 *= biginv; pkm1 *= biginv; qkm2 *= biginv; qkm1 *= biginv; } } while( t > MACHEP ); return( ans * ax ); } /* left tail of incomplete gamma function: * * inf. k * a -x - x * x e > ---------- * - - * k=0 | (a+k+1) * */ double igam( a, x ) double a, x; { double ans, ax, c, r; if( (x <= 0) || ( a <= 0) ) return( 0.0 ); if( (x > 1.0) && (x > a ) ) return( 1.e0 - igamc(a,x) ); /* Compute x**a * exp(-x) / gamma(a) */ ax = a * log(x) - x - lgam(a); if( ax < -MAXLOG ) { mtherr( "igam", UNDERFLOW ); return( 0.0 ); } ax = exp(ax); /* power series */ r = a; c = 1.0; ans = 1.0; do { r += 1.0; c *= x/r; ans += c; } while( c/ans > MACHEP ); return( ans * ax/a ); } /* gamma.c * * Gamma function * * * * SYNOPSIS: * * double x, y, gamma(); * extern int sgngam; * * y = gamma( x ); * * * * DESCRIPTION: * * Returns gamma function of the argument. The result is * correctly signed, and the sign (+1 or -1) is also * returned in a global (extern) variable named sgngam. * This variable is also filled in by the logarithmic gamma * function lgam(). * * 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. * */ /* lgam() * * Natural logarithm of gamma function * * * * SYNOPSIS: * * double x, y, lgam(); * extern int sgngam; * * y = lgam( x ); * * * * DESCRIPTION: * * Returns the base e (2.718...) logarithm of the absolute * value of the gamma function of the argument. * The sign (+1 or -1) of the gamma function is returned in a * global (extern) variable named 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 * */ /* gamma.c */ /* gamma function */ /* Cephes Math Library Release 2.2: July, 1992 Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */ /* 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 gamma(x) double x; { 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 ) { #ifdef NANS gamnan: mtherr( "gamma", DOMAIN ); return (NAN); #else goto goverf; #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 ) { #ifdef INFINITIES return( sgngam * INFINITY); #else goverf: mtherr( "gamma", OVERFLOW ); 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 ) { #ifdef INFINITIES #ifdef NANS goto gamnan; #else return( INFINITY ); #endif #else mtherr( "gamma", SING ); return( MAXNUM ); #endif } else return( z/((1.0 + 0.5772156649015329 * x) * x) ); } /* 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 #endif #ifdef DEC static unsigned short A[] = { 0035524,0141201,0034633,0031405, 0135433,0176755,0126007,0045030, 0035520,0006371,0003342,0172730, 0136066,0005540,0132605,0026407, 0037252,0125252,0125252,0125132 }; static 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 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 unsigned short LS2P[] = {040153,037616,041445,0172645,}; #define LS2PI *(double *)LS2P #define MAXLGM 2.035093e36 #endif #ifdef IBMPC static unsigned short A[] = { 0x6661,0x2733,0x9850,0x3f4a, 0xe943,0xb580,0x7fbd,0xbf43, 0x5ebb,0x20dc,0x019f,0x3f4a, 0xa5a1,0x16b0,0xc16c,0xbf66, 0x554b,0x5555,0x5555,0x3fb5 }; static 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 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 unsigned short LS2P[] = { 0xbeb5,0xc864,0x67f1,0x3fed }; #define LS2PI *(double *)LS2P #define MAXLGM 2.556348e305 #endif #ifdef MIEEE static unsigned short A[] = { 0x3f4a,0x9850,0x2733,0x6661, 0xbf43,0x7fbd,0xb580,0xe943, 0x3f4a,0x019f,0x20dc,0x5ebb, 0xbf66,0xc16c,0x16b0,0xa5a1, 0x3fb5,0x5555,0x5555,0x554b }; static 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 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 unsigned short LS2P[] = { 0x3fed,0x67f1,0xc864,0xbeb5 }; #define LS2PI *(double *)LS2P #define MAXLGM 2.556348e305 #endif /* Logarithm of gamma function */ double lgam(x) double x; { 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 = lgam(q); /* note this modifies sgngam! */ p = floor(q); if( p == q ) { lgsing: #ifdef INFINITIES mtherr( "lgam", SING ); return (INFINITY); #else goto loverf; #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 ) { #ifdef INFINITIES return( sgngam * INFINITY ); #else loverf: mtherr( "lgam", OVERFLOW ); 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 ); } /* mtherr.c * * Library common error handling routine * * * * SYNOPSIS: * * char *fctnam; * int code; * int mtherr(); * * mtherr( fctnam, code ); * * * * DESCRIPTION: * * This routine may be called to report one of the following * error conditions (in the include file mconf.h). * * Mnemonic Value Significance * * DOMAIN 1 argument domain error * SING 2 function singularity * OVERFLOW 3 overflow range error * UNDERFLOW 4 underflow range error * TLOSS 5 total loss of precision * PLOSS 6 partial loss of precision * EDOM 33 Unix domain error code * ERANGE 34 Unix range error code * * The default version of the file prints the function name, * passed to it by the pointer fctnam, followed by the * error condition. The display is directed to the standard * output device. The routine then returns to the calling * program. Users may wish to modify the program to abort by * calling exit() under severe error conditions such as domain * errors. * * Since all error conditions pass control to this function, * the display may be easily changed, eliminated, or directed * to an error logging device. * * SEE ALSO: * * mconf.h * */ /* Cephes Math Library Release 2.0: April, 1987 Copyright 1984, 1987 by Stephen L. Moshier Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */ int mtherr( name, code ) char *name; int code; { /* Display string passed by calling program, * which is supposed to be the name of the * function in which the error occurred: */ printf("\n%s ", name); /* Set global error message word */ merror = code; /* Display error message defined * by the code argument. */ if( (code <= 0) || (code >= 7) ) code = 0; printf( "%s error\n", ermsg[code] ); /* Return to calling * program */ return( 0 ); } /* polevl.c * p1evl.c * * Evaluate polynomial * * * * SYNOPSIS: * * int N; * double x, y, coef[N+1], polevl[]; * * y = polevl( 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 */ double polevl( x, coef, N ) double x; double coef[]; int N; { double ans; int i; double *p; p = coef; ans = *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. */ double p1evl( x, coef, N ) double x; double coef[]; int N; { double ans; double *p; int i; p = coef; ans = x + *p++; i = N-1; do ans = ans * x + *p++; while( --i ); return( ans ); }