/*							cgamma
 *
 *	Complex gamma function
 *
 *
 *
 * SYNOPSIS:
 *
 * #include <complex.h>
 * double complex x, y, cgamma();
 *
 * y = cgamma( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns complex-valued gamma function of the complex argument.
 * This variable is also filled in by the logarithmic gamma
 * function clgam().
 *
 * Arguments |x| < 18 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      -20,20      500000      2.0e-14     2.7e-15
 *    IEEE     -100,100     100000      1.4e-13     1.5e-14
 *
 * Error for arguments outside the test range will be larger
 * owing to error amplification by the exponential function.
 *
 */
/*							clgam
 *
 *	Natural logarithm of complex gamma function
 *
 *
 *
 * SYNOPSIS:
 *
 * #include <complex.h>
 * double complex x, y, clgam();
 *
 * y = clgam( 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 14 are increased by recurrence.
 * The cosecant reflection formula is employed for arguments
 * having real part less than -14.
 *
 * Arguments greater than MAXLGM return MAXNUM and an error
 * message.  MAXLGM = 2.556348e305 for IEEE arithmetic.
 *
 *
 *
 * ACCURACY:
 *
 *
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -20,20      500000      1.4e-14     4.5e-16
 *    IEEE     -100,100     100000                  1.6e-16
 * 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 double LOGPI = 1.14472988584940017414;

/* Stirling's formula for the gamma function */
#define NSTIR 7
static double STIR[ NSTIR ] = {
#if 0
 7.20489541602001055909E-5,
 8.39498720672087279993E-4,
-5.17179090826059219337E-5,
#endif
-5.92166437353693882865E-4,
 6.97281375836585777429E-5,
 7.84039221720066627474E-4,
-2.29472093621399176955E-4,
-2.68132716049382716049E-3,
 3.47222222222222222222E-3,
 8.33333333333333333333E-2
};
#define MAXSTIR 143.01608
static double SQTPI = 2.50662827463100050242E0;

extern double MAXLOG, MAXNUM, PI;
#ifdef ANSIPROT
extern double sinh (double);
extern double cosh (double);
extern double sin (double);
extern double cos (double);
extern double fabs (double);
extern double cabs (double complex);
#else
double log(), sin(), polevl(), p1evl(), floor(), fabs();
double sinh(), cosh(), cos();
double complex cpow(), cexp(), cabs();
#endif

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

/* static double complex cstirf(x) */
double complex cstirf(x)
double complex x;
{
double complex y, w;
int i;

w = 1.0/x;

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

w = 1.0 + w * y;
#if 1
y = cpow( x, x - 0.5 ) * cexp(-x);
#else
y = (x - 0.5) * clog(x) - x;
y = cexp(y);
#endif
y = SQTPI * y * w;
return( y );
}



double complex _RTL_FUNC cgamma(x)
double complex x;
{
double p, q;
double complex c, u;
int k;

if (fabs(creal(x)) > 18.0)
	{
	if( creal(x) < 0.0 )
		{
		q = creal(x);
		p = floor(q);
		if((p == q) && (cimag(x) == 0.0))
			{
			mtherr( "cgamma", OVERFLOW );
			return( MAXNUM + I * MAXNUM);
			}
		/*		c = csin( PI * x ); */
		/* Compute sin(pi x)  */
		k = q - 2.0 * floor (0.5 * q);
		q = PI * (q - p);
		p = PI * cimag(x);
		c = sin(q) * cosh(p) + cos(q) * sinh(p) * I;
		if (k & 1)
		  c = -c;
		/* Reflection formula.  */
		c = PI/(c * cgamma(1.0 - x) );
		}
	else
		{
		c = cstirf(x);
		}
	return( c );
	}

c = 1.0;
p = 0.0;
u = x;
while( creal(u) < 18.0 )
	{
	if ((fabs(creal(u)) < 1.0e-9) && (fabs(cimag(u)) < 1.0e-9))
		goto small;
	c *= u;
	p += 1.0;
	u = x + p;
	}
u = cstirf(u);
return( u / c );


small:
if((creal(x) == 0.0) && (cimag(x) == 0.0))
	{
	mtherr( "cgamma", SING );
	return( MAXNUM + MAXNUM * I );
	}
else
	return( 1.0/(((1.0 + 0.5772156649015329 * u) * u)*c) );
}



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



/* Logarithm of gamma function */

double complex _RTL_FUNC clgam(x)
double complex x;
{
double p, q, a;
double complex c, w, u, v;
int i, cj;

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

/* Reflection formula -z gamma(-z) gamma(z) = pi / sin(pi z) */
if((creal(x) < -14.0) || (cimag(x) < -14.0))
	{
	q = creal(x);
	p = floor(q);
	if( p == q )
		goto loverf;
	if (fabs(cimag(x)) > 36.7)
	  {
	    /* 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 = PI * cimag(x) - 0.6931471805599453094
	      + I * PI * (0.5 - q);
	    c = LOGPI - c - clgam(1.0 - x);
	  }
	else
	  {
	    /* Reduce sine arg mod pi.  */
	    u = csin( PI * (x - p) );
	    if( u == 0.0 )
	      goto loverf;
	    w = clgam(1.0 - x);
	    c = LOGPI - clog( u ) - w;
	    /* Adjust for reduced sine arg.  */
	    cimag(c) += PI * p;
	  }
	goto ldone;
	}
w = 0.0;
if(creal(x) < 14.0 )
	{
	  /* To satisfy Im {clgam(z)} = arg cgamma(z), accumulate
	     arg u during the recurrence.  */
	  a = 0.0;
	  w = 1.0;
	  p = 0.0;
	  u = x;
	  while( creal(u) < 14.0 )
		{
		if( u == 0.0 )
			goto loverf;
		w *= u;
		a += carg(u);
		p += 1.0;
		u = x + p;
		}
	x = u;
	w = -log(cabs(w)) - I * a;
	}

if( creal(x) > MAXLGM )
	{
loverf:
	mtherr( "clgam", OVERFLOW );
	c = MAXNUM + MAXNUM * I;
	goto ldone;
	}

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

if( cabs(x) > 1.0e8 )
  goto ldone;

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

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