diff options
Diffstat (limited to 'nist/cephes.c')
| -rw-r--r-- | nist/cephes.c | 1327 |
1 files changed, 1327 insertions, 0 deletions
diff --git a/nist/cephes.c b/nist/cephes.c new file mode 100644 index 0000000..6fc1500 --- /dev/null +++ b/nist/cephes.c @@ -0,0 +1,1327 @@ +#include <stdio.h> +#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 ); +} |