/* cgammaf
 *
 * Complex gamma function, single precision
 *
 *
 *
 * SYNOPSIS:
 *
 * #include <complex.h>
 * float complex x, y, cgammaf();
 *
 * y = cgammaf( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns complex-valued gamma function of the complex argument.
 *
 * Arguments |Re(x)| <= 14 are increased by recurrence.
 * Large arguments are handled by Stirling's formula. Large negative
 * arguments are made positive using the reflection formula.  
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -10,10      100000      5.2e-6      5.8e-7
 *    IEEE      -22,22      100000      1.1e-5      1.2e-6
 *    IEEE      -27,27      100000      4.6e-5      1.5e-6
 *
 * Error for arguments outside the test range will be larger
 * owing to error amplification by the exponential function.
 *
 */
/* clgamf
 *
 * Natural logarithm of complex gamma function, single precision
 *
 *
 *
 * SYNOPSIS:
 *
 * #include <complex.h>
 * float complex x, y, clgamf();
 *
 * y = clgamf( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the base e (2.718...) logarithm of the complex gamma
 * function of the argument.
 *
 * The logarithm of the gamma function is approximated by the
 * logarithmic version of Stirling's asymptotic formula.
 * Arguments of real part less than 8 are increased by recurrence.
 * The cosecant reflection formula is employed for arguments
 * having real part less than -8.
 *
 * Arguments greater than MAXLGM return MAXNUM and an error
 * message.  MAXLGM = 2.035093e36.
 *
 *
 *
 * ACCURACY:
 *
 *
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -20,20     9 000 000    2.6e-6       1.1e-7
 *    IEEE     -100,100      500 000                 7.5e-8
 * The error criterion was relative when the function magnitude
 * was greater than one but absolute when it was less than one.
 */

/*
Cephes Math Library Release 2.7:  March, 1998
Copyright 1984, 1998 Stephen L. Moshier
*/

#include "complex.h"
#include "mconf.h"

#define MAXGAM 171.624376956302725
static float LOGPI = 1.14472988584940017414f;

/* Stirling's formula for the gamma function */
#define NSTIR 7
static float STIR[ NSTIR ] = {
 -5.92166437353693882865E-4,
 6.97281375836585777429E-5,
 7.84039221720066627474E-4,
-2.29472093621399176955E-4,
-2.68132716049382716049E-3,
 3.47222222222222222222E-3,
 8.33333333333333333333E-2
};
#define MAXSTIR 26.0
static float SQTPI = 2.50662827463100050242f;

extern float MAXLOGF, MAXNUMF, PIF;
#ifdef ANSIPROT
extern float sinhf ( float x );
extern float coshf ( float x );
extern float sinf ( float x );
extern float cosf ( float x );
#else
float logf(), sinf(), polevlf(), p1evlf(), floorf(), fabsf();
float cosf(), sinhf(), coshf();
float complex cpowf(), cexpf(), cabsf();
#endif

/* Gamma function computed by Stirling's formula.  */

static float complex cstirf(x)
float complex x;
{
float complex y, v, w;
int i;

w = 1.0f/x;

y = STIR[0];
for (i = 1; i < NSTIR; i++)
  {
    y = y * w + STIR[i];
  }

w = 1.0f + w * y;
if( cabsf(x) > MAXSTIR )
 { /* Avoid overflow in pow() */
 v = cpowf( x, 0.5f * x - 0.25f );
 y = v * (v * cexpf(-x));
 }
else
 {
 y = cpowf( x, x - 0.5f ) * cexpf(-x);
 }
/* y = cpowf( x, x - 0.5f ) * cexpf(-x); */

y = SQTPI * y * w;
return( y );
}



