CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

IncompleteGamma.cc
Go to the documentation of this file.
1// -*- C++ -*-
2// $Id: IncompleteGamma.cc,v 1.4 2003/10/10 17:40:39 garren Exp $
3
5#include <assert.h>
6#include <cmath>
7using namespace std;
8
9namespace Genfun {
10FUNCTION_OBJECT_IMP(IncompleteGamma)
11
12const int IncompleteGamma::ITMAX =100;
13const double IncompleteGamma::EPS =3.0E-7;
14const double IncompleteGamma::FPMIN =1.0e-30;
15
16
18 _a("a", 1.0, 0,10)
19{}
20
22AbsFunction(right), _a(right._a) {
23}
24
27
28double IncompleteGamma::operator() (double x) const {
29
30 assert(x>=0.0 && _a.getValue() > 0.0);
31
32 if (x < (_a.getValue()+1.0))
33 return _gamser(_a.getValue(),x,_logGamma(_a.getValue()));
34 else
35 return 1.0-_gammcf(_a.getValue(),x,_logGamma(_a.getValue()));
36}
37
39 return _a;
40}
41
42/* ------------------Incomplete gamma function-----------------*/
43/* ------------------via its series representation-------------*/
44double IncompleteGamma::_gamser(double xa, double x, double logGamma) const {
45 double n;
46 double ap,del,sum;
47
48 ap=xa;
49 del=sum=1.0/xa;
50 for (n=1;n<ITMAX;n++) {
51 ++ap;
52 del *= x/ap;
53 sum += del;
54 if (fabs(del) < fabs(sum)*EPS) return sum*exp(-x + xa*log(x) - logGamma);
55 }
56 assert(0);
57 return 0;
58}
59
60/* ------------------Incomplete gamma function complement------*/
61/* ------------------via its continued fraction representation-*/
62
63double IncompleteGamma::_gammcf(double xa, double x, double logGamma) const {
64
65 double an,b,c,d,del,h;
66 int i;
67
68 b = x + 1.0 -xa;
69 c = 1.0/FPMIN;
70 d = 1.0/b;
71 h = d;
72 for (i=1;i<ITMAX;i++) {
73 an = -i*(i-xa);
74 b+=2.0;
75 d=an*d+b;
76 if (fabs(d) < FPMIN) d = FPMIN;
77 c = b+an/c;
78 if (fabs(c) < FPMIN) c = FPMIN;
79 d = 1.0/d;
80 del=d*c;
81 h *= del;
82 if (fabs(del-1.0) < EPS) return exp(-x+xa*log(x)-logGamma)*h;
83 }
84 assert(0);
85 return 0;
86}
87
88
89
90
91} // namespace Genfun
#define FUNCTION_OBJECT_IMP(classname)
virtual double operator()(double argument) const
virtual double getValue() const
Definition Parameter.cc:27