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, 0 insertions, 423 deletions
diff --git a/libm/float/gammaf.c b/libm/float/gammaf.c
deleted file mode 100644
index e8c4694c4..000000000
--- a/libm/float/gammaf.c
+++ /dev/null
@@ -1,423 +0,0 @@
-/* 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 );
-}
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-26 03:37:10 +0000