Statistics/src/GammaContinuedFraction.cc

0001 #include "CommonTools/Statistics/src/GammaContinuedFraction.h"
0002 #include <cmath>
0003 #include <iostream>
0004
0005 #define ITMAX 100      // maximum allowed number of iterations
0006 #define EPS 3.0e-7     // relative accuracy
0007 #define FPMIN 1.0e-30  // number near the smallest representable floating-point number
0008
0009 float GammaContinuedFraction(float a, float x) {
0010   int i;
0011   float an, del;
0012
0013   /* Set up for evaluating continued fraction by modified Lentz's method (par.5.2
0014      in Numerical Recipes in C) with b_0 = 0 */
0015   double b = x + 1.0 - a;
0016   double c = 1.0 / FPMIN;
0017   double d = 1.0 / b;
0018   double h = d;
0019   for (i = 1; i <= ITMAX; i++) {
0020     an = -i * (i - a);
0021     b += 2.0;
0022     d = an * d + b;
0023     if (fabs(d) < FPMIN)
0024       d = FPMIN;
0025     c = b + an / c;
0026     if (fabs(c) < FPMIN)
0027       c = FPMIN;
0028     d = 1.0 / d;
0029     del = d * c;
0030     h *= del;
0031     if (fabs(del - 1.0) < EPS)
0032       break;
0033   }
0034   if (i > ITMAX)
0035     std::cerr << "GammaContinuedFraction::a too large, "
0036               << "ITMAX too small" << std::endl;
0037   return h;
0038 }
0039 #undef ITMAX
0040 #undef EPS
0041 #undef FPMIN
0042 /* (C) Copr. 1986-92 Numerical Recipes Software B2.. */