#include #include #define LogRootTwoPi_ 0.9189385332046727418 #define FFMAX(a,b) ((a) > (b) ? (a) : (b)) #define FFMIN(a,b) ((a) > (b) ? (b) : (a)) /* coefficients for gamma=7, kmax=8 Lanczos method */ double gamma_lanczos_rational(double x) { double Ag, denom = (40320+x*(109584+x*(118124+x*(67284+x*(22449+x*(4536+x*(546+x*(36+x)))))))); Ag = (((((((24.45833333334874451493*x + 1053.042534722144725)*x + 19629.51066503033455781)*x + 205786.215499718336944265)*x + 1312212.26587868165589257)*x + 5092745.01776084957)*x + 11137345.4846599424858181)*x + 10579290.0990759986696); Ag = 0.99999999999980993 + (Ag/denom); return Ag; } static const double p[7] = { 1.0e0,0.553413866010467e+00,0.279648642639792e+00,0.730981088720487e-01, 0.204493667594920e-01,0.261939260042690e-02,0.539637273585445e-03, }; static const double q[7] = { 1.0e0,0.113062953091122e+01,-0.567902761974940e-01,-0.170458969313360e+00, 0.225211131035340e-01,0.470059485860584e-02,-0.832979206704073e-03, }; /* Lanczos method for real x > 0; * gamma=7, truncated at 1/(z+8) * [J. SIAM Numer. Anal, Ser. B, 1 (1964) 86] */ static double sf_lngamma(double x) { /* Lanczos writes z! instead of Gamma(z) */ double Ag = gamma_lanczos_rational(x-1.0); /* (x-0.5)*log(x+6.5) - (x+6.5) + LogRootTwoPi_ + log(Ag(x-1.0)) */ double term1 = (x-0.5)*log(x+6.5) - (x+6.5) + LogRootTwoPi_; return term1 + log(Ag); } static const double w1 = 8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */ w2 = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */ w3 = 7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */ w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */ w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */ w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */ double sf_lngamma_asymp(double a) { double ainv = 1.0/a; double z = ainv*ainv; double w = ainv*(w1+z*(w2+z*(w3+z*(w4+z*(w5+z*w6))))); return w; } /* COMPUTE 1/GAMMA(1 + A) FOR 0 .LE. A .LT. 1 */ static double gamma1(double a) { double n = p[0]+a*(p[1]+a*(p[2]+a*(p[3]+a*(p[4]+a*(p[5]+p[6]*a))))); double d = q[0]+a*(q[1]+a*(q[2]+a*(q[3]+a*(q[4]+a*(q[5]+q[6]*a))))); return d/n; } static const double as[4] = {3.31125922108741, 11.6616720288968, 4.28342155967104, 0.213623493715853}; static const double bs[5] = {1, 6.61053765625462, 6.40691597760039, 1.27364489782223, 0.03611708101884203}; double find_inverse_s(double p) { /* * Computation of the Incomplete Gamma Function Ratios and their Inverse * ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR. * ACM Transactions on Mathematical Software, Vol. 12, No. 4, * December 1986, Pages 377-393. * * See equation 32. */ double n, d, s, t = ((p<0.5)?p:(1.0-p)); t = sqrt(-2 * log(t)); n = (as[0]+t*(as[1]+t*(as[2]+as[3]*t))); d = bs[0]+t*(bs[1]+t*(bs[2]+t*(bs[3]+bs[4]*t))); s = t - n/d; if(p < 0.5) s = -s; return s; } /* The dominant part, * D(a,x) := x^a / Gamma(a+1) */ #if 0 double igamma_D(double a, double x) { double val = 0.0; if (a < 1.0) { val = exp(a*log(x))*gamma1(a); } else if (a < 10.0) { val = exp(a*log(x) - sf_lngamma(a+1)); } else { val = exp(a*log(M_E*x/a) - sf_lngamma_asymp(a))/sqrt(2*M_PI*a); } return val; } #endif /* * Computation of the Incomplete Gamma Function Ratios and their Inverse, Equation 34. */ static double didonato_SN(double a, double x, unsigned int N, double tolerance) { double partial = x / (a + 1), sum = 1.0+partial; unsigned int i; for(i = 2; i <= N; ++i) { partial *= x/(a + i); sum += partial; if(partial < tolerance) { break; } } return sum; } double igamma_inverse_initial_guess(double a, double p, double q) { /* * In order to understand what's going on here, you will * need to refer to: * * Computation of the Incomplete Gamma Function Ratios and their Inverse * ARMIDO R. DIDONATO and ALFRED H. MORRIS, JR. * ACM Transactions on Mathematical Software, Vol. 12, No. 4, * December 1986, Pages 377-393. */ double am1 = (a-1.0); double result; if (a == 1) { if (q > 0.9) { result = -log1p(-p); } else { result = -log(q); } } else if (a < 1) { double g = a/gamma1(a); double b = q * g; if ((b > 0.6) || ((b >= 0.45) && (a >= 0.3))) { /* DiDonato & Morris Eq 21: * * There is a slight variation from DiDonato and Morris here: * the first form given here is unstable when p is close to 1, * making it impossible to compute the inverse of Q(a,x) for small * q. Fortunately the second form works perfectly well in this case. */ double u = exp((-q/a) - 0.5772156649015329); result = u/(1-(u/(a + 1))); } else if ((a < 0.3) && (b >= 0.35)) { /* DiDonato & Morris Eq 22: */ double t = exp(-(0.5772156649015329+b)); result = t * exp(t * exp(t)); } else { double y = -log(b); double c1 = am1*log(y); if ((b > 0.15) || (a >= 0.3)) { /* DiDonato & Morris Eq 23: */ double u = y + c1; result = y + am1*log(u) - log(1 + (1-a)/(1+u)); } else if (b > 0.1) { /* DiDonato & Morris Eq 24: */ double u = y + c1; result = y + am1*log(u) - log((u*u + 2*(3-a)*u + (2-a)*(3-a)) / (u*u + (5-a)*u + 2)); } else { /* DiDonato & Morris Eq 25: */ double c1_2 = c1*c1; double c2 = am1 * (1+c1); double c3 = am1 * ((-0.5*c1 + (a-2))*c1 + 0.5*(3*a-5)); double c4 = am1 * (c1*(c1*((c1/3) - 0.5*(3*a-5)) + ((a-6)*a+7)) + ((11*a-46)*a + 47)/6); double c5 = am1 * (c1_2*(-0.25*c1_2 + (11*a-17)*c1/6 + ((-3*a+13)*a-13)) + ((( 2*a- 25)*a+ 72)*a - 61)*c1/2 + (((25*a-195)*a+477)*a - 379)/12); double yi = 1.0/y; result = y + c1 + yi*(c2 + yi*(c3 + yi*(c4 + c5*yi))); } } } else { /* DiDonato and Morris Eq 31: */ double s = find_inverse_s(p), s2 = s*s; double w, rta = sqrt(a); w = a + s*rta + (s2-1) / 3; w += s*(s2-7)/(36*rta); w -= ((3*s2+7)*s2-16)/(810*a); w += s*((9*s2+256)*s2-433)/(38880*a*rta); if ((a >= 500) && (fabs(1-w/a) < 1e-06)) { result = w; } else if (p > 0.5) { if (w < 3 * a) { result = w; } else { double D = FFMAX(2, a*am1); double lg = sf_lngamma(a); double lb = log(q) + lg; if (lb < -D * 2.3) { /* DiDonato and Morris Eq 25: */ double y = -lb; double c1 = am1*log(y); double c2 = am1*(1+c1); double c1_2 = c1 * c1; double c3 = am1*(-(0.5*c1_2) + (a-2)*c1 + (3*a-5)/2); double c4 = am1*(c1_2*(c1/3 - 0.5*(3*a-5)) + ((a-6)*a+7)*c1 + ((11*a-46)*a+47)/6); double c5 = am1 * (-(0.25*c1_2*c1_2) + c1*(11*a-17)*c1_2/6 + ((-3*a+13)*a-13)*c1_2 + (((2*a-25)*a+72)*a-61)*c1/2 + (((25*a-195)*a+477)*a-379)/12); double yi = 1.0/y; result = y + c1 + yi*(c2 + yi*(c3 + yi*(c4 + c5*yi))); } else { /* DiDonato and Morris Eq 33: */ double u = -lb + am1*log(w) - log1p(-am1/(1+w)); result = -lb + am1*log(u) - log1p(-am1/(1+u)); } } } else { double z = w; double ap1 = a + 1; double ap2 = a + 2; double ap3 = a + 3; double v = log(p) + sf_lngamma(ap1); if (w < 0.15 * ap1) { /* DiDonato and Morris Eq 35: */ //double v = log(p) + sf_lngamma(ap1); z = exp((v + w) / a); s = log1p(z / ap1 * (1 + z / ap2)); z = exp((v + z - s) / a); s = log1p(z / ap1 * (1 + z / ap2)); z = exp((v + z - s) / a); s = log1p(z / ap1 * (1 + z / ap2 * (1 + z / ap3))); z = exp((v + z - s) / a); } //result = z; if ((z <= 0.01 * ap1) || (z > 0.7 * ap1)) { result = z; } else { /* DiDonato and Morris Eq 36: */ double z0 = z, ls = log(didonato_SN(a, z, 27, 1e-4)); z = exp((v + z0 - ls) / a); result = z * (1 - ((z0-z)/(a-z))); //result = z * (1 - (a*log(z) - z - v + ls) / (a-z)); } } } return result; }