diff options
Diffstat (limited to 'libm/double/rgamma.c')
| -rw-r--r-- | libm/double/rgamma.c | 209 |
1 files changed, 209 insertions, 0 deletions
diff --git a/libm/double/rgamma.c b/libm/double/rgamma.c new file mode 100644 index 000000000..1d6ff3840 --- /dev/null +++ b/libm/double/rgamma.c @@ -0,0 +1,209 @@ +/* rgamma.c + * + * Reciprocal gamma function + * + * + * + * SYNOPSIS: + * + * double x, y, rgamma(); + * + * y = rgamma( x ); + * + * + * + * DESCRIPTION: + * + * Returns one divided by the gamma function of the argument. + * + * The function is approximated by a Chebyshev expansion in + * the interval [0,1]. Range reduction is by recurrence + * for arguments between -34.034 and +34.84425627277176174. + * 1/MAXNUM is returned for positive arguments outside this + * range. For arguments less than -34.034 the cosecant + * reflection formula is applied; lograrithms are employed + * to avoid unnecessary overflow. + * + * The reciprocal gamma function has no singularities, + * but overflow and underflow may occur for large arguments. + * These conditions return either MAXNUM or 1/MAXNUM with + * appropriate sign. + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -30,+30 4000 1.2e-16 1.8e-17 + * IEEE -30,+30 30000 1.1e-15 2.0e-16 + * For arguments less than -34.034 the peak error is on the + * order of 5e-15 (DEC), excepting overflow or underflow. + */ + +/* +Cephes Math Library Release 2.8: June, 2000 +Copyright 1985, 1987, 2000 by Stephen L. Moshier +*/ + +#include <math.h> + +/* Chebyshev coefficients for reciprocal gamma function + * in interval 0 to 1. Function is 1/(x gamma(x)) - 1 + */ + +#ifdef UNK +static double R[] = { + 3.13173458231230000000E-17, +-6.70718606477908000000E-16, + 2.20039078172259550000E-15, + 2.47691630348254132600E-13, +-6.60074100411295197440E-12, + 5.13850186324226978840E-11, + 1.08965386454418662084E-9, +-3.33964630686836942556E-8, + 2.68975996440595483619E-7, + 2.96001177518801696639E-6, +-8.04814124978471142852E-5, + 4.16609138709688864714E-4, + 5.06579864028608725080E-3, +-6.41925436109158228810E-2, +-4.98558728684003594785E-3, + 1.27546015610523951063E-1 +}; +#endif + +#ifdef DEC +static unsigned short R[] = { +0022420,0066376,0176751,0071636, +0123501,0051114,0042104,0131153, +0024036,0107013,0126504,0033361, +0025613,0070040,0035174,0162316, +0126750,0037060,0077775,0122202, +0027541,0177143,0037675,0105150, +0030625,0141311,0075005,0115436, +0132017,0067714,0125033,0014721, +0032620,0063707,0105256,0152643, +0033506,0122235,0072757,0170053, +0134650,0144041,0015617,0016143, +0035332,0066125,0000776,0006215, +0036245,0177377,0137173,0131432, +0137203,0073541,0055645,0141150, +0136243,0057043,0026226,0017362, +0037402,0115554,0033441,0012310 +}; +#endif + +#ifdef IBMPC +static unsigned short R[] = { +0x2e74,0xdfbd,0x0d9f,0x3c82, +0x964d,0x8888,0x2a49,0xbcc8, +0x86de,0x75a8,0xd1c1,0x3ce3, +0x9c9a,0x074f,0x6e04,0x3d51, +0xb490,0x0fff,0x07c6,0xbd9d, +0xb14d,0x67f7,0x3fcc,0x3dcc, +0xb364,0x2f40,0xb859,0x3e12, +0x633a,0x9543,0xedf9,0xbe61, +0xdab4,0xf155,0x0cf8,0x3e92, +0xfe05,0xaebd,0xd493,0x3ec8, +0xe38c,0x2371,0x1904,0xbf15, +0xc192,0xa03f,0x4d8a,0x3f3b, +0x7663,0xf7cf,0xbfdf,0x3f74, +0xb84d,0x2b74,0x6eec,0xbfb0, +0xc3de,0x6592,0x6bc4,0xbf74, +0x2299,0x86e4,0x536d,0x3fc0 +}; +#endif + +#ifdef MIEEE +static unsigned short R[] = { +0x3c82,0x0d9f,0xdfbd,0x2e74, +0xbcc8,0x2a49,0x8888,0x964d, +0x3ce3,0xd1c1,0x75a8,0x86de, +0x3d51,0x6e04,0x074f,0x9c9a, +0xbd9d,0x07c6,0x0fff,0xb490, +0x3dcc,0x3fcc,0x67f7,0xb14d, +0x3e12,0xb859,0x2f40,0xb364, +0xbe61,0xedf9,0x9543,0x633a, +0x3e92,0x0cf8,0xf155,0xdab4, +0x3ec8,0xd493,0xaebd,0xfe05, +0xbf15,0x1904,0x2371,0xe38c, +0x3f3b,0x4d8a,0xa03f,0xc192, +0x3f74,0xbfdf,0xf7cf,0x7663, +0xbfb0,0x6eec,0x2b74,0xb84d, +0xbf74,0x6bc4,0x6592,0xc3de, +0x3fc0,0x536d,0x86e4,0x2299 +}; +#endif + +static char name[] = "rgamma"; + +#ifdef ANSIPROT +extern double chbevl ( double, void *, int ); +extern double exp ( double ); +extern double log ( double ); +extern double sin ( double ); +extern double lgam ( double ); +#else +double chbevl(), exp(), log(), sin(), lgam(); +#endif +extern double PI, MAXLOG, MAXNUM; + + +double rgamma(x) +double x; +{ +double w, y, z; +int sign; + +if( x > 34.84425627277176174) + { + mtherr( name, UNDERFLOW ); + return(1.0/MAXNUM); + } +if( x < -34.034 ) + { + w = -x; + z = sin( PI*w ); + if( z == 0.0 ) + return(0.0); + if( z < 0.0 ) + { + sign = 1; + z = -z; + } + else + sign = -1; + + y = log( w * z ) - log(PI) + lgam(w); + if( y < -MAXLOG ) + { + mtherr( name, UNDERFLOW ); + return( sign * 1.0 / MAXNUM ); + } + if( y > MAXLOG ) + { + mtherr( name, OVERFLOW ); + return( sign * MAXNUM ); + } + return( sign * exp(y)); + } +z = 1.0; +w = x; + +while( w > 1.0 ) /* Downward recurrence */ + { + w -= 1.0; + z *= w; + } +while( w < 0.0 ) /* Upward recurrence */ + { + z /= w; + w += 1.0; + } +if( w == 0.0 ) /* Nonpositive integer */ + return(0.0); +if( w == 1.0 ) /* Other integer */ + return( 1.0/z ); + +y = w * ( 1.0 + chbevl( 4.0*w-2.0, R, 16 ) ) / z; +return(y); +} |