GammaSeries.cc

Go to the documentation of this file.

#include "CommonTools/Statistics/src/GammaSeries.h"
#include <iostream>
#include <cmath>

#define ITMAX 100             // maximum allowed number of iterations
#define EPS 3.0e-7            // relative accuracy

float GammaSeries( float a, float x )
{
  if( x < 0.0 ) 
    std::cerr << "GammaSeries::negative argument x" << std::endl;

  if( x == 0. )
    return 0.;

  if( a == 0. ) // this happens at the end, but save all the iterations
    return 0.;

  // coefficient c_n of x^n is Gamma(a)/Gamma(a+1+n), which leads to the
  // recurrence relation c_n = c_(n-1) / (a+n-1) with c_0 = 1/a
  double term = 1/a;
  double sum = term;
  double aplus = a;
  for( int index = 1; index <= ITMAX; index++) {
    ++aplus;
    term *= x/aplus;
    sum += term;
    if( fabs(term) < fabs(sum)*EPS )
      // global coefficient e^-x * x^a / Gamma(a)
      return sum;
  }
  std::cerr << "GammaSeries::a too large, ITMAX too small" << std::endl;
  return 0.;
}