/**
 * 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"

#ifdef UNK
static const uD P[7] = {
  { { 1.60119522476751861407E-4 } },
  { { 1.19135147006586384913E-3 } },
  { { 1.04213797561761569935E-2 } },
  { { 4.76367800457137231464E-2 } },
  { { 2.07448227648435975150E-1 } },
  { { 4.94214826801497100753E-1 } },
  { { 9.99999999999999996796E-1 } }
};
static const uD Q[8] = {
  { { -2.31581873324120129819E-5 } },
  { { 5.39605580493303397842E-4 } },
  { { -4.45641913851797240494E-3 } },
  { { 1.18139785222060435552E-2 } },
  { { 3.58236398605498653373E-2 } },
  { { -2.34591795718243348568E-1 } },
  { { 7.14304917030273074085E-2 } },
  { { 1.00000000000000000320E0 } }
};
#define MAXGAM 171.624376956302725
static const double LOGPI = 1.14472988584940017414;
#endif

#ifdef IBMPC
static const uD P[7] = {
  { { 0x2153,0x3998,0xfcb8,0x3f24 } },
  { { 0xbfab,0xe686,0x84e3,0x3f53 } },
  { { 0x14b0,0xe9db,0x57cd,0x3f85 } },
  { { 0x23d3,0x18c4,0x63d9,0x3fa8 } },
  { { 0x7d31,0xdcae,0x8da9,0x3fca } },
  { { 0xe312,0x3993,0xa137,0x3fdf } },
  { { 0x0000,0x0000,0x0000,0x3ff0 } }
};
static const uD Q[8] = {
  { { 0xd3af,0x8400,0x487a,0xbef8 } },
  { { 0x2573,0x2915,0xae8a,0x3f41 } },
  { { 0xb44a,0xe750,0x40e4,0xbf72 } },
  { { 0xb117,0x5b1b,0x31ed,0x3f88 } },
  { { 0xde67,0xe33f,0x5779,0x3fa2 } },
  { { 0x87c2,0x9d42,0x071a,0xbfce } },
  { { 0x3c51,0xc9cd,0x4944,0x3fb2 } },
  { { 0x0000,0x0000,0x0000,0x3ff0 } }
};
#define MAXGAM 171.624376956302725
#endif 

#ifdef MIEEE
static const uD P[7] = {
  { { 0x3f24,0xfcb8,0x3998,0x2153 } },
  { { 0x3f53,0x84e3,0xe686,0xbfab } },
  { { 0x3f85,0x57cd,0xe9db,0x14b0 } },
  { { 0x3fa8,0x63d9,0x18c4,0x23d3 } },
  { { 0x3fca,0x8da9,0xdcae,0x7d31 } },
  { { 0x3fdf,0xa137,0x3993,0xe312 } },
  { { 0x3ff0,0x0000,0x0000,0x0000 } }
};
static const unsigned short Q[8] = {
  { { 0xbef8,0x487a,0x8400,0xd3af } },
  { { 0x3f41,0xae8a,0x2915,0x2573 } },
  { { 0xbf72,0x40e4,0xe750,0xb44a } },
  { { 0x3f88,0x31ed,0x5b1b,0xb117 } },
  { { 0x3fa2,0x5779,0xe33f,0xde67 } },
  { { 0xbfce,0x071a,0x9d42,0x87c2 } },
  { { 0x3fb2,0x4944,0xc9cd,0x3c51 } },
  { { 0x3ff0,0x0000,0x0000,0x0000 } }
};
#define MAXGAM 171.624376956302725
#endif 

/* Stirling's formula for the gamma function */
#if UNK
static const uD STIR[5] = {
  { { 7.87311395793093628397E-4 } },
  { { -2.29549961613378126380E-4 } },
  { { -2.68132617805781232825E-3 } },
  { { 3.47222221605458667310E-3 } },
  { { 8.33333333333482257126E-2 } }
};
#define MAXSTIR 143.01608
static const double SQTPI = 2.50662827463100050242E0;
#endif
#if IBMPC
static const uD STIR[5] = {
  { { 0x7293,0x592d,0xcc72,0x3f49 } },
  { { 0x1d7c,0x27e6,0x166b,0xbf2e } },
  { { 0x4fd7,0x07d4,0xf726,0xbf65 } },
  { { 0xc5fd,0x1b98,0x71c7,0x3f6c } },
  { { 0x5986,0x5555,0x5555,0x3fb5 } }
};
#define MAXSTIR 143.01608

static const union
{
  unsigned short s[4];
  double d;
} sqt = {{0x2706,0x1ff6,0x0d93,0x4004}};
#define SQTPI (sqt.d)
#endif
#if MIEEE
static const uD STIR[5] = {
  { { 0x3f49,0xcc72,0x592d,0x7293 } },
  { { 0xbf2e,0x166b,0x27e6,0x1d7c } },
  { { 0xbf65,0xf726,0x07d4,0x4fd7 } },
  { { 0x3f6c,0x71c7,0x1b98,0xc5fd } },
  { { 0x3fb5,0x5555,0x5555,0x5986 } }
};
#define MAXSTIR 143.01608
static const uD SQT = {
  { { 0x4004,0x0d93,0x1ff6,0x2706 } }
};
#define SQTPI SQT.d
#endif

static double stirf ( double );

/* Gamma function computed by Stirling's formula.
 * The polynomial STIR is valid for 33 <= x <= 172.
 */
static double stirf(double x)
{
  double y, w, v;

  w = 1.0/x;
  w = 1.0 + w * polevl( w, STIR, 4 );
  y = exp(x);
  if( x > MAXSTIR )
  { /* Avoid overflow in pow() */
    v = pow( x, 0.5 * x - 0.25 );
    y = v * (v / y);
  }
  else
  {
    y = pow( x, x - 0.5 ) / y;
  }
  y = SQTPI * y * w;
  return( y );
}


double __tgamma_r(double x, int *sgngam)
{
  double p, q, z;
  int i;

  *sgngam = 1;
#ifdef NANS
  if( isnan(x) )
    return(x);
#endif
#ifdef INFINITIES
#ifdef NANS
  if( x == INFINITY )
    return(x);
  if( x == -INFINITY )
    return(NAN);
#else
  if( !isfinite(x) )
    return(x);
#endif
#endif
  q = fabs(x);

  if( q > 33.0 )
  {
    if( x < 0.0 )
    {
      p = floor(q);
      if( p == q )
      {
gsing:
	_SET_ERRNO(EDOM);
	mtherr( "tgamma", SING );
#ifdef INFINITIES
	return (INFINITY);
#else
	return (MAXNUM);
#endif
      }
      i = p;
      if( (i & 1) == 0 )
	*sgngam = -1;
      z = q - p;
      if( z > 0.5 )
      {
	p += 1.0;
	z = q - p;
      }
      z = q * sin( PI * z );
      if( z == 0.0 )
      {
	_SET_ERRNO(ERANGE);
	mtherr( "tgamma", OVERFLOW );
#ifdef INFINITIES
	return( *sgngam * INFINITY);
#else
	return( *sgngam * MAXNUM);
#endif
      }
      z = fabs(z);
      z = PI/(z * stirf(q) );
    }
    else
    {
      z = stirf(x);
    }
    return( *sgngam * z );
  }

  z = 1.0;
  while( x >= 3.0 )
  {
    x -= 1.0;
    z *= x;
  }

  while( x < 0.0 )
  {
    if( x > -1.E-9 )
      goto Small;
    z /= x;
    x += 1.0;
  }

  while( x < 2.0 )
  {
    if( x < 1.e-9 )
      goto Small;
    z /= x;
    x += 1.0;
  }

  if( x == 2.0 )
    return(z);

  x -= 2.0;
  p = polevl( x, P, 6 );
  q = polevl( x, Q, 7 );
  return( z * p / q );

Small:
  if( x == 0.0 )
  {
    goto gsing;
  }
  else
    return( z/((1.0 + 0.5772156649015329 * x) * x) );
}

/* This is the C99 version */

double tgamma(double x)
{
  int local_sgngam=0;
  return (__tgamma_r(x, &local_sgngam));
}