lib/libm/rgammaf.c

   1 /*                                              rgammaf.c
   2  *
   3  *      Reciprocal gamma function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, rgammaf();
  10  *
  11  * y = rgammaf( 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/MAXNUMF is returned for positive arguments outside this
  23  * range.
  24  *
  25  * The reciprocal gamma function has no singularities,
  26  * but overflow and underflow may occur for large arguments.
  27  * These conditions return either MAXNUMF or 1/MAXNUMF with
  28  * appropriate sign.
  29  *
  30  * ACCURACY:
  31  *
  32  *                      Relative error:
  33  * arithmetic   domain     # trials      peak         rms
  34  *    IEEE     -34,+34      100000      8.9e-7      1.1e-7
  35  */
  36 \f
  37 /*
  38 Cephes Math Library Release 2.2:  June, 1992
  39 Copyright 1985, 1987, 1992 by Stephen L. Moshier
  40 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  41 */
  42
  43 #include "mconf.h"
  44
  45 /* Chebyshev coefficients for reciprocal gamma function
  46  * in interval 0 to 1.  Function is 1/(x gamma(x)) - 1
  47  */
  48
  49 static float R[] = {
  50  1.08965386454418662084E-9,
  51 -3.33964630686836942556E-8,
  52  2.68975996440595483619E-7,
  53  2.96001177518801696639E-6,
  54 -8.04814124978471142852E-5,
  55  4.16609138709688864714E-4,
  56  5.06579864028608725080E-3,
  57 -6.41925436109158228810E-2,
  58 -4.98558728684003594785E-3,
  59  1.27546015610523951063E-1
  60 };
  61
  62
  63 static char name[] = "rgammaf";
  64
  65 extern float PIF, MAXLOGF, MAXNUMF;
  66
  67
  68
  69 #ifdef ANSIC
  70 float chbevlf(float, float *, int);
  71 float expf(float), logf(float), sinf(float), lgamf(float);
  72
  73 float rgammaf(float xx)
  74 #else
  75 float chbevlf(), expf(), logf(), sinf(), lgamf();
  76
  77 float rgammaf(xx)
  78 double xx;
  79 #endif
  80 {
  81 float x, w, y, z;
  82 int sign;
  83
  84 x = xx;
  85 if( x > 34.84425627277176174)
  86         {
  87         mtherr( name, UNDERFLOW );
  88         return(1.0/MAXNUMF);
  89         }
  90 if( x < -34.034 )
  91         {
  92         w = -x;
  93         z = sinf( PIF*w );
  94         if( z == 0.0 )
  95                 return(0.0);
  96         if( z < 0.0 )
  97                 {
  98                 sign = 1;
  99                 z = -z;
 100                 }
 101         else
 102                 sign = -1;
 103
 104         y = logf( w * z / PIF ) + lgamf(w);
 105         if( y < -MAXLOGF )
 106                 {
 107                 mtherr( name, UNDERFLOW );
 108                 return( sign * 1.0 / MAXNUMF );
 109                 }
 110         if( y > MAXLOGF )
 111                 {
 112                 mtherr( name, OVERFLOW );
 113                 return( sign * MAXNUMF );
 114                 }
 115         return( sign * expf(y));
 116         }
 117 z = 1.0;
 118 w = x;
 119
 120 while( w > 1.0 )        /* Downward recurrence */
 121         {
 122         w -= 1.0;
 123         z *= w;
 124         }
 125 while( w < 0.0 )        /* Upward recurrence */
 126         {
 127         z /= w;
 128         w += 1.0;
 129         }
 130 if( w == 0.0 )          /* Nonpositive integer */
 131         return(0.0);
 132 if( w == 1.0 )          /* Other integer */
 133         return( 1.0/z );
 134
 135 y = w * ( 1.0 + chbevlf( 4.0*w-2.0, R, 10 ) ) / z;
 136 return(y);
 137 }