mingw-w64-crt/math/lgammaf.c

/**
 * This file has no copyright assigned and is placed in the Public Domain.
 * This file is part of the w64 mingw-runtime package.
 * No warranty is given; refer to the file DISCLAIMER within this package.
 */
/* 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;

#include "cephes_mconf.h"

/* 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));
}