CMSSW Reference Manual

00001 #include "CommonTools/Statistics/src/GammaContinuedFraction.h"
00002 #include "CommonTools/Statistics/src/GammaLn.h"
00003 #include <cmath>
00004 #include <iostream>
00005 
00006 #define ITMAX 100        // maximum allowed number of iterations
00007 #define EPS 3.0e-7       // relative accuracy
00008 #define FPMIN 1.0e-30    // number near the smallest representable floating-point number
00009 
00010 float
00011 GammaContinuedFraction( float a, float x )
00012 {
00013   int i;
00014   float an,del;
00015 
00016   /* Set up for evaluating continued fraction by modified Lentz's method (par.5.2
00017      in Numerical Recipes in C) with b_0 = 0 */
00018   double b = x+1.0-a;
00019   double c = 1.0/FPMIN;
00020   double d = 1.0/b;
00021   double h = d;
00022   for (i=1;i<=ITMAX;i++) {
00023     an = -i*(i-a);
00024     b += 2.0;
00025     d=an*d+b;
00026     if (fabs(d) < FPMIN) d=FPMIN;
00027     c=b+an/c;
00028     if (fabs(c) < FPMIN) c=FPMIN;
00029     d=1.0/d;
00030     del=d*c;
00031     h *= del;
00032     if (fabs(del-1.0) < EPS) break;
00033   }
00034   if( i > ITMAX ) std::cerr << "GammaContinuedFraction::a too large, "
00035                        << "ITMAX too small" << std::endl;
00036   return h;
00037 }
00038 #undef ITMAX
00039 #undef EPS
00040 #undef FPMIN
00041 /* (C) Copr. 1986-92 Numerical Recipes Software B2.. */
GammaContinuedFraction.cc