float complex cgammaf(x)
float complex x;
{
float p, q;
float complex c, u;
int k;

if( fabsf(crealf(x)) > 13.0f )
 {
 if( crealf(x) < 0.0f )
  {
  q = crealf(x);
  p = floorf(q);
  if((p == q ) && (cimag(x) == 0.0f))
   {
cgoverf:
   mtherr( "cgammaf", OVERFLOW );
   return( MAXNUMF + I * MAXNUMF);
   }
  /* Compute sin(pi x)  */
  k = q - 2.0f * floorf (0.5f * q);
  q = PIF * (q - p);
  p = PIF * cimag(x);
  c = sinf(q) * coshf(p) + cosf(q) * sinhf(p) * I;
  if (k & 1)
    c = -c;
  /* Reflection formula.
     ??? The denominator has been observed to overflow,
     producing a NaN result on a SPARCstation instead of zero.
     Example:
   cgammaf(-19.425196335160404f - 18.188121066880587f*I) */
  c = PIF/(c * cgammaf(1.0f - x));
  }
 else
  {
  c = cstirf(x);
  }
 return( c );
 }

c = 1.0f;
p = 0.0f;
u = x;
while( crealf(u) < 13.0f )
 {
 if ((crealf(u) == 0.0f) && (cimagf(u) == 0.0f))
  goto cgoverf;
 c *= u;
 p += 1.0f;
 u = x + p;
 }
u = cstirf(u);
return( u / c );
}



/* Asymptotic expansion of log gamma  */
static float A[] = {
#if 0
-1.3924322169059011164274322169059011164274E0,
 1.7964437236883057316493849001588939669435E-1,
-2.9550653594771241830065359477124183006536E-2,
 6.4102564102564102564102564102564102564103E-3,
-1.9175269175269175269175269175269175269175E-3,
 8.4175084175084175084175084175084175084175E-4,
#endif
-5.9523809523809523809523809523809523809524E-4,
 7.9365079365079365079365079365079365079365E-4,
-2.7777777777777777777777777777777777777778E-3,
 8.3333333333333333333333333333333333333333E-2
};
/* log( sqrt( 2*pi ) ) */
static float LS2PI  =  0.91893853320467274178f;
#define MAXLGM 2.035093e36f



/* Logarithm of gamma function */

float complex clgamf(x)
float complex x;
{
float complex c, w, u, v;
float p, q, a;
int i, cj;

cj = 0;
if (cimagf(x) < 0)
  {
    cj = 1;
    x = conj(x);
  }

/* Reflection formula -z gamma(-z) gamma(z) = pi / sin(pi z)  */
if((crealf(x) < -7.0f) || (cimagf(x) < -7.0f))
 {
 q = crealf(x);
 p = floorf(q);
 if( p == q )
  goto loverf;
 if (fabsf(cimag(x)) > 18.4f)
   {
     /* sin z grows exponentially with Im(z).  Find ln sin(pi z)
        from |sin z| = sqrt( sin^2 x + sinh^2 y),
               arg sin z = arctan(tanh y / tan x).  */
     c = PIF * cimagf(x) - 0.6931471805599453094f
       + I * PIF * (0.5f - q);
     c = LOGPI - c - clgamf(1.0f - x);
   }
 else
   {
     /* Reduce sine arg mod pi.  */
     u = csinf( PIF * (x - p) );
     if( u == 0.0f )
       goto loverf;
     w = clgamf(1.0f - x);
     c = LOGPI - clogf( u ) - w;
     /* Adjust for reduced sine arg.  */
     cimagf(c) += PIF * p;
   }
 goto ldone;
 }
w = 0.0f;
if( crealf(x) < 7.0f )
 {
   /* To satisfy Im {clgam(z)} = arg cgamma(z), accumulate
     arg u during the recurrence.  */
   a = 0.0f;
   w = 1.0f;
   p = 0.0f;
   u = x;
   while( crealf(u) < 7.0f )
  {
  if( u == 0.0f )
   goto loverf;
  w *= u;
  a += cargf(u);
  p += 1.0f;
  u = x + p;
  }
 x = u;
 w = -logf(cabsf(w)) - I * a;
 }

if( crealf(x) > MAXLGM )
 {
loverf:
 mtherr( "clgamf", OVERFLOW );
 c = MAXNUMF + MAXNUMF * I;
 goto ldone;
 }

c = ( x - 0.5f ) * clogf(x) - x + LS2PI + w;

if( cabsf(x) > 1.0e8f )
  goto ldone;

v = 1.0f/(x*x);
u = A[0];
for (i = 1; i < 4; i++)
  {
    u = u * v + A[i];
  }
c = c + u / x;

ldone:
if (cj)
  c = conj(c);
return( c );
}