Statistics/src/GammaSeries.cc

0001 #include "CommonTools/Statistics/src/GammaSeries.h"
0002 #include <iostream>
0003 #include <cmath>
0004
0005 #define ITMAX 100   // maximum allowed number of iterations
0006 #define EPS 3.0e-7  // relative accuracy
0007
0008 float GammaSeries(float a, float x) {
0009   if (x < 0.0)
0010     std::cerr << "GammaSeries::negative argument x" << std::endl;
0011
0012   if (x == 0.)
0013     return 0.;
0014
0015   if (a == 0.)  // this happens at the end, but save all the iterations
0016     return 0.;
0017
0018   // coefficient c_n of x^n is Gamma(a)/Gamma(a+1+n), which leads to the
0019   // recurrence relation c_n = c_(n-1) / (a+n-1) with c_0 = 1/a
0020   double term = 1 / a;
0021   double sum = term;
0022   double aplus = a;
0023   for (int index = 1; index <= ITMAX; index++) {
0024     ++aplus;
0025     term *= x / aplus;
0026     sum += term;
0027     if (fabs(term) < fabs(sum) * EPS)
0028       // global coefficient e^-x * x^a / Gamma(a)
0029       return sum;
0030   }
0031   std::cerr << "GammaSeries::a too large, ITMAX too small" << std::endl;
0032   return 0.;
0033 }