#include "../rng.h" #include "./incomplete_gamma.h" double igam(double a, double x) { double maxlog = 7.09782712893383996732E2; double machep = 1.11022302462515654042E-16; double big = 4.503599627370496e15; double biginv = 2.22044604925031308085e-16; double ans, ax, c, r; if( (x<=0 || a<=0)){ return 0.0; } // if (x>1.0 && x>a){ // return 1.0 - igamc(a,x); // } ax = a * std::log(x) - x - std::lgamma(a); if (ax < -maxlog){ // std::cout << "max log error on incomplete gamma function\n"; return 0.0; } ax = std::exp(ax); r = a; c = 1.0; ans = 1.0; do { r+=1; c*= x/r; ans +=c; }while(c/ans > machep); return ans * ax/a; } // regularized left tail of incomplete gamma // Q(s, x) = \frac{1}{\Gamma(s)} \int_x^\infty t^{s-1} e^{-t} \, dt double igamc(double a, double x){ return 1-igam(a,x); // double maxlog = 7.09782712893383996732E2; // double machep = 1.11022302462515654042E-16; // double big = 4.503599627370496e15; // double biginv = 2.22044604925031308085e-16; // if(x<=0 || a<=0){ // return 1.0; // } // if( (x < 1.0 ) || ( x < a)){ // std::cout << "bad\n"; // return 1.0 - igam(a,x); // } // double ax = a * std::log(x) - std::lgamma(a); // // maxlog? // if(ax < -maxlog){ // std::cout << "max log error on incomplete gamma function\n"; // return 0.0; // } // ax = std::exp(ax); // double y = 1.0 - a; // double z = x + y + 1.0; // double c = 0.0; // double pkm2 = 1.0; // double qkm2 = x; // double pkm1 = x + 1.0; // double qkm1 = z * x; // double ans = pkm1/qkm1; // double pk,qk,r,t,yc; // t = 1.0; // do { // c+=1.0; // y+=1.0; // z+=2.0; // yc = y*c; // pk = pkm1 * z - pkm2 * yc; // qk = qkm1 * z - qkm2 * yc; // if (qk !=0){ // r = pk/qk; // t = std::fabs((ans-r)/r); // ans = r; // }else{ // t = 1.0; // } // pkm2 = pkm1; // pkm1 = pk; // qkm2 = qkm1; // qkm1 = qk; // if(std::fabs(pk) > big){ // pkm2 *= biginv; // pkm1 *= biginv; // qkm2 *= biginv; // qkm1 *= biginv; // } // std::cout << ans << "-\n"; // } while (t>machep); // return ans * ax; }