/**
 * 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.
 */
#include "cephes_mconf.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 const float STIR[] = {
-2.705194986674176E-003,
 3.473255786154910E-003,
 8.333331788340907E-002,
};
static const float MAXSTIR = 26.77;
static const float SQTPIF = 2.50662827463100050242; /* sqrt( 2 pi ) */

static float stirf(float);

/* 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 x )
{
float  y, w, v;

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 const 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 __tgammaf_r( float x, int* sgngamf)
{
float p, q, z, nz;
int i, direction, negative;

#ifdef NANS
if( isnan(x) )
	return(x);
#endif
#ifdef INFINITIES
#ifdef NANS
if( x == INFINITYF )
	return(x);
if( x == -INFINITYF )
	return(NANF);
#else
if( !isfinite(x) )
	return(x);
#endif
#endif

*sgngamf = 1;
negative = 0;
nz = 0.0;
if( x < 0.0 )
	{
	negative = 1;
	q = -x;
	p = floorf(q);
	if( p == q )
		{
gsing:
		_SET_ERRNO(EDOM);
		mtherr( "tgammaf", SING );
#ifdef INFINITIES
		return (INFINITYF);
#else
		return (MAXNUMF);
#endif
			}
	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 )
		{
		_SET_ERRNO(ERANGE);
		mtherr( "tgamma", OVERFLOW );
#ifdef INFINITIES
		return( *sgngamf * INFINITYF);
#else
		return( *sgngamf * MAXNUMF);
#endif
		}
	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 )
	{
	goto gsing;
	}
else
	{
	p = z / ((1.0 + 0.5772156649015329 * x) * x);
	goto gdone;
	}
}

/* This is the C99 version */

float tgammaf(float x)
{
  int local_sgngamf=0;
  return (__tgammaf_r(x, &local_sgngamf));
}