summaryrefslogtreecommitdiff
path: root/libm/float/rgammaf.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/float/rgammaf.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/float/rgammaf.c')
-rw-r--r--libm/float/rgammaf.c130
1 files changed, 0 insertions, 130 deletions
diff --git a/libm/float/rgammaf.c b/libm/float/rgammaf.c
deleted file mode 100644
index 5afa25e91..000000000
--- a/libm/float/rgammaf.c
+++ /dev/null
@@ -1,130 +0,0 @@
-/* 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 <math.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;
-
-
-
-float chbevlf(float, float *, int);
-float expf(float), logf(float), sinf(float), lgamf(float);
-
-float rgammaf(float xx)
-{
-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);
-}