1 files changed, 130 insertions, 0 deletions
diff --git a/libm/float/rgammaf.c b/libm/float/rgammaf.c
new file mode 100644
index 000000000..5afa25e91
--- /dev/null
+++ b/libm/float/rgammaf.c
@@ -0,0 +1,130 @@
+/*						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);
+}