/* igam.c
*
* Incomplete md_gamma integral
*
*
*
* SYNOPSIS:
*
* double a, x, y, igam();
*
* y = igam( a, x );
*
* DESCRIPTION:
*
* The function is defined by
*
* x
* -
* 1 | | -t a-1
* igam(a,x) = ----- | e t dt.
* - | |
* | (a) -
* 0
*
*
* In this implementation both arguments must be positive.
* The integral is evaluated by either a power series or
* continued fraction expansion, depending on the relative
* values of a and x.
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,30 200000 3.6e-14 2.9e-15
* IEEE 0,100 300000 9.9e-14 1.5e-14
*/
/* igamc()
*
* Complemented incomplete md_gamma integral
*
*
*
* SYNOPSIS:
*
* double a, x, y, igamc();
*
* y = igamc( a, x );
*
* DESCRIPTION:
*
* The function is defined by
*
*
* igamc(a,x) = 1 - igam(a,x)
*
* inf.
* -
* 1 | | -t a-1
* = ----- | e t dt.
* - | |
* | (a) -
* x
*
*
* In this implementation both arguments must be positive.
* The integral is evaluated by either a power series or
* continued fraction expansion, depending on the relative
* values of a and x.
*
* ACCURACY:
*
* Tested at random a, x.
* a x Relative error:
* arithmetic domain domain # trials peak rms
* IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15
* IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1985, 1987, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef ANSIPROT
extern double lgam ( double );
extern double md_exp ( double );
extern double md_log ( double );
extern double md_fabs ( double );
extern double igam ( double, double );
extern double igamc ( double, double );
#else
double lgam(), md_exp(), md_log(), md_fabs(), igam(), igamc();
#endif
extern double MACHEP, MAXLOG;
static double big = 4.503599627370496e15;
static double biginv = 2.22044604925031308085e-16;
double igamc( a, x )
double a, x;
{
double ans, ax, c, yc, r, t, y, z;
double pk, pkm1, pkm2, qk, qkm1, qkm2;
if( (x <= 0) || ( a <= 0) )
return( 1.0 );
if( (x < 1.0) || (x < a) )
return( 1.0 - igam(a,x) );
ax = a * md_log(x) - x - lgam(a);
if( ax < -MAXLOG )
{
mtherr( "igamc", UNDERFLOW );
return( 0.0 );
}
ax = md_exp(ax);
/* continued fraction */
y = 1.0 - a;
z = x + y + 1.0;
c = 0.0;
pkm2 = 1.0;
qkm2 = x;
pkm1 = x + 1.0;
qkm1 = z * x;
ans = pkm1/qkm1;
do
{
c += 1.0;
y += 1.0;
z += 2.0;
yc = y * c;
pk = pkm1 * z - pkm2 * yc;
qk = qkm1 * z - qkm2 * yc;
if( qk != 0 )
{
r = pk/qk;
t = md_fabs( (ans - r)/r );
ans = r;
}
else
t = 1.0;
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
if( md_fabs(pk) > big )
{
pkm2 *= biginv;
pkm1 *= biginv;
qkm2 *= biginv;
qkm1 *= biginv;
}
}
while( t > MACHEP );
return( ans * ax );
}
/* left tail of incomplete md_gamma function:
*
* inf. k
* a -x - x
* x e > ----------
* - -
* k=0 | (a+k+1)
*
*/
double igam( a, x )
double a, x;
{
double ans, ax, c, r;
if( (x <= 0) || ( a <= 0) )
return( 0.0 );
if( (x > 1.0) && (x > a ) )
return( 1.0 - igamc(a,x) );
/* Compute x**a * md_exp(-x) / md_gamma(a) */
ax = a * md_log(x) - x - lgam(a);
if( ax < -MAXLOG )
{
mtherr( "igam", UNDERFLOW );
return( 0.0 );
}
ax = md_exp(ax);
/* power series */
r = a;
c = 1.0;
ans = 1.0;
do
{
r += 1.0;
c *= x/r;
ans += c;
}
while( c/ans > MACHEP );
return( ans * ax/a );
}