Jason Evans | b1941c6 | 2013-12-09 23:36:37 -0800 | [diff] [blame] | 1 | #ifndef JEMALLOC_ENABLE_INLINE |
| 2 | double ln_gamma(double x); |
| 3 | double i_gamma(double x, double p, double ln_gamma_p); |
| 4 | double pt_norm(double p); |
| 5 | double pt_chi2(double p, double df, double ln_gamma_df_2); |
| 6 | double pt_gamma(double p, double shape, double scale, double ln_gamma_shape); |
| 7 | #endif |
| 8 | |
| 9 | #if (defined(JEMALLOC_ENABLE_INLINE) || defined(MATH_C_)) |
| 10 | /* |
| 11 | * Compute the natural log of Gamma(x), accurate to 10 decimal places. |
| 12 | * |
| 13 | * This implementation is based on: |
| 14 | * |
| 15 | * Pike, M.C., I.D. Hill (1966) Algorithm 291: Logarithm of Gamma function |
| 16 | * [S14]. Communications of the ACM 9(9):684. |
| 17 | */ |
| 18 | JEMALLOC_INLINE double |
| 19 | ln_gamma(double x) |
| 20 | { |
| 21 | double f, z; |
| 22 | |
| 23 | assert(x > 0.0); |
| 24 | |
| 25 | if (x < 7.0) { |
| 26 | f = 1.0; |
| 27 | z = x; |
| 28 | while (z < 7.0) { |
| 29 | f *= z; |
| 30 | z += 1.0; |
| 31 | } |
| 32 | x = z; |
| 33 | f = -log(f); |
| 34 | } else |
| 35 | f = 0.0; |
| 36 | |
| 37 | z = 1.0 / (x * x); |
| 38 | |
| 39 | return (f + (x-0.5) * log(x) - x + 0.918938533204673 + |
| 40 | (((-0.000595238095238 * z + 0.000793650793651) * z - |
| 41 | 0.002777777777778) * z + 0.083333333333333) / x); |
| 42 | } |
| 43 | |
| 44 | /* |
| 45 | * Compute the incomplete Gamma ratio for [0..x], where p is the shape |
| 46 | * parameter, and ln_gamma_p is ln_gamma(p). |
| 47 | * |
| 48 | * This implementation is based on: |
| 49 | * |
| 50 | * Bhattacharjee, G.P. (1970) Algorithm AS 32: The incomplete Gamma integral. |
| 51 | * Applied Statistics 19:285-287. |
| 52 | */ |
| 53 | JEMALLOC_INLINE double |
| 54 | i_gamma(double x, double p, double ln_gamma_p) |
| 55 | { |
| 56 | double acu, factor, oflo, gin, term, rn, a, b, an, dif; |
| 57 | double pn[6]; |
| 58 | unsigned i; |
| 59 | |
| 60 | assert(p > 0.0); |
| 61 | assert(x >= 0.0); |
| 62 | |
| 63 | if (x == 0.0) |
| 64 | return (0.0); |
| 65 | |
| 66 | acu = 1.0e-10; |
| 67 | oflo = 1.0e30; |
| 68 | gin = 0.0; |
| 69 | factor = exp(p * log(x) - x - ln_gamma_p); |
| 70 | |
| 71 | if (x <= 1.0 || x < p) { |
| 72 | /* Calculation by series expansion. */ |
| 73 | gin = 1.0; |
| 74 | term = 1.0; |
| 75 | rn = p; |
| 76 | |
| 77 | while (true) { |
| 78 | rn += 1.0; |
| 79 | term *= x / rn; |
| 80 | gin += term; |
| 81 | if (term <= acu) { |
| 82 | gin *= factor / p; |
| 83 | return (gin); |
| 84 | } |
| 85 | } |
| 86 | } else { |
| 87 | /* Calculation by continued fraction. */ |
| 88 | a = 1.0 - p; |
| 89 | b = a + x + 1.0; |
| 90 | term = 0.0; |
| 91 | pn[0] = 1.0; |
| 92 | pn[1] = x; |
| 93 | pn[2] = x + 1.0; |
| 94 | pn[3] = x * b; |
| 95 | gin = pn[2] / pn[3]; |
| 96 | |
| 97 | while (true) { |
| 98 | a += 1.0; |
| 99 | b += 2.0; |
| 100 | term += 1.0; |
| 101 | an = a * term; |
| 102 | for (i = 0; i < 2; i++) |
| 103 | pn[i+4] = b * pn[i+2] - an * pn[i]; |
| 104 | if (pn[5] != 0.0) { |
| 105 | rn = pn[4] / pn[5]; |
| 106 | dif = fabs(gin - rn); |
| 107 | if (dif <= acu && dif <= acu * rn) { |
| 108 | gin = 1.0 - factor * gin; |
| 109 | return (gin); |
| 110 | } |
| 111 | gin = rn; |
| 112 | } |
| 113 | for (i = 0; i < 4; i++) |
| 114 | pn[i] = pn[i+2]; |
| 115 | |
| 116 | if (fabs(pn[4]) >= oflo) { |
| 117 | for (i = 0; i < 4; i++) |
| 118 | pn[i] /= oflo; |
| 119 | } |
| 120 | } |
| 121 | } |
| 122 | } |
| 123 | |
| 124 | /* |
| 125 | * Given a value p in [0..1] of the lower tail area of the normal distribution, |
| 126 | * compute the limit on the definite integral from [-inf..z] that satisfies p, |
| 127 | * accurate to 16 decimal places. |
| 128 | * |
| 129 | * This implementation is based on: |
| 130 | * |
| 131 | * Wichura, M.J. (1988) Algorithm AS 241: The percentage points of the normal |
| 132 | * distribution. Applied Statistics 37(3):477-484. |
| 133 | */ |
| 134 | JEMALLOC_INLINE double |
| 135 | pt_norm(double p) |
| 136 | { |
| 137 | double q, r, ret; |
| 138 | |
| 139 | assert(p > 0.0 && p < 1.0); |
| 140 | |
| 141 | q = p - 0.5; |
| 142 | if (fabs(q) <= 0.425) { |
| 143 | /* p close to 1/2. */ |
| 144 | r = 0.180625 - q * q; |
| 145 | return (q * (((((((2.5090809287301226727e3 * r + |
| 146 | 3.3430575583588128105e4) * r + 6.7265770927008700853e4) * r |
| 147 | + 4.5921953931549871457e4) * r + 1.3731693765509461125e4) * |
| 148 | r + 1.9715909503065514427e3) * r + 1.3314166789178437745e2) |
| 149 | * r + 3.3871328727963666080e0) / |
| 150 | (((((((5.2264952788528545610e3 * r + |
| 151 | 2.8729085735721942674e4) * r + 3.9307895800092710610e4) * r |
| 152 | + 2.1213794301586595867e4) * r + 5.3941960214247511077e3) * |
| 153 | r + 6.8718700749205790830e2) * r + 4.2313330701600911252e1) |
| 154 | * r + 1.0)); |
| 155 | } else { |
| 156 | if (q < 0.0) |
| 157 | r = p; |
| 158 | else |
| 159 | r = 1.0 - p; |
| 160 | assert(r > 0.0); |
| 161 | |
| 162 | r = sqrt(-log(r)); |
| 163 | if (r <= 5.0) { |
| 164 | /* p neither close to 1/2 nor 0 or 1. */ |
| 165 | r -= 1.6; |
| 166 | ret = ((((((((7.74545014278341407640e-4 * r + |
| 167 | 2.27238449892691845833e-2) * r + |
| 168 | 2.41780725177450611770e-1) * r + |
| 169 | 1.27045825245236838258e0) * r + |
| 170 | 3.64784832476320460504e0) * r + |
| 171 | 5.76949722146069140550e0) * r + |
| 172 | 4.63033784615654529590e0) * r + |
| 173 | 1.42343711074968357734e0) / |
| 174 | (((((((1.05075007164441684324e-9 * r + |
| 175 | 5.47593808499534494600e-4) * r + |
| 176 | 1.51986665636164571966e-2) |
| 177 | * r + 1.48103976427480074590e-1) * r + |
| 178 | 6.89767334985100004550e-1) * r + |
| 179 | 1.67638483018380384940e0) * r + |
| 180 | 2.05319162663775882187e0) * r + 1.0)); |
| 181 | } else { |
| 182 | /* p near 0 or 1. */ |
| 183 | r -= 5.0; |
| 184 | ret = ((((((((2.01033439929228813265e-7 * r + |
| 185 | 2.71155556874348757815e-5) * r + |
| 186 | 1.24266094738807843860e-3) * r + |
| 187 | 2.65321895265761230930e-2) * r + |
| 188 | 2.96560571828504891230e-1) * r + |
| 189 | 1.78482653991729133580e0) * r + |
| 190 | 5.46378491116411436990e0) * r + |
| 191 | 6.65790464350110377720e0) / |
| 192 | (((((((2.04426310338993978564e-15 * r + |
| 193 | 1.42151175831644588870e-7) * r + |
| 194 | 1.84631831751005468180e-5) * r + |
| 195 | 7.86869131145613259100e-4) * r + |
| 196 | 1.48753612908506148525e-2) * r + |
| 197 | 1.36929880922735805310e-1) * r + |
| 198 | 5.99832206555887937690e-1) |
| 199 | * r + 1.0)); |
| 200 | } |
| 201 | if (q < 0.0) |
| 202 | ret = -ret; |
| 203 | return (ret); |
| 204 | } |
| 205 | } |
| 206 | |
| 207 | /* |
| 208 | * Given a value p in [0..1] of the lower tail area of the Chi^2 distribution |
| 209 | * with df degrees of freedom, where ln_gamma_df_2 is ln_gamma(df/2.0), compute |
| 210 | * the upper limit on the definite integral from [0..z] that satisfies p, |
| 211 | * accurate to 12 decimal places. |
| 212 | * |
| 213 | * This implementation is based on: |
| 214 | * |
| 215 | * Best, D.J., D.E. Roberts (1975) Algorithm AS 91: The percentage points of |
| 216 | * the Chi^2 distribution. Applied Statistics 24(3):385-388. |
| 217 | * |
| 218 | * Shea, B.L. (1991) Algorithm AS R85: A remark on AS 91: The percentage |
| 219 | * points of the Chi^2 distribution. Applied Statistics 40(1):233-235. |
| 220 | */ |
| 221 | JEMALLOC_INLINE double |
| 222 | pt_chi2(double p, double df, double ln_gamma_df_2) |
| 223 | { |
| 224 | double e, aa, xx, c, ch, a, q, p1, p2, t, x, b, s1, s2, s3, s4, s5, s6; |
| 225 | unsigned i; |
| 226 | |
| 227 | assert(p >= 0.0 && p < 1.0); |
| 228 | assert(df > 0.0); |
| 229 | |
| 230 | e = 5.0e-7; |
| 231 | aa = 0.6931471805; |
| 232 | |
| 233 | xx = 0.5 * df; |
| 234 | c = xx - 1.0; |
| 235 | |
| 236 | if (df < -1.24 * log(p)) { |
| 237 | /* Starting approximation for small Chi^2. */ |
| 238 | ch = pow(p * xx * exp(ln_gamma_df_2 + xx * aa), 1.0 / xx); |
| 239 | if (ch - e < 0.0) |
| 240 | return (ch); |
| 241 | } else { |
| 242 | if (df > 0.32) { |
| 243 | x = pt_norm(p); |
| 244 | /* |
| 245 | * Starting approximation using Wilson and Hilferty |
| 246 | * estimate. |
| 247 | */ |
| 248 | p1 = 0.222222 / df; |
| 249 | ch = df * pow(x * sqrt(p1) + 1.0 - p1, 3.0); |
| 250 | /* Starting approximation for p tending to 1. */ |
| 251 | if (ch > 2.2 * df + 6.0) { |
| 252 | ch = -2.0 * (log(1.0 - p) - c * log(0.5 * ch) + |
| 253 | ln_gamma_df_2); |
| 254 | } |
| 255 | } else { |
| 256 | ch = 0.4; |
| 257 | a = log(1.0 - p); |
| 258 | while (true) { |
| 259 | q = ch; |
| 260 | p1 = 1.0 + ch * (4.67 + ch); |
| 261 | p2 = ch * (6.73 + ch * (6.66 + ch)); |
| 262 | t = -0.5 + (4.67 + 2.0 * ch) / p1 - (6.73 + ch |
| 263 | * (13.32 + 3.0 * ch)) / p2; |
| 264 | ch -= (1.0 - exp(a + ln_gamma_df_2 + 0.5 * ch + |
| 265 | c * aa) * p2 / p1) / t; |
| 266 | if (fabs(q / ch - 1.0) - 0.01 <= 0.0) |
| 267 | break; |
| 268 | } |
| 269 | } |
| 270 | } |
| 271 | |
| 272 | for (i = 0; i < 20; i++) { |
| 273 | /* Calculation of seven-term Taylor series. */ |
| 274 | q = ch; |
| 275 | p1 = 0.5 * ch; |
| 276 | if (p1 < 0.0) |
| 277 | return (-1.0); |
| 278 | p2 = p - i_gamma(p1, xx, ln_gamma_df_2); |
| 279 | t = p2 * exp(xx * aa + ln_gamma_df_2 + p1 - c * log(ch)); |
| 280 | b = t / ch; |
| 281 | a = 0.5 * t - b * c; |
| 282 | s1 = (210.0 + a * (140.0 + a * (105.0 + a * (84.0 + a * (70.0 + |
| 283 | 60.0 * a))))) / 420.0; |
| 284 | s2 = (420.0 + a * (735.0 + a * (966.0 + a * (1141.0 + 1278.0 * |
| 285 | a)))) / 2520.0; |
| 286 | s3 = (210.0 + a * (462.0 + a * (707.0 + 932.0 * a))) / 2520.0; |
| 287 | s4 = (252.0 + a * (672.0 + 1182.0 * a) + c * (294.0 + a * |
| 288 | (889.0 + 1740.0 * a))) / 5040.0; |
| 289 | s5 = (84.0 + 264.0 * a + c * (175.0 + 606.0 * a)) / 2520.0; |
| 290 | s6 = (120.0 + c * (346.0 + 127.0 * c)) / 5040.0; |
| 291 | ch += t * (1.0 + 0.5 * t * s1 - b * c * (s1 - b * (s2 - b * (s3 |
| 292 | - b * (s4 - b * (s5 - b * s6)))))); |
| 293 | if (fabs(q / ch - 1.0) <= e) |
| 294 | break; |
| 295 | } |
| 296 | |
| 297 | return (ch); |
| 298 | } |
| 299 | |
| 300 | /* |
| 301 | * Given a value p in [0..1] and Gamma distribution shape and scale parameters, |
Jason Evans | e12eaf9 | 2014-12-08 14:40:14 -0800 | [diff] [blame^] | 302 | * compute the upper limit on the definite integral from [0..z] that satisfies |
Jason Evans | b1941c6 | 2013-12-09 23:36:37 -0800 | [diff] [blame] | 303 | * p. |
| 304 | */ |
| 305 | JEMALLOC_INLINE double |
| 306 | pt_gamma(double p, double shape, double scale, double ln_gamma_shape) |
| 307 | { |
| 308 | |
| 309 | return (pt_chi2(p, shape * 2.0, ln_gamma_shape) * 0.5 * scale); |
| 310 | } |
| 311 | #endif |