#include "CommonTools/Statistics/src/GammaContinuedFraction.h"
#include "CommonTools/Statistics/src/GammaLn.h"
#include <cmath>
#include <iostream>
Go to the source code of this file.
Defines | |
#define | EPS 3.0e-7 |
#define | FPMIN 1.0e-30 |
#define | ITMAX 100 |
Functions | |
float | GammaContinuedFraction (float a, float x) |
Returns the continued fraction summation of the (complement of the) incomplete gamma function P(a,x) evaluated by its continued fraction representation. |
#define EPS 3.0e-7 |
Definition at line 7 of file GammaContinuedFraction.cc.
Referenced by GammaContinuedFraction(), GammaSeries(), RectangularPixelTopology::localPosition(), reco::PreshowerCluster::operator==(), and SingleParticleEvent::subtractEloss().
#define FPMIN 1.0e-30 |
#define ITMAX 100 |
Definition at line 6 of file GammaContinuedFraction.cc.
Referenced by GammaContinuedFraction(), and GammaSeries().
float GammaContinuedFraction | ( | float | a, | |
float | x | |||
) |
Returns the continued fraction summation of the (complement of the) incomplete gamma function P(a,x) evaluated by its continued fraction representation.
source: Numerical Recipes
Definition at line 11 of file GammaContinuedFraction.cc.
References b, c, TestMuL1L2Filter_cff::cerr, d, lat::endl(), EPS, FPMIN, h, i, and ITMAX.
Referenced by IncompleteGammaComplement::ln(), and IncompleteGammaComplement::value().
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 }