summaryrefslogtreecommitdiff
path: root/libm/float/gammaf.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/float/gammaf.c')
-rw-r--r--libm/float/gammaf.c423
1 files changed, 423 insertions, 0 deletions
diff --git a/libm/float/gammaf.c b/libm/float/gammaf.c
new file mode 100644
index 000000000..e8c4694c4
--- /dev/null
+++ b/libm/float/gammaf.c
@@ -0,0 +1,423 @@
+/* gammaf.c
+ *
+ * Gamma function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, gammaf();
+ * extern int sgngamf;
+ *
+ * y = gammaf( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns gamma function of the argument. The result is
+ * correctly signed, and the sign (+1 or -1) is also
+ * returned in a global (extern) variable named sgngamf.
+ * This same variable is also filled in by the logarithmic
+ * gamma function lgam().
+ *
+ * Arguments between 0 and 10 are reduced by recurrence and the
+ * function is approximated by a polynomial function covering
+ * the interval (2,3). Large arguments are handled by Stirling's
+ * formula. Negative arguments are made positive using
+ * a reflection formula.
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0,-33 100,000 5.7e-7 1.0e-7
+ * IEEE -33,0 100,000 6.1e-7 1.2e-7
+ *
+ *
+ */
+/* lgamf()
+ *
+ * Natural logarithm of gamma function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, lgamf();
+ * extern int sgngamf;
+ *
+ * y = lgamf( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base e (2.718...) logarithm of the absolute
+ * value of the gamma function of the argument.
+ * The sign (+1 or -1) of the gamma function is returned in a
+ * global (extern) variable named sgngamf.
+ *
+ * For arguments greater than 6.5, the logarithm of the gamma
+ * function is approximated by the logarithmic version of
+ * Stirling's formula. Arguments between 0 and +6.5 are reduced by
+ * by recurrence to the interval [.75,1.25] or [1.5,2.5] of a rational
+ * approximation. The cosecant reflection formula is employed for
+ * arguments less than zero.
+ *
+ * Arguments greater than MAXLGM = 2.035093e36 return MAXNUM and an
+ * error message.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *
+ *
+ * arithmetic domain # trials peak rms
+ * IEEE -100,+100 500,000 7.4e-7 6.8e-8
+ * The error criterion was relative when the function magnitude
+ * was greater than one but absolute when it was less than one.
+ * The routine has low relative error for positive arguments.
+ *
+ * The following test used the relative error criterion.
+ * IEEE -2, +3 100000 4.0e-7 5.6e-8
+ *
+ */
+
+/* gamma.c */
+/* gamma function */
+
+/*
+Cephes Math Library Release 2.7: July, 1998
+Copyright 1984, 1987, 1989, 1992, 1998 by Stephen L. Moshier
+*/
+
+
+#include <math.h>
+
+/* define MAXGAM 34.84425627277176174 */
+
+/* Stirling's formula for the gamma function
+ * gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) ( 1 + 1/x P(1/x) )
+ * .028 < 1/x < .1
+ * relative error < 1.9e-11
+ */
+static float STIR[] = {
+-2.705194986674176E-003,
+ 3.473255786154910E-003,
+ 8.333331788340907E-002,
+};
+static float MAXSTIR = 26.77;
+static float SQTPIF = 2.50662827463100050242; /* sqrt( 2 pi ) */
+
+int sgngamf = 0;
+extern int sgngamf;
+extern float MAXLOGF, MAXNUMF, PIF;
+
+#ifdef ANSIC
+float expf(float);
+float logf(float);
+float powf( float, float );
+float sinf(float);
+float gammaf(float);
+float floorf(float);
+static float stirf(float);
+float polevlf( float, float *, int );
+float p1evlf( float, float *, int );
+#else
+float expf(), logf(), powf(), sinf(), floorf();
+float polevlf(), p1evlf();
+static float stirf();
+#endif
+
+/* Gamma function computed by Stirling's formula,
+ * sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
+ * The polynomial STIR is valid for 33 <= x <= 172.
+ */
+static float stirf( float xx )
+{
+float x, y, w, v;
+
+x = xx;
+w = 1.0/x;
+w = 1.0 + w * polevlf( w, STIR, 2 );
+y = expf( -x );
+if( x > MAXSTIR )
+ { /* Avoid overflow in pow() */
+ v = powf( x, 0.5 * x - 0.25 );
+ y *= v;
+ y *= v;
+ }
+else
+ {
+ y = powf( x, x - 0.5 ) * y;
+ }
+y = SQTPIF * y * w;
+return( y );
+}
+
+
+/* gamma(x+2), 0 < x < 1 */
+static float P[] = {
+ 1.536830450601906E-003,
+ 5.397581592950993E-003,
+ 4.130370201859976E-003,
+ 7.232307985516519E-002,
+ 8.203960091619193E-002,
+ 4.117857447645796E-001,
+ 4.227867745131584E-001,
+ 9.999999822945073E-001,
+};
+
+float gammaf( float xx )
+{
+float p, q, x, z, nz;
+int i, direction, negative;
+
+x = xx;
+sgngamf = 1;
+negative = 0;
+nz = 0.0;
+if( x < 0.0 )
+ {
+ negative = 1;
+ q = -x;
+ p = floorf(q);
+ if( p == q )
+ goto goverf;
+ i = p;
+ if( (i & 1) == 0 )
+ sgngamf = -1;
+ nz = q - p;
+ if( nz > 0.5 )
+ {
+ p += 1.0;
+ nz = q - p;
+ }
+ nz = q * sinf( PIF * nz );
+ if( nz == 0.0 )
+ {
+goverf:
+ mtherr( "gamma", OVERFLOW );
+ return( sgngamf * MAXNUMF);
+ }
+ if( nz < 0 )
+ nz = -nz;
+ x = q;
+ }
+if( x >= 10.0 )
+ {
+ z = stirf(x);
+ }
+if( x < 2.0 )
+ direction = 1;
+else
+ direction = 0;
+z = 1.0;
+while( x >= 3.0 )
+ {
+ x -= 1.0;
+ z *= x;
+ }
+/*
+while( x < 0.0 )
+ {
+ if( x > -1.E-4 )
+ goto small;
+ z *=x;
+ x += 1.0;
+ }
+*/
+while( x < 2.0 )
+ {
+ if( x < 1.e-4 )
+ goto small;
+ z *=x;
+ x += 1.0;
+ }
+
+if( direction )
+ z = 1.0/z;
+
+if( x == 2.0 )
+ return(z);
+
+x -= 2.0;
+p = z * polevlf( x, P, 7 );
+
+gdone:
+
+if( negative )
+ {
+ p = sgngamf * PIF/(nz * p );
+ }
+return(p);
+
+small:
+if( x == 0.0 )
+ {
+ mtherr( "gamma", SING );
+ return( MAXNUMF );
+ }
+else
+ {
+ p = z / ((1.0 + 0.5772156649015329 * x) * x);
+ goto gdone;
+ }
+}
+
+
+
+
+/* log gamma(x+2), -.5 < x < .5 */
+static float B[] = {
+ 6.055172732649237E-004,
+-1.311620815545743E-003,
+ 2.863437556468661E-003,
+-7.366775108654962E-003,
+ 2.058355474821512E-002,
+-6.735323259371034E-002,
+ 3.224669577325661E-001,
+ 4.227843421859038E-001
+};
+
+/* log gamma(x+1), -.25 < x < .25 */
+static float C[] = {
+ 1.369488127325832E-001,
+-1.590086327657347E-001,
+ 1.692415923504637E-001,
+-2.067882815621965E-001,
+ 2.705806208275915E-001,
+-4.006931650563372E-001,
+ 8.224670749082976E-001,
+-5.772156501719101E-001
+};
+
+/* log( sqrt( 2*pi ) ) */
+static float LS2PI = 0.91893853320467274178;
+#define MAXLGM 2.035093e36
+static float PIINV = 0.318309886183790671538;
+
+/* Logarithm of gamma function */
+
+
+float lgamf( float xx )
+{
+float p, q, w, z, x;
+float nx, tx;
+int i, direction;
+
+sgngamf = 1;
+
+x = xx;
+if( x < 0.0 )
+ {
+ q = -x;
+ w = lgamf(q); /* note this modifies sgngam! */
+ p = floorf(q);
+ if( p == q )
+ goto loverf;
+ i = p;
+ if( (i & 1) == 0 )
+ sgngamf = -1;
+ else
+ sgngamf = 1;
+ z = q - p;
+ if( z > 0.5 )
+ {
+ p += 1.0;
+ z = p - q;
+ }
+ z = q * sinf( PIF * z );
+ if( z == 0.0 )
+ goto loverf;
+ z = -logf( PIINV*z ) - w;
+ return( z );
+ }
+
+if( x < 6.5 )
+ {
+ direction = 0;
+ z = 1.0;
+ tx = x;
+ nx = 0.0;
+ if( x >= 1.5 )
+ {
+ while( tx > 2.5 )
+ {
+ nx -= 1.0;
+ tx = x + nx;
+ z *=tx;
+ }
+ x += nx - 2.0;
+iv1r5:
+ p = x * polevlf( x, B, 7 );
+ goto cont;
+ }
+ if( x >= 1.25 )
+ {
+ z *= x;
+ x -= 1.0; /* x + 1 - 2 */
+ direction = 1;
+ goto iv1r5;
+ }
+ if( x >= 0.75 )
+ {
+ x -= 1.0;
+ p = x * polevlf( x, C, 7 );
+ q = 0.0;
+ goto contz;
+ }
+ while( tx < 1.5 )
+ {
+ if( tx == 0.0 )
+ goto loverf;
+ z *=tx;
+ nx += 1.0;
+ tx = x + nx;
+ }
+ direction = 1;
+ x += nx - 2.0;
+ p = x * polevlf( x, B, 7 );
+
+cont:
+ if( z < 0.0 )
+ {
+ sgngamf = -1;
+ z = -z;
+ }
+ else
+ {
+ sgngamf = 1;
+ }
+ q = logf(z);
+ if( direction )
+ q = -q;
+contz:
+ return( p + q );
+ }
+
+if( x > MAXLGM )
+ {
+loverf:
+ mtherr( "lgamf", OVERFLOW );
+ return( sgngamf * MAXNUMF );
+ }
+
+/* Note, though an asymptotic formula could be used for x >= 3,
+ * there is cancellation error in the following if x < 6.5. */
+q = LS2PI - x;
+q += ( x - 0.5 ) * logf(x);
+
+if( x <= 1.0e4 )
+ {
+ z = 1.0/x;
+ p = z * z;
+ q += (( 6.789774945028216E-004 * p
+ - 2.769887652139868E-003 ) * p
+ + 8.333316229807355E-002 ) * z;
+ }
+return( q );
+}