6 * AS245, 2nd algorithm, http://lib.stat.cmu.edu/apstat/245
8 double kf_lgamma(double z)
11 x += 0.1659470187408462e-06 / (z+7);
12 x += 0.9934937113930748e-05 / (z+6);
13 x -= 0.1385710331296526 / (z+5);
14 x += 12.50734324009056 / (z+4);
15 x -= 176.6150291498386 / (z+3);
16 x += 771.3234287757674 / (z+2);
17 x -= 1259.139216722289 / (z+1);
18 x += 676.5203681218835 / z;
19 x += 0.9999999999995183;
20 return log(x) - 5.58106146679532777 - z + (z-0.5) * log(z+6.5);
23 /* complementary error function
24 * \frac{2}{\sqrt{\pi}} \int_x^{\infty} e^{-t^2} dt
25 * AS66, 2nd algorithm, http://lib.stat.cmu.edu/apstat/66
27 double kf_erfc(double x)
29 const double p0 = 220.2068679123761;
30 const double p1 = 221.2135961699311;
31 const double p2 = 112.0792914978709;
32 const double p3 = 33.912866078383;
33 const double p4 = 6.37396220353165;
34 const double p5 = .7003830644436881;
35 const double p6 = .03526249659989109;
36 const double q0 = 440.4137358247522;
37 const double q1 = 793.8265125199484;
38 const double q2 = 637.3336333788311;
39 const double q3 = 296.5642487796737;
40 const double q4 = 86.78073220294608;
41 const double q5 = 16.06417757920695;
42 const double q6 = 1.755667163182642;
43 const double q7 = .08838834764831844;
45 z = fabs(x) * M_SQRT2;
46 if (z > 37.) return x > 0.? 0. : 2.;
47 expntl = exp(z * z * - .5);
48 if (z < 10. / M_SQRT2) // for small z
49 p = expntl * ((((((p6 * z + p5) * z + p4) * z + p3) * z + p2) * z + p1) * z + p0)
50 / (((((((q7 * z + q6) * z + q5) * z + q4) * z + q3) * z + q2) * z + q1) * z + q0);
51 else p = expntl / 2.506628274631001 / (z + 1. / (z + 2. / (z + 3. / (z + 4. / (z + .65)))));
52 return x > 0.? 2. * p : 2. * (1. - p);
55 /* The following computes regularized incomplete gamma functions.
56 * Formulas are taken from Wiki, with additional input from Numerical
57 * Recipes in C (for modified Lentz's algorithm) and AS245
58 * (http://lib.stat.cmu.edu/apstat/245).
61 #define KF_GAMMA_EPS 1e-14
62 #define KF_TINY 1e-290
64 // regularized lower incomplete gamma function, by series expansion
65 static double _kf_gammap(double s, double z)
69 for (k = 1, sum = x = 1.; k < 100; ++k) {
70 sum += (x *= z / (s + k));
71 if (x / sum < KF_GAMMA_EPS) break;
73 return exp(s * log(z) - z - kf_lgamma(s + 1.) + log(sum));
75 // regularized upper incomplete gamma function, by continued fraction
76 static double _kf_gammaq(double s, double z)
80 f = 1. + z - s; C = f; D = 0.;
81 // Modified Lentz's algorithm for computing continued fraction
82 // See Numerical Recipes in C, 2nd edition, page 172
83 for (j = 1; j < 100; ++j) {
84 double a = j * (s - j), b = (j<<1) + 1 + z - s, d;
86 if (D < KF_TINY) D = KF_TINY;
88 if (C < KF_TINY) C = KF_TINY;
92 if (fabs(d - 1.) < KF_GAMMA_EPS) break;
94 return exp(s * log(z) - z - kf_lgamma(s) - log(f));
97 double kf_gammap(double s, double z)
99 return z <= 1. || z < s? _kf_gammap(s, z) : 1. - _kf_gammaq(s, z);
102 double kf_gammaq(double s, double z)
104 return z <= 1. || z < s? 1. - _kf_gammap(s, z) : _kf_gammaq(s, z);
109 int main(int argc, char *argv[])
111 double x = 5.5, y = 3;
112 printf("erfc(%lg): %lg, %lg\n", x, erfc(x), kf_erfc(x));
113 printf("lower-gamma(%lg,%lg): %lg\n", x, y, (1.0-kf_gammap(y, x))*tgamma(y));