26 double x, y, tmp, ser;
27 static double cof[6]={76.18009172947146,-86.50532032941677,
28 24.01409824083091,-1.231739572450155,
29 0.1208650973866179e-2,-0.5395239384953e-5};
33 tmp -= ( x + 0.5 ) * log( tmp );
34 ser = 1.000000000190015;
36 for (
int j = 0; j <= 5; j++)
37 ser += cof[j] / (++y);
39 return -tmp + log( 2.5066282746310005 * ser / x);
43 void gser(
double *gamser,
double a,
double x,
double *gln)
52 if (x < 0.0) cerr <<
"Error: x less than 0 in routine gser" << endl;
60 for (n = 1; n <=
ITMAX; n++)
66 if (abs(del) < abs(sum)*
EPS)
68 *gamser = sum * exp( -x + a * log( x ) -( *gln ) );
73 cerr <<
"ERROR: a too large, ITMAX too small in routine gser" << endl;
78 void gcf(
double *gammcf,
double a,
double x,
double *gln )
81 double an,b,c,d,del,h;
89 for( i = 1; i <=
ITMAX; i++ )
95 if ( abs( d ) <
FPMIN )
100 if ( abs( c ) <
FPMIN )
107 if ( abs( del - 1.0 ) <
EPS)
111 if (i >
ITMAX) cerr <<
"a too large, ITMAX too small in gcf" << endl;
113 *gammcf = exp( -x + a * log( x ) -( *gln ) ) * h;
118 void gcf(
double *gammcf,
double a,
double x,
double *gln);
120 double gamser,gammcf,gln;
122 if (x < 0.0 || a <= 0.0)
124 cerr <<
"a = " << a <<
" x = " << x << endl;
125 cerr <<
"Invalid arguments in routine GAMMQ" << endl;
131 gser(&gamser,a,x,&gln);
136 gcf(&gammcf,a,x,&gln);