Quelle  cgamma.c

/* 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 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 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 );
}