115 lines
2.4 KiB
C++
115 lines
2.4 KiB
C++
#include "../rng.h"
|
|
|
|
#pragma once
|
|
|
|
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;
|
|
}
|
|
|