blob: 771b55c779125aacb097e0f0f8b702c713153c7d [file] [log] [blame]
/**
* This file has no copyright assigned and is placed in the Public Domain.
* This file is part of the mingw-w64 runtime package.
* No warranty is given; refer to the file DISCLAIMER.PD 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 __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));
}