text: /*						rgammaf.c
   *
   *	Reciprocal gamma function
   *
   *
   *
   * SYNOPSIS:
   *
   * float x, y, rgammaf();
   *
   * y = rgammaf( 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/MAXNUMF is returned for positive arguments outside this
   * range.
   *
   * The reciprocal gamma function has no singularities,
   * but overflow and underflow may occur for large arguments.
   * These conditions return either MAXNUMF or 1/MAXNUMF with
   * appropriate sign.
   *
   * ACCURACY:
   *
   *                      Relative error:
   * arithmetic   domain     # trials      peak         rms
   *    IEEE     -34,+34      100000      8.9e-7      1.1e-7
   */
  
  /*
  Cephes Math Library Release 2.2:  June, 1992
  Copyright 1985, 1987, 1992 by Stephen L. Moshier
  Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  */
  
  #include "mconf.h"
  
  /* Chebyshev coefficients for reciprocal gamma function
   * in interval 0 to 1.  Function is 1/(x gamma(x)) - 1
   */
  
  static float R[] = {
   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
  };
  
  
  static char name[] = "rgammaf";
  
  extern float PIF, MAXLOGF, MAXNUMF;
  
  
  
  #ifdef ANSIC
  float chbevlf(float, float *, int);
  float expf(float), logf(float), sinf(float), lgamf(float);
  
  float rgammaf(float xx)
  #else
  float chbevlf(), expf(), logf(), sinf(), lgamf();
  
  float rgammaf(xx)
  double xx;
  #endif
  {
  float x, w, y, z;
  int sign;
  
  x = xx;
  if( x > 34.84425627277176174)
  	{
  	mtherr( name, UNDERFLOW );
  	return(1.0/MAXNUMF);
  	}
  if( x < -34.034 )
  	{
  	w = -x;
  	z = sinf( PIF*w );
  	if( z == 0.0 )
  		return(0.0);
  	if( z < 0.0 )
  		{
  		sign = 1;
  		z = -z;
  		}
  	else
  		sign = -1;
  
  	y = logf( w * z / PIF ) + lgamf(w);
  	if( y < -MAXLOGF )
  		{
  		mtherr( name, UNDERFLOW );
  		return( sign * 1.0 / MAXNUMF );
  		}
  	if( y > MAXLOGF )
  		{
  		mtherr( name, OVERFLOW );
  		return( sign * MAXNUMF );
  		}
  	return( sign * expf(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 + chbevlf( 4.0*w-2.0, R, 10 ) ) / z;
  return(y);
  }
  