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>
7 using namespace std;
8 
9 namespace Genfun {
10 FUNCTION_OBJECT_IMP(IncompleteGamma)
11 
12 const int IncompleteGamma::ITMAX =100;
13 const double IncompleteGamma::EPS =3.0E-7;
14 const double IncompleteGamma::FPMIN =1.0e-30;
15 
16 
17 IncompleteGamma::IncompleteGamma():
18  _a("a", 1.0, 0,10)
19 {}
20 
22 AbsFunction(right), _a(right._a) {
23 }
24 
26 }
27 
28 double 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-------------*/
44 double 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 
63 double 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
STL namespace.
virtual double getValue() const
Definition: Parameter.cc:27
#define FUNCTION_OBJECT_IMP(classname)
virtual double operator()(double argument) const