]>
git.proxmox.com Git - rustc.git/blob - src/jemalloc/test/include/test/math.h
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
);
9 #if (defined(JEMALLOC_ENABLE_INLINE) || defined(MATH_C_))
11 * Compute the natural log of Gamma(x), accurate to 10 decimal places.
13 * This implementation is based on:
15 * Pike, M.C., I.D. Hill (1966) Algorithm 291: Logarithm of Gamma function
16 * [S14]. Communications of the ACM 9(9):684.
18 JEMALLOC_INLINE
double
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
);
45 * Compute the incomplete Gamma ratio for [0..x], where p is the shape
46 * parameter, and ln_gamma_p is ln_gamma(p).
48 * This implementation is based on:
50 * Bhattacharjee, G.P. (1970) Algorithm AS 32: The incomplete Gamma integral.
51 * Applied Statistics 19:285-287.
53 JEMALLOC_INLINE
double
54 i_gamma(double x
, double p
, double ln_gamma_p
)
56 double acu
, factor
, oflo
, gin
, term
, rn
, a
, b
, an
, dif
;
69 factor
= exp(p
* log(x
) - x
- ln_gamma_p
);
71 if (x
<= 1.0 || x
< p
) {
72 /* Calculation by series expansion. */
87 /* Calculation by continued fraction. */
102 for (i
= 0; i
< 2; i
++)
103 pn
[i
+4] = b
* pn
[i
+2] - an
* pn
[i
];
106 dif
= fabs(gin
- rn
);
107 if (dif
<= acu
&& dif
<= acu
* rn
) {
108 gin
= 1.0 - factor
* gin
;
113 for (i
= 0; i
< 4; i
++)
116 if (fabs(pn
[4]) >= oflo
) {
117 for (i
= 0; i
< 4; i
++)
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.
129 * This implementation is based on:
131 * Wichura, M.J. (1988) Algorithm AS 241: The percentage points of the normal
132 * distribution. Applied Statistics 37(3):477-484.
134 JEMALLOC_INLINE
double
139 assert(p
> 0.0 && p
< 1.0);
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
)
164 /* p neither close to 1/2 nor 0 or 1. */
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));
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)
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.
213 * This implementation is based on:
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.
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.
221 JEMALLOC_INLINE
double
222 pt_chi2(double p
, double df
, double ln_gamma_df_2
)
224 double e
, aa
, xx
, c
, ch
, a
, q
, p1
, p2
, t
, x
, b
, s1
, s2
, s3
, s4
, s5
, s6
;
227 assert(p
>= 0.0 && p
< 1.0);
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
);
245 * Starting approximation using Wilson and Hilferty
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
) +
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)
272 for (i
= 0; i
< 20; i
++) {
273 /* Calculation of seven-term Taylor series. */
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
));
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 *
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
)
301 * Given a value p in [0..1] and Gamma distribution shape and scale parameters,
302 * compute the upper limit on the definite integral from [0..z] that satisfies
305 JEMALLOC_INLINE
double
306 pt_gamma(double p
, double shape
, double scale
, double ln_gamma_shape
)
309 return (pt_chi2(p
, shape
* 2.0, ln_gamma_shape
) * 0.5 * scale
);