CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: IncompleteGamma.cc,v 1.4 2003/10/10 17:40:39 garren Exp $ 00003 00004 #include "CLHEP/GenericFunctions/IncompleteGamma.hh" 00005 #include <assert.h> 00006 #include <cmath> 00007 using namespace std; 00008 00009 namespace Genfun { 00010 FUNCTION_OBJECT_IMP(IncompleteGamma) 00011 00012 const int IncompleteGamma::ITMAX =100; 00013 const double IncompleteGamma::EPS =3.0E-7; 00014 const double IncompleteGamma::FPMIN =1.0e-30; 00015 00016 00017 IncompleteGamma::IncompleteGamma(): 00018 _a("a", 1.0, 0,10) 00019 {} 00020 00021 IncompleteGamma::IncompleteGamma(const IncompleteGamma & right): 00022 _a(right._a) { 00023 } 00024 00025 IncompleteGamma::~IncompleteGamma() { 00026 } 00027 00028 double IncompleteGamma::operator() (double x) const { 00029 00030 assert(x>=0.0 && _a.getValue() > 0.0); 00031 00032 if (x < (_a.getValue()+1.0)) 00033 return _gamser(_a.getValue(),x,_logGamma(_a.getValue())); 00034 else 00035 return 1.0-_gammcf(_a.getValue(),x,_logGamma(_a.getValue())); 00036 } 00037 00038 Parameter & IncompleteGamma::a() { 00039 return _a; 00040 } 00041 00042 /* ------------------Incomplete gamma function-----------------*/ 00043 /* ------------------via its series representation-------------*/ 00044 double IncompleteGamma::_gamser(double a, double x, double logGamma) const { 00045 double n; 00046 double ap,del,sum; 00047 00048 ap=a; 00049 del=sum=1.0/a; 00050 for (n=1;n<ITMAX;n++) { 00051 ++ap; 00052 del *= x/ap; 00053 sum += del; 00054 if (fabs(del) < fabs(sum)*EPS) return sum*exp(-x + a*log(x) - logGamma); 00055 } 00056 assert(0); 00057 return 0; 00058 } 00059 00060 /* ------------------Incomplete gamma function complement------*/ 00061 /* ------------------via its continued fraction representation-*/ 00062 00063 double IncompleteGamma::_gammcf(double a, double x, double logGamma) const { 00064 00065 double an,b,c,d,del,h; 00066 int i; 00067 00068 b = x + 1.0 -a; 00069 c = 1.0/FPMIN; 00070 d = 1.0/b; 00071 h = d; 00072 for (i=1;i<ITMAX;i++) { 00073 an = -i*(i-a); 00074 b+=2.0; 00075 d=an*d+b; 00076 if (fabs(d) < FPMIN) d = FPMIN; 00077 c = b+an/c; 00078 if (fabs(c) < FPMIN) c = FPMIN; 00079 d = 1.0/d; 00080 del=d*c; 00081 h *= del; 00082 if (fabs(del-1.0) < EPS) return exp(-x+a*log(x)-logGamma)*h; 00083 } 00084 assert(0); 00085 return 0; 00086 } 00087 00088 00089 00090 00091 } // namespace Genfun