/* lgamf()
*
* Natural logarithm of gamma function
*
*
*
* SYNOPSIS:
*
* float x, y, __lgammaf_r();
* int* sgngamf;
* y = __lgammaf_r( x, sgngamf );
*
* float x, y, lgammaf();
* y = lgammaf( x);
*
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of the absolute
* value of the gamma function of the argument. In the reentrant
* version the sign (+1 or -1) of the gamma function is returned in
* variable referenced by 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
*
*/
/*
Cephes Math Library Release 2.7: July, 1998
Copyright 1984, 1987, 1989, 1992, 1998 by Stephen L. Moshier
*/
/*
26-11-2002 Modified for mingw.
Danny Smith <[email protected]>
*/
/* log gamma(x+2), -.5 < x < .5 */
static const 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 const 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 const float LS2PI = 0.91893853320467274178;
#define MAXLGM 2.035093e36
static const float PIINV = 0.318309886183790671538;
#ifndef __MINGW32__
#include "mconf.h"
float floorf(float);
float polevlf( float, float *, int );
float p1evlf( float, float *, int );
#else
#include "cephes_mconf.h"
#endif
/* Reentrant version */
/* Logarithm of gamma function */
float __lgammaf_r( float x, int* sgngamf )
{
float p, q, w, z;
float nx, tx;
int i, direction;
*sgngamf = 1;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
#ifdef INFINITIES
if( !isfinite(x) )
return(x);
#endif
if( x < 0.0 )
{
q = -x;
w = __lgammaf_r(q, sgngamf); /* note this modifies sgngam! */
p = floorf(q);
if( p == q )
{
lgsing:
_SET_ERRNO(EDOM);
mtherr( "lgamf", SING );
#ifdef INFINITIES
return (INFINITYF);
#else
return( *sgngamf * MAXNUMF );
#endif
}
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 lgsing;
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 lgsing;
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 )
{
_SET_ERRNO(ERANGE);
mtherr( "lgamf", OVERFLOW );
#ifdef INFINITIES
return( *sgngamf * INFINITYF );
#else
return( *sgngamf * MAXNUMF );
#endif
}
/* 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 );
}
/* This is the C99 version */
float lgammaf(float x)
{
int local_sgngamf=0;
return (__lgammaf_r(x, &local_sgngamf));
}