00001
00010 #include "Gammaq.h"
00011
00012 #include <iostream>
00013
00014 #include <cassert>
00015 #include <cmath>
00016
00017 using std::cerr;
00018 using std::endl;
00019 using std::abs;
00020
00021 namespace hippodraw {
00022 namespace Numeric {
00023
00024 double gammln( double xx )
00025 {
00026 double x, y, tmp, ser;
00027 static double cof[6]={76.18009172947146,-86.50532032941677,
00028 24.01409824083091,-1.231739572450155,
00029 0.1208650973866179e-2,-0.5395239384953e-5};
00030
00031 y = x = xx;
00032 tmp = x + 5.5;
00033 tmp -= ( x + 0.5 ) * log( tmp );
00034 ser = 1.000000000190015;
00035
00036 for ( int j = 0; j <= 5; j++)
00037 ser += cof[j] / (++y);
00038
00039 return -tmp + log( 2.5066282746310005 * ser / x);
00040 }
00041
00042
00043 void gser(double *gamser, double a, double x, double *gln)
00044 {
00045 int n;
00046 double sum, del, ap;
00047
00048 *gln = gammln( a );
00049
00050 if ( x <= 0.0 )
00051 {
00052 if (x < 0.0) cerr << "Error: x less than 0 in routine gser" << endl;
00053 *gamser=0.0;
00054 return;
00055 }
00056 else
00057 {
00058 ap = a;
00059 del = sum =1.0 / a;
00060 for (n = 1; n <= ITMAX; n++)
00061 {
00062 ++ap;
00063 del *= x/ap;
00064 sum += del;
00065
00066 if (abs(del) < abs(sum)*EPS)
00067 {
00068 *gamser = sum * exp( -x + a * log( x ) -( *gln ) );
00069 return;
00070 }
00071 }
00072
00073 cerr << "ERROR: a too large, ITMAX too small in routine gser" << endl;
00074 return;
00075 }
00076 }
00077
00078 void gcf( double *gammcf, double a, double x, double *gln )
00079 {
00080 int i;
00081 double an,b,c,d,del,h;
00082
00083 *gln = gammln( a );
00084 b = x + 1.0 - a;
00085 c = 1.0 / FPMIN;
00086 d = 1.0 / b;
00087 h = d;
00088
00089 for( i = 1; i <= ITMAX; i++ )
00090 {
00091 an = -i*(i-a);
00092 b += 2.0;
00093 d=an*d+b;
00094
00095 if ( abs( d ) < FPMIN )
00096 d = FPMIN;
00097
00098 c = b + an / c;
00099
00100 if ( abs( c ) < FPMIN )
00101 c = FPMIN;
00102
00103 d = 1.0 / d;
00104 del = d * c;
00105 h *= del;
00106
00107 if ( abs( del - 1.0 ) < EPS)
00108 break;
00109 }
00110
00111 if (i > ITMAX) cerr << "a too large, ITMAX too small in gcf" << endl;
00112
00113 *gammcf = exp( -x + a * log( x ) -( *gln ) ) * h;
00114 }
00115
00116 double gammq(double a, double x)
00117 {
00118 void gcf(double *gammcf, double a, double x, double *gln);
00119
00120 double gamser,gammcf,gln;
00121
00122 if (x < 0.0 || a <= 0.0)
00123 {
00124 cerr << "a = " << a << " x = " << x << endl;
00125 cerr << "Invalid arguments in routine GAMMQ" << endl;
00126 assert(0);
00127 }
00128
00129 if (x < (a+1.0))
00130 {
00131 gser(&gamser,a,x,&gln);
00132 return 1.0-gamser;
00133 }
00134 else
00135 {
00136 gcf(&gammcf,a,x,&gln);
00137 return gammcf;
00138 }
00139 }
00140
00141 }
00142 }
00143
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186 #undef ITMAX
00187 #undef EPS
00188 #undef FPMIN
00189 #undef MAXSTR