Statistics/src/IncompleteGammaComplement.cc

0001 #include "CommonTools/Statistics/src/IncompleteGammaComplement.h"
0002 #include "CommonTools/Statistics/src/GammaContinuedFraction.h"
0003 #include "CommonTools/Statistics/src/GammaSeries.h"
0004 #include "CommonTools/Statistics/src/GammaLn.h"
0005
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include <iostream>
0009 #include <cmath>
0010
0011 float IncompleteGammaComplement::value(float a, float x) {
0012   if (x < 0.0 || a <= 0.0)
0013     edm::LogInfo("IncompleteGammaComplement") << "IncompleteGammaComplement::invalid arguments";
0014   if (x < (a + 1.0))
0015     // take the complement of the series representation
0016     return 1. - GammaSeries(a, x) * (exp(-x + a * log(x) - GammaLn(a)));
0017   else
0018     // use the continued fraction representation
0019     return GammaContinuedFraction(a, x) * (exp(-x + a * log(x) - GammaLn(a)));
0020 }
0021
0022 float IncompleteGammaComplement::ln(float a, float x) {
0023   if (x < 0.0 || a <= 0.0)
0024     edm::LogInfo("IncompleteGammaComplement") << "IncompleteGammaComplement::invalid arguments";
0025   if (x < (a + 1.0))
0026     // take the complement of the series representation
0027     return log(1. - GammaSeries(a, x) * (exp(-x + a * log(x) - GammaLn(a))));
0028   else
0029     // use the continued fraction representation
0030     return log(GammaContinuedFraction(a, x)) - x + a * log(x) - GammaLn(a);
0031 }