lib/libm/rgamma.c

   1 /*                                              rgamma.c
   2  *
   3  *      Reciprocal gamma function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, rgamma();
  10  *
  11  * y = rgamma( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns one divided by the gamma function of the argument.
  18  *
  19  * The function is approximated by a Chebyshev expansion in
  20  * the interval [0,1].  Range reduction is by recurrence
  21  * for arguments between -34.034 and +34.84425627277176174.
  22  * 1/MAXNUM is returned for positive arguments outside this
  23  * range.  For arguments less than -34.034 the cosecant
  24  * reflection formula is applied; lograrithms are employed
  25  * to avoid unnecessary overflow.
  26  *
  27  * The reciprocal gamma function has no singularities,
  28  * but overflow and underflow may occur for large arguments.
  29  * These conditions return either MAXNUM or 1/MAXNUM with
  30  * appropriate sign.
  31  *
  32  * ACCURACY:
  33  *
  34  *                      Relative error:
  35  * arithmetic   domain     # trials      peak         rms
  36  *    DEC      -30,+30       4000       1.2e-16     1.8e-17
  37  *    IEEE     -30,+30      30000       1.1e-15     2.0e-16
  38  * For arguments less than -34.034 the peak error is on the
  39  * order of 5e-15 (DEC), excepting overflow or underflow.
  40  */
  41 \f
  42 /*
  43 Cephes Math Library Release 2.8:  June, 2000
  44 Copyright 1985, 1987, 2000 by Stephen L. Moshier
  45 */
  46
  47 #include "mconf.h"
  48
  49 /* Chebyshev coefficients for reciprocal gamma function
  50  * in interval 0 to 1.  Function is 1/(x gamma(x)) - 1
  51  */
  52
  53 #ifdef UNK
  54 static double R[] = {
  55  3.13173458231230000000E-17,
  56 -6.70718606477908000000E-16,
  57  2.20039078172259550000E-15,
  58  2.47691630348254132600E-13,
  59 -6.60074100411295197440E-12,
  60  5.13850186324226978840E-11,
  61  1.08965386454418662084E-9,
  62 -3.33964630686836942556E-8,
  63  2.68975996440595483619E-7,
  64  2.96001177518801696639E-6,
  65 -8.04814124978471142852E-5,
  66  4.16609138709688864714E-4,
  67  5.06579864028608725080E-3,
  68 -6.41925436109158228810E-2,
  69 -4.98558728684003594785E-3,
  70  1.27546015610523951063E-1
  71 };
  72 #endif
  73
  74 #ifdef DEC
  75 static unsigned short R[] = {
  76 0022420,0066376,0176751,0071636,
  77 0123501,0051114,0042104,0131153,
  78 0024036,0107013,0126504,0033361,
  79 0025613,0070040,0035174,0162316,
  80 0126750,0037060,0077775,0122202,
  81 0027541,0177143,0037675,0105150,
  82 0030625,0141311,0075005,0115436,
  83 0132017,0067714,0125033,0014721,
  84 0032620,0063707,0105256,0152643,
  85 0033506,0122235,0072757,0170053,
  86 0134650,0144041,0015617,0016143,
  87 0035332,0066125,0000776,0006215,
  88 0036245,0177377,0137173,0131432,
  89 0137203,0073541,0055645,0141150,
  90 0136243,0057043,0026226,0017362,
  91 0037402,0115554,0033441,0012310
  92 };
  93 #endif
  94
  95 #ifdef IBMPC
  96 static unsigned short R[] = {
  97 0x2e74,0xdfbd,0x0d9f,0x3c82,
  98 0x964d,0x8888,0x2a49,0xbcc8,
  99 0x86de,0x75a8,0xd1c1,0x3ce3,
 100 0x9c9a,0x074f,0x6e04,0x3d51,
 101 0xb490,0x0fff,0x07c6,0xbd9d,
 102 0xb14d,0x67f7,0x3fcc,0x3dcc,
 103 0xb364,0x2f40,0xb859,0x3e12,
 104 0x633a,0x9543,0xedf9,0xbe61,
 105 0xdab4,0xf155,0x0cf8,0x3e92,
 106 0xfe05,0xaebd,0xd493,0x3ec8,
 107 0xe38c,0x2371,0x1904,0xbf15,
 108 0xc192,0xa03f,0x4d8a,0x3f3b,
 109 0x7663,0xf7cf,0xbfdf,0x3f74,
 110 0xb84d,0x2b74,0x6eec,0xbfb0,
 111 0xc3de,0x6592,0x6bc4,0xbf74,
 112 0x2299,0x86e4,0x536d,0x3fc0
 113 };
 114 #endif
 115
 116 #ifdef MIEEE
 117 static unsigned short R[] = {
 118 0x3c82,0x0d9f,0xdfbd,0x2e74,
 119 0xbcc8,0x2a49,0x8888,0x964d,
 120 0x3ce3,0xd1c1,0x75a8,0x86de,
 121 0x3d51,0x6e04,0x074f,0x9c9a,
 122 0xbd9d,0x07c6,0x0fff,0xb490,
 123 0x3dcc,0x3fcc,0x67f7,0xb14d,
 124 0x3e12,0xb859,0x2f40,0xb364,
 125 0xbe61,0xedf9,0x9543,0x633a,
 126 0x3e92,0x0cf8,0xf155,0xdab4,
 127 0x3ec8,0xd493,0xaebd,0xfe05,
 128 0xbf15,0x1904,0x2371,0xe38c,
 129 0x3f3b,0x4d8a,0xa03f,0xc192,
 130 0x3f74,0xbfdf,0xf7cf,0x7663,
 131 0xbfb0,0x6eec,0x2b74,0xb84d,
 132 0xbf74,0x6bc4,0x6592,0xc3de,
 133 0x3fc0,0x536d,0x86e4,0x2299
 134 };
 135 #endif
 136
 137 static char name[] = "rgamma";
 138
 139 #ifdef ANSIPROT
 140 extern double chbevl ( double, void *, int );
 141 extern double exp ( double );
 142 extern double log ( double );
 143 extern double sin ( double );
 144 extern double lgam ( double );
 145 #else
 146 double chbevl(), exp(), log(), sin(), lgam();
 147 #endif
 148 extern double PI, MAXLOG, MAXNUM;
 149
 150
 151 double rgamma(x)
 152 double x;
 153 {
 154 double w, y, z;
 155 int sign;
 156
 157 if( x > 34.84425627277176174)
 158         {
 159         mtherr( name, UNDERFLOW );
 160         return(1.0/MAXNUM);
 161         }
 162 if( x < -34.034 )
 163         {
 164         w = -x;
 165         z = sin( PI*w );
 166         if( z == 0.0 )
 167                 return(0.0);
 168         if( z < 0.0 )
 169                 {
 170                 sign = 1;
 171                 z = -z;
 172                 }
 173         else
 174                 sign = -1;
 175
 176         y = log( w * z ) - log(PI) + lgam(w);
 177         if( y < -MAXLOG )
 178                 {
 179                 mtherr( name, UNDERFLOW );
 180                 return( sign * 1.0 / MAXNUM );
 181                 }
 182         if( y > MAXLOG )
 183                 {
 184                 mtherr( name, OVERFLOW );
 185                 return( sign * MAXNUM );
 186                 }
 187         return( sign * exp(y));
 188         }
 189 z = 1.0;
 190 w = x;
 191
 192 while( w > 1.0 )        /* Downward recurrence */
 193         {
 194         w -= 1.0;
 195         z *= w;
 196         }
 197 while( w < 0.0 )        /* Upward recurrence */
 198         {
 199         z /= w;
 200         w += 1.0;
 201         }
 202 if( w == 0.0 )          /* Nonpositive integer */
 203         return(0.0);
 204 if( w == 1.0 )          /* Other integer */
 205         return( 1.0/z );
 206
 207 y = w * ( 1.0 + chbevl( 4.0*w-2.0, R, 16 ) ) / z;
 208 return(y);
 209 }