1 /* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 /* This Source Code Form is subject to the terms of the Mozilla Public
3 * License, v. 2.0. If a copy of the MPL was not distributed with this
4 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
5
6 #ifndef nsIncompleteGamma_h__
7 #define nsIncompleteGamma_h__
8
9 /* An implementation of the incomplete gamma functions for real
10 arguments. P is defined as
11
12 x
13 /
14 1 [ a - 1 - t
15 P(a, x) = -------- I t e dt
16 Gamma(a) ]
17 /
18 0
19
20 and
21
22 infinity
23 /
24 1 [ a - 1 - t
25 Q(a, x) = -------- I t e dt
26 Gamma(a) ]
27 /
28 x
29
30 so that P(a,x) + Q(a,x) = 1.
31
32 Both a series expansion and a continued fraction exist. This
33 implementation uses the more efficient method based on the arguments.
34
35 Either case involves calculating a multiplicative term:
36 e^(-x)*x^a/Gamma(a).
37 Here we calculate the log of this term. Most math libraries have a
38 "lgamma" function but it is not re-entrant. Some libraries have a
39 "lgamma_r" which is re-entrant. Use it if possible. I have included a
40 simple replacement but it is certainly not as accurate.
41
42 Relative errors are almost always < 1e-10 and usually < 1e-14. Very
43 small and very large arguments cause trouble.
44
45 The region where a < 0.5 and x < 0.5 has poor error properties and is
46 not too stable. Get a better routine if you need results in this
47 region.
48
49 The error argument will be set negative if there is a domain error or
50 positive for an internal calculation error, currently lack of
51 convergence. A value is always returned, though.
52
53 */
54
55 #include <math.h>
56 #include <float.h>
57
58 // the main routine
59 static double nsIncompleteGammaP(double a, double x, int* error);
60
61 // nsLnGamma(z): either a wrapper around lgamma_r or the internal function.
62 // C_m = B[2*m]/(2*m*(2*m-1)) where B is a Bernoulli number
63 static const double C_1 = 1.0 / 12.0;
64 static const double C_2 = -1.0 / 360.0;
65 static const double C_3 = 1.0 / 1260.0;
66 static const double C_4 = -1.0 / 1680.0;
67 static const double C_5 = 1.0 / 1188.0;
68 static const double C_6 = -691.0 / 360360.0;
69 static const double C_7 = 1.0 / 156.0;
70 static const double C_8 = -3617.0 / 122400.0;
71 static const double C_9 = 43867.0 / 244188.0;
72 static const double C_10 = -174611.0 / 125400.0;
73 static const double C_11 = 77683.0 / 5796.0;
74
75 // truncated asymptotic series in 1/z
lngamma_asymp(double z)76 static inline double lngamma_asymp(double z) {
77 double w, w2, sum;
78 w = 1.0 / z;
79 w2 = w * w;
80 sum =
81 w *
82 (w2 * (w2 * (w2 * (w2 * (w2 * (w2 * (w2 * (w2 * (w2 * (C_11 * w2 + C_10) +
83 C_9) +
84 C_8) +
85 C_7) +
86 C_6) +
87 C_5) +
88 C_4) +
89 C_3) +
90 C_2) +
91 C_1);
92
93 return sum;
94 }
95
96 struct fact_table_s {
97 double fact;
98 double lnfact;
99 };
100
101 // for speed and accuracy
102 static const struct fact_table_s FactTable[] = {
103 {1.000000000000000, 0.0000000000000000000000e+00},
104 {1.000000000000000, 0.0000000000000000000000e+00},
105 {2.000000000000000, 6.9314718055994530942869e-01},
106 {6.000000000000000, 1.7917594692280550007892e+00},
107 {24.00000000000000, 3.1780538303479456197550e+00},
108 {120.0000000000000, 4.7874917427820459941458e+00},
109 {720.0000000000000, 6.5792512120101009952602e+00},
110 {5040.000000000000, 8.5251613610654142999881e+00},
111 {40320.00000000000, 1.0604602902745250228925e+01},
112 {362880.0000000000, 1.2801827480081469610995e+01},
113 {3628800.000000000, 1.5104412573075515295248e+01},
114 {39916800.00000000, 1.7502307845873885839769e+01},
115 {479001600.0000000, 1.9987214495661886149228e+01},
116 {6227020800.000000, 2.2552163853123422886104e+01},
117 {87178291200.00000, 2.5191221182738681499610e+01},
118 {1307674368000.000, 2.7899271383840891566988e+01},
119 {20922789888000.00, 3.0671860106080672803835e+01},
120 {355687428096000.0, 3.3505073450136888885825e+01},
121 {6402373705728000., 3.6395445208033053576674e+01}};
122 #define FactTableLength (int)(sizeof(FactTable) / sizeof(FactTable[0]))
123
124 // for speed
125 static const double ln_2pi_2 = 0.918938533204672741803; // log(2*PI)/2
126
127 /* A simple lgamma function, not very robust.
128
129 Valid for z_in > 0 ONLY.
130
131 For z_in > 8 precision is quite good, relative errors < 1e-14 and
132 usually better. For z_in < 8 relative errors increase but are usually
133 < 1e-10. In two small regions, 1 +/- .001 and 2 +/- .001 errors
134 increase quickly.
135 */
nsLnGamma(double z_in,int * gsign)136 static double nsLnGamma(double z_in, int* gsign) {
137 double scale, z, sum, result;
138 *gsign = 1;
139
140 int zi = (int)z_in;
141 if (z_in == (double)zi) {
142 if (0 < zi && zi <= FactTableLength)
143 return FactTable[zi - 1].lnfact; // gamma(z) = (z-1)!
144 }
145
146 for (scale = 1.0, z = z_in; z < 8.0; ++z) scale *= z;
147
148 sum = lngamma_asymp(z);
149 result = (z - 0.5) * log(z) - z + ln_2pi_2 - log(scale);
150 result += sum;
151 return result;
152 }
153
154 // log( e^(-x)*x^a/Gamma(a) )
lnPQfactor(double a,double x)155 static inline double lnPQfactor(double a, double x) {
156 int gsign; // ignored because a > 0
157 return a * log(x) - x - nsLnGamma(a, &gsign);
158 }
159
Pseries(double a,double x,int * error)160 static double Pseries(double a, double x, int* error) {
161 double sum, term;
162 const double eps = 2.0 * DBL_EPSILON;
163 const int imax = 5000;
164 int i;
165
166 sum = term = 1.0 / a;
167 for (i = 1; i < imax; ++i) {
168 term *= x / (a + i);
169 sum += term;
170 if (fabs(term) < eps * fabs(sum)) break;
171 }
172
173 if (i >= imax) *error = 1;
174
175 return sum;
176 }
177
Qcontfrac(double a,double x,int * error)178 static double Qcontfrac(double a, double x, int* error) {
179 double result, D, C, e, f, term;
180 const double eps = 2.0 * DBL_EPSILON;
181 const double small = DBL_EPSILON * DBL_EPSILON * DBL_EPSILON * DBL_EPSILON;
182 const int imax = 5000;
183 int i;
184
185 // modified Lentz method
186 f = x - a + 1.0;
187 if (fabs(f) < small) f = small;
188 C = f + 1.0 / small;
189 D = 1.0 / f;
190 result = D;
191 for (i = 1; i < imax; ++i) {
192 e = i * (a - i);
193 f += 2.0;
194 D = f + e * D;
195 if (fabs(D) < small) D = small;
196 D = 1.0 / D;
197 C = f + e / C;
198 if (fabs(C) < small) C = small;
199 term = C * D;
200 result *= term;
201 if (fabs(term - 1.0) < eps) break;
202 }
203
204 if (i >= imax) *error = 1;
205 return result;
206 }
207
nsIncompleteGammaP(double a,double x,int * error)208 static double nsIncompleteGammaP(double a, double x, int* error) {
209 double result, dom, ldom;
210 // domain errors. the return values are meaningless but have
211 // to return something.
212 *error = -1;
213 if (a <= 0.0) return 1.0;
214 if (x < 0.0) return 0.0;
215 *error = 0;
216 if (x == 0.0) return 0.0;
217
218 ldom = lnPQfactor(a, x);
219 dom = exp(ldom);
220 // might need to adjust the crossover point
221 if (a <= 0.5) {
222 if (x < a + 1.0)
223 result = dom * Pseries(a, x, error);
224 else
225 result = 1.0 - dom * Qcontfrac(a, x, error);
226 } else {
227 if (x < a)
228 result = dom * Pseries(a, x, error);
229 else
230 result = 1.0 - dom * Qcontfrac(a, x, error);
231 }
232
233 // not clear if this can ever happen
234 if (result > 1.0) result = 1.0;
235 if (result < 0.0) result = 0.0;
236 return result;
237 }
238
239 #endif
240