statrs/function/
gamma.rs

1//! Provides the [gamma](https://en.wikipedia.org/wiki/Gamma_function) and
2//! related functions
3
4use crate::consts;
5use crate::error::StatsError;
6use crate::is_zero;
7use crate::prec;
8use crate::Result;
9use std::f64;
10
11/// Auxiliary variable when evaluating the `gamma_ln` function
12const GAMMA_R: f64 = 10.900511;
13
14/// Polynomial coefficients for approximating the `gamma_ln` function
15const GAMMA_DK: &[f64] = &[
16    2.48574089138753565546e-5,
17    1.05142378581721974210,
18    -3.45687097222016235469,
19    4.51227709466894823700,
20    -2.98285225323576655721,
21    1.05639711577126713077,
22    -1.95428773191645869583e-1,
23    1.70970543404441224307e-2,
24    -5.71926117404305781283e-4,
25    4.63399473359905636708e-6,
26    -2.71994908488607703910e-9,
27];
28
29/// Computes the logarithm of the gamma function
30/// with an accuracy of 16 floating point digits.
31/// The implementation is derived from
32/// "An Analysis of the Lanczos Gamma Approximation",
33/// Glendon Ralph Pugh, 2004 p. 116
34pub fn ln_gamma(x: f64) -> f64 {
35    if x < 0.5 {
36        let s = GAMMA_DK
37            .iter()
38            .enumerate()
39            .skip(1)
40            .fold(GAMMA_DK[0], |s, t| s + t.1 / (t.0 as f64 - x));
41
42        consts::LN_PI
43            - (f64::consts::PI * x).sin().ln()
44            - s.ln()
45            - consts::LN_2_SQRT_E_OVER_PI
46            - (0.5 - x) * ((0.5 - x + GAMMA_R) / f64::consts::E).ln()
47    } else {
48        let s = GAMMA_DK
49            .iter()
50            .enumerate()
51            .skip(1)
52            .fold(GAMMA_DK[0], |s, t| s + t.1 / (x + t.0 as f64 - 1.0));
53
54        s.ln()
55            + consts::LN_2_SQRT_E_OVER_PI
56            + (x - 0.5) * ((x - 0.5 + GAMMA_R) / f64::consts::E).ln()
57    }
58}
59
60/// Computes the gamma function with an accuracy
61/// of 16 floating point digits. The implementation
62/// is derived from "An Analysis of the Lanczos Gamma Approximation",
63/// Glendon Ralph Pugh, 2004 p. 116
64pub fn gamma(x: f64) -> f64 {
65    if x < 0.5 {
66        let s = GAMMA_DK
67            .iter()
68            .enumerate()
69            .skip(1)
70            .fold(GAMMA_DK[0], |s, t| s + t.1 / (t.0 as f64 - x));
71
72        f64::consts::PI
73            / ((f64::consts::PI * x).sin()
74                * s
75                * consts::TWO_SQRT_E_OVER_PI
76                * ((0.5 - x + GAMMA_R) / f64::consts::E).powf(0.5 - x))
77    } else {
78        let s = GAMMA_DK
79            .iter()
80            .enumerate()
81            .skip(1)
82            .fold(GAMMA_DK[0], |s, t| s + t.1 / (x + t.0 as f64 - 1.0));
83
84        s * consts::TWO_SQRT_E_OVER_PI * ((x - 0.5 + GAMMA_R) / f64::consts::E).powf(x - 0.5)
85    }
86}
87
88/// Computes the upper incomplete gamma function
89/// `Gamma(a,x) = int(exp(-t)t^(a-1), t=0..x) for a > 0, x > 0`
90/// where `a` is the argument for the gamma function and
91/// `x` is the lower intergral limit.
92///
93/// # Panics
94///
95/// if `a` or `x` are not in `(0, +inf)`
96pub fn gamma_ui(a: f64, x: f64) -> f64 {
97    checked_gamma_ui(a, x).unwrap()
98}
99
100/// Computes the upper incomplete gamma function
101/// `Gamma(a,x) = int(exp(-t)t^(a-1), t=0..x) for a > 0, x > 0`
102/// where `a` is the argument for the gamma function and
103/// `x` is the lower intergral limit.
104///
105/// # Errors
106///
107/// if `a` or `x` are not in `(0, +inf)`
108pub fn checked_gamma_ui(a: f64, x: f64) -> Result<f64> {
109    checked_gamma_ur(a, x).map(|x| x * gamma(a))
110}
111
112/// Computes the lower incomplete gamma function
113/// `gamma(a,x) = int(exp(-t)t^(a-1), t=0..x) for a > 0, x > 0`
114/// where `a` is the argument for the gamma function and `x`
115/// is the upper integral limit.
116///
117///
118/// # Panics
119///
120/// if `a` or `x` are not in `(0, +inf)`
121pub fn gamma_li(a: f64, x: f64) -> f64 {
122    checked_gamma_li(a, x).unwrap()
123}
124
125/// Computes the lower incomplete gamma function
126/// `gamma(a,x) = int(exp(-t)t^(a-1), t=0..x) for a > 0, x > 0`
127/// where `a` is the argument for the gamma function and `x`
128/// is the upper integral limit.
129///
130///
131/// # Errors
132///
133/// if `a` or `x` are not in `(0, +inf)`
134pub fn checked_gamma_li(a: f64, x: f64) -> Result<f64> {
135    checked_gamma_lr(a, x).map(|x| x * gamma(a))
136}
137
138/// Computes the upper incomplete regularized gamma function
139/// `Q(a,x) = 1 / Gamma(a) * int(exp(-t)t^(a-1), t=0..x) for a > 0, x > 0`
140/// where `a` is the argument for the gamma function and
141/// `x` is the lower integral limit.
142///
143/// # Remarks
144///
145/// Returns `f64::NAN` if either argument is `f64::NAN`
146///
147/// # Panics
148///
149/// if `a` or `x` are not in `(0, +inf)`
150pub fn gamma_ur(a: f64, x: f64) -> f64 {
151    checked_gamma_ur(a, x).unwrap()
152}
153
154/// Computes the upper incomplete regularized gamma function
155/// `Q(a,x) = 1 / Gamma(a) * int(exp(-t)t^(a-1), t=0..x) for a > 0, x > 0`
156/// where `a` is the argument for the gamma function and
157/// `x` is the lower integral limit.
158///
159/// # Remarks
160///
161/// Returns `f64::NAN` if either argument is `f64::NAN`
162///
163/// # Errors
164///
165/// if `a` or `x` are not in `(0, +inf)`
166pub fn checked_gamma_ur(a: f64, x: f64) -> Result<f64> {
167    if a.is_nan() || x.is_nan() {
168        return Ok(f64::NAN);
169    }
170    if a <= 0.0 || a == f64::INFINITY {
171        return Err(StatsError::ArgIntervalExcl("a", 0.0, f64::INFINITY));
172    }
173    if x <= 0.0 || x == f64::INFINITY {
174        return Err(StatsError::ArgIntervalExcl("x", 0.0, f64::INFINITY));
175    }
176
177    let eps = 0.000000000000001;
178    let big = 4503599627370496.0;
179    let big_inv = 2.22044604925031308085e-16;
180
181    if x < 1.0 || x <= a {
182        return Ok(1.0 - gamma_lr(a, x));
183    }
184
185    let mut ax = a * x.ln() - x - ln_gamma(a);
186    if ax < -709.78271289338399 {
187        return if a < x { Ok(0.0) } else { Ok(1.0) };
188    }
189
190    ax = ax.exp();
191    let mut y = 1.0 - a;
192    let mut z = x + y + 1.0;
193    let mut c = 0.0;
194    let mut pkm2 = 1.0;
195    let mut qkm2 = x;
196    let mut pkm1 = x + 1.0;
197    let mut qkm1 = z * x;
198    let mut ans = pkm1 / qkm1;
199    loop {
200        y += 1.0;
201        z += 2.0;
202        c += 1.0;
203        let yc = y * c;
204        let pk = pkm1 * z - pkm2 * yc;
205        let qk = qkm1 * z - qkm2 * yc;
206
207        pkm2 = pkm1;
208        pkm1 = pk;
209        qkm2 = qkm1;
210        qkm1 = qk;
211
212        if pk.abs() > big {
213            pkm2 *= big_inv;
214            pkm1 *= big_inv;
215            qkm2 *= big_inv;
216            qkm1 *= big_inv;
217        }
218
219        if !is_zero(qk) {
220            let r = pk / qk;
221            let t = ((ans - r) / r).abs();
222            ans = r;
223
224            if t <= eps {
225                break;
226            }
227        }
228    }
229    Ok(ans * ax)
230}
231
232/// Computes the lower incomplete regularized gamma function
233/// `P(a,x) = 1 / Gamma(a) * int(exp(-t)t^(a-1), t=0..x) for real a > 0, x > 0`
234/// where `a` is the argument for the gamma function and `x` is the upper
235/// integral limit.
236///
237/// # Remarks
238///
239/// Returns `f64::NAN` if either argument is `f64::NAN`
240///
241/// # Panics
242///
243/// if `a` or `x` are not in `(0, +inf)`
244pub fn gamma_lr(a: f64, x: f64) -> f64 {
245    checked_gamma_lr(a, x).unwrap()
246}
247
248/// Computes the lower incomplete regularized gamma function
249/// `P(a,x) = 1 / Gamma(a) * int(exp(-t)t^(a-1), t=0..x) for real a > 0, x > 0`
250/// where `a` is the argument for the gamma function and `x` is the upper
251/// integral limit.
252///
253/// # Remarks
254///
255/// Returns `f64::NAN` if either argument is `f64::NAN`
256///
257/// # Errors
258///
259/// if `a` or `x` are not in `(0, +inf)`
260pub fn checked_gamma_lr(a: f64, x: f64) -> Result<f64> {
261    if a.is_nan() || x.is_nan() {
262        return Ok(f64::NAN);
263    }
264    if a <= 0.0 || a == f64::INFINITY {
265        return Err(StatsError::ArgIntervalExcl("a", 0.0, f64::INFINITY));
266    }
267    if x <= 0.0 || x == f64::INFINITY {
268        return Err(StatsError::ArgIntervalExcl("x", 0.0, f64::INFINITY));
269    }
270
271    let eps = 0.000000000000001;
272    let big = 4503599627370496.0;
273    let big_inv = 2.22044604925031308085e-16;
274
275    if prec::almost_eq(a, 0.0, prec::DEFAULT_F64_ACC) {
276        return Ok(1.0);
277    }
278    if prec::almost_eq(x, 0.0, prec::DEFAULT_F64_ACC) {
279        return Ok(0.0);
280    }
281
282    let ax = a * x.ln() - x - ln_gamma(a);
283    if ax < -709.78271289338399 {
284        if a < x {
285            return Ok(1.0);
286        }
287        return Ok(0.0);
288    }
289    if x <= 1.0 || x <= a {
290        let mut r2 = a;
291        let mut c2 = 1.0;
292        let mut ans2 = 1.0;
293        loop {
294            r2 += 1.0;
295            c2 *= x / r2;
296            ans2 += c2;
297
298            if c2 / ans2 <= eps {
299                break;
300            }
301        }
302        return Ok(ax.exp() * ans2 / a);
303    }
304
305    let mut y = 1.0 - a;
306    let mut z = x + y + 1.0;
307    let mut c = 0;
308
309    let mut p3 = 1.0;
310    let mut q3 = x;
311    let mut p2 = x + 1.0;
312    let mut q2 = z * x;
313    let mut ans = p2 / q2;
314
315    loop {
316        y += 1.0;
317        z += 2.0;
318        c += 1;
319        let yc = y * f64::from(c);
320
321        let p = p2 * z - p3 * yc;
322        let q = q2 * z - q3 * yc;
323
324        p3 = p2;
325        p2 = p;
326        q3 = q2;
327        q2 = q;
328
329        if p.abs() > big {
330            p3 *= big_inv;
331            p2 *= big_inv;
332            q3 *= big_inv;
333            q2 *= big_inv;
334        }
335
336        if q != 0.0 {
337            let nextans = p / q;
338            let error = ((ans - nextans) / nextans).abs();
339            ans = nextans;
340
341            if error <= eps {
342                break;
343            }
344        }
345    }
346    Ok(1.0 - ax.exp() * ans)
347}
348
349/// Computes the Digamma function which is defined as the derivative of
350/// the log of the gamma function. The implementation is based on
351/// "Algorithm AS 103", Jose Bernardo, Applied Statistics, Volume 25, Number 3
352/// 1976, pages 315 - 317
353pub fn digamma(x: f64) -> f64 {
354    let c = 12.0;
355    let d1 = -0.57721566490153286;
356    let d2 = 1.6449340668482264365;
357    let s = 1e-6;
358    let s3 = 1.0 / 12.0;
359    let s4 = 1.0 / 120.0;
360    let s5 = 1.0 / 252.0;
361    let s6 = 1.0 / 240.0;
362    let s7 = 1.0 / 132.0;
363
364    if x == f64::NEG_INFINITY || x.is_nan() {
365        return f64::NAN;
366    }
367    if x <= 0.0 && ulps_eq!(x.floor(), x) {
368        return f64::NEG_INFINITY;
369    }
370    if x < 0.0 {
371        return digamma(1.0 - x) + f64::consts::PI / (-f64::consts::PI * x).tan();
372    }
373    if x <= s {
374        return d1 - 1.0 / x + d2 * x;
375    }
376
377    let mut result = 0.0;
378    let mut z = x;
379    while z < c {
380        result -= 1.0 / z;
381        z += 1.0;
382    }
383
384    if z >= c {
385        let mut r = 1.0 / z;
386        result += z.ln() - 0.5 * r;
387        r *= r;
388
389        result -= r * (s3 - r * (s4 - r * (s5 - r * (s6 - r * s7))));
390    }
391    result
392}
393
394pub fn inv_digamma(x: f64) -> f64 {
395    if x.is_nan() {
396        return f64::NAN;
397    }
398    if x == f64::NEG_INFINITY {
399        return 0.0;
400    }
401    if x == f64::INFINITY {
402        return f64::INFINITY;
403    }
404    let mut y = x.exp();
405    let mut i = 1.0;
406    while i > 1e-15 {
407        y += i * signum(x - digamma(y));
408        i /= 2.0;
409    }
410    y
411}
412
413// modified signum that returns 0.0 if x == 0.0. Used
414// by inv_digamma, may consider extracting into a public
415// method
416fn signum(x: f64) -> f64 {
417    if x == 0.0 {
418        0.0
419    } else {
420        x.signum()
421    }
422}
423
424/// Computes the [multivariate gamma function](https://en.wikipedia.org/wiki/Multivariate_gamma_function).
425pub fn mvgamma(p: i64, a: f64) -> f64 {
426    let mut res = std::f64::consts::PI.powf((p * (p - 1)) as f64 / 4.0);
427    for ii in 1..=p {
428        res *= gamma(a + (1 - ii) as f64 / 2.0);
429    }
430    res
431}
432
433/// Computes the log [multivariate gamma function](https://en.wikipedia.org/wiki/Multivariate_gamma_function).
434pub fn mvlgamma(p: i64, a: f64) -> f64 {
435    let mut res = (p * (p - 1)) as f64 * consts::LN_PI / 4.0;
436    for ii in 1..=p {
437        res += ln_gamma(a + (1 - ii) as f64 / 2.0);
438    }
439    res
440}
441
442#[rustfmt::skip]
443#[cfg(test)]
444mod tests {
445    use std::f64::{self, consts};
446
447    #[test]
448    fn test_gamma() {
449        assert!(super::gamma(f64::NAN).is_nan());
450        assert_almost_eq!(super::gamma(1.000001e-35), 9.9999900000099999900000099999899999522784235098567139293e+34, 1e20);
451        assert_almost_eq!(super::gamma(1.000001e-10), 9.99998999943278432519738283781280989934496494539074049002e+9, 1e-5);
452        assert_almost_eq!(super::gamma(1.000001e-5), 99999.32279432557746387178953902739303931424932435387031653234, 1e-10);
453        assert_almost_eq!(super::gamma(1.000001e-2), 99.43248512896257405886134437203369035261893114349805309870831, 1e-13);
454        assert_almost_eq!(super::gamma(-4.8), -0.06242336135475955314181664931547009890495158793105543559676, 1e-13);
455        assert_almost_eq!(super::gamma(-1.5), 2.363271801207354703064223311121526910396732608163182837618410, 1e-13);
456        assert_almost_eq!(super::gamma(-0.5), -3.54490770181103205459633496668229036559509891224477425642761, 1e-13);
457        assert_almost_eq!(super::gamma(1.0e-5 + 1.0e-16), 99999.42279322556767360213300482199406241771308740302819426480, 1e-9);
458        assert_almost_eq!(super::gamma(0.1), 9.513507698668731836292487177265402192550578626088377343050000, 1e-14);
459        assert_eq!(super::gamma(1.0 - 1.0e-14), 1.000000000000005772156649015427511664653698987042926067639529);
460        assert_almost_eq!(super::gamma(1.0), 1.0, 1e-15);
461        assert_almost_eq!(super::gamma(1.0 + 1.0e-14), 0.99999999999999422784335098477029953441189552403615306268023, 1e-15);
462        assert_almost_eq!(super::gamma(1.5), 0.886226925452758013649083741670572591398774728061193564106903, 1e-14);
463        assert_almost_eq!(super::gamma(consts::PI/2.0), 0.890560890381539328010659635359121005933541962884758999762766, 1e-15);
464        assert_eq!(super::gamma(2.0), 1.0);
465        assert_almost_eq!(super::gamma(2.5), 1.329340388179137020473625612505858887098162092091790346160355, 1e-13);
466        assert_almost_eq!(super::gamma(3.0), 2.0, 1e-14);
467        assert_almost_eq!(super::gamma(consts::PI), 2.288037795340032417959588909060233922889688153356222441199380, 1e-13);
468        assert_almost_eq!(super::gamma(3.5), 3.323350970447842551184064031264647217745405230229475865400889, 1e-14);
469        assert_almost_eq!(super::gamma(4.0), 6.0, 1e-13);
470        assert_almost_eq!(super::gamma(4.5), 11.63172839656744892914422410942626526210891830580316552890311, 1e-12);
471        assert_almost_eq!(super::gamma(5.0 - 1.0e-14), 23.99999999999963853175957637087420162718107213574617032780374, 1e-13);
472        assert_almost_eq!(super::gamma(5.0), 24.0, 1e-12);
473        assert_almost_eq!(super::gamma(5.0 + 1.0e-14), 24.00000000000036146824042363510111050137786752408660789873592, 1e-12);
474        assert_almost_eq!(super::gamma(5.5), 52.34277778455352018114900849241819367949013237611424488006401, 1e-12);
475        assert_almost_eq!(super::gamma(10.1), 454760.7514415859508673358368319076190405047458218916492282448, 1e-7);
476        assert_almost_eq!(super::gamma(150.0 + 1.0e-12), 3.8089226376496421386707466577615064443807882167327097140e+260, 1e248);
477    }
478
479    #[test]
480    fn test_ln_gamma() {
481        assert!(super::ln_gamma(f64::NAN).is_nan());
482        assert_eq!(super::ln_gamma(1.000001e-35), 80.59047725479209894029636783061921392709972287131139201585211);
483        assert_almost_eq!(super::ln_gamma(1.000001e-10), 23.02584992988323521564308637407936081168344192865285883337793, 1e-14);
484        assert_almost_eq!(super::ln_gamma(1.000001e-5), 11.51291869289055371493077240324332039045238086972508869965363, 1e-14);
485        assert_eq!(super::ln_gamma(1.000001e-2), 4.599478872433667224554543378460164306444416156144779542513592);
486        assert_almost_eq!(super::ln_gamma(0.1), 2.252712651734205959869701646368495118615627222294953765041739, 1e-14);
487        assert_almost_eq!(super::ln_gamma(1.0 - 1.0e-14), 5.772156649015410852768463312546533565566459794933360600e-15, 1e-15);
488        assert_almost_eq!(super::ln_gamma(1.0), 0.0, 1e-15);
489        assert_almost_eq!(super::ln_gamma(1.0 + 1.0e-14), -5.77215664901524635936177848990288632404978978079827014e-15, 1e-15);
490        assert_almost_eq!(super::ln_gamma(1.5), -0.12078223763524522234551844578164721225185272790259946836386, 1e-14);
491        assert_almost_eq!(super::ln_gamma(consts::PI/2.0), -0.11590380084550241329912089415904874214542604767006895, 1e-14);
492        assert_eq!(super::ln_gamma(2.0), 0.0);
493        assert_almost_eq!(super::ln_gamma(2.5), 0.284682870472919159632494669682701924320137695559894729250145, 1e-13);
494        assert_almost_eq!(super::ln_gamma(3.0), 0.693147180559945309417232121458176568075500134360255254120680, 1e-14);
495        assert_almost_eq!(super::ln_gamma(consts::PI), 0.82769459232343710152957855845235995115350173412073715, 1e-13);
496        assert_almost_eq!(super::ln_gamma(3.5), 1.200973602347074224816021881450712995770238915468157197042113, 1e-14);
497        assert_almost_eq!(super::ln_gamma(4.0), 1.791759469228055000812477358380702272722990692183004705855374, 1e-14);
498        assert_almost_eq!(super::ln_gamma(4.5), 2.453736570842442220504142503435716157331823510689763131380823, 1e-13);
499        assert_almost_eq!(super::ln_gamma(5.0 - 1.0e-14), 3.178053830347930558470257283303394288448414225994179545985931, 1e-14);
500        assert_almost_eq!(super::ln_gamma(5.0), 3.178053830347945619646941601297055408873990960903515214096734, 1e-14);
501        assert_almost_eq!(super::ln_gamma(5.0 + 1.0e-14), 3.178053830347960680823625919312848824873279228348981287761046, 1e-13);
502        assert_almost_eq!(super::ln_gamma(5.5), 3.957813967618716293877400855822590998551304491975006780729532, 1e-14);
503        assert_almost_eq!(super::ln_gamma(10.1), 13.02752673863323795851370097886835481188051062306253294740504, 1e-14);
504        assert_almost_eq!(super::ln_gamma(150.0 + 1.0e-12), 600.0094705553324354062157737572509902987070089159051628001813, 1e-12);
505        assert_almost_eq!(super::ln_gamma(1.001e+7), 1.51342135323817913130119829455205139905331697084416059779e+8, 1e-13);
506    }
507
508    #[test]
509    fn test_gamma_lr() {
510        assert!(super::gamma_lr(f64::NAN, f64::NAN).is_nan());
511        assert_almost_eq!(super::gamma_lr(0.1, 1.0), 0.97587265627367222115949155252812057714751052498477013, 1e-14);
512        assert_eq!(super::gamma_lr(0.1, 2.0), 0.99432617602018847196075251078067514034772764693462125);
513        assert_eq!(super::gamma_lr(0.1, 8.0), 0.99999507519205198048686442150578226823401842046310854);
514        assert_almost_eq!(super::gamma_lr(1.5, 1.0), 0.42759329552912016600095238564127189392715996802703368, 1e-13);
515        assert_almost_eq!(super::gamma_lr(1.5, 2.0), 0.73853587005088937779717792402407879809718939080920993, 1e-15);
516        assert_eq!(super::gamma_lr(1.5, 8.0), 0.99886601571021467734329986257903021041757398191304284);
517        assert_almost_eq!(super::gamma_lr(2.5, 1.0), 0.15085496391539036377410688601371365034788861473418704, 1e-13);
518        assert_almost_eq!(super::gamma_lr(2.5, 2.0), 0.45058404864721976739416885516693969548484517509263197, 1e-14);
519        assert_almost_eq!(super::gamma_lr(2.5, 8.0), 0.99315592607757956900093935107222761316136944145439676, 1e-15);
520        assert_almost_eq!(super::gamma_lr(5.5, 1.0), 0.0015041182825838038421585211353488839717739161316985392, 1e-15);
521        assert_almost_eq!(super::gamma_lr(5.5, 2.0), 0.030082976121226050615171484772387355162056796585883967, 1e-14);
522        assert_almost_eq!(super::gamma_lr(5.5, 8.0), 0.85886911973294184646060071855669224657735916933487681, 1e-14);
523        assert_almost_eq!(super::gamma_lr(100.0, 0.5), 0.0, 1e-188);
524        assert_almost_eq!(super::gamma_lr(100.0, 1.5), 0.0, 1e-141);
525        assert_almost_eq!(super::gamma_lr(100.0, 90.0), 0.1582209891864301681049696996709105316998233457433473, 1e-13);
526        assert_almost_eq!(super::gamma_lr(100.0, 100.0), 0.5132987982791486648573142565640291634709251499279450, 1e-13);
527        assert_almost_eq!(super::gamma_lr(100.0, 110.0), 0.8417213299399129061982996209829688531933500308658222, 1e-13);
528        assert_almost_eq!(super::gamma_lr(100.0, 200.0), 1.0, 1e-14);
529        assert_eq!(super::gamma_lr(500.0, 0.5), 0.0);
530        assert_eq!(super::gamma_lr(500.0, 1.5), 0.0);
531        assert_almost_eq!(super::gamma_lr(500.0, 200.0), 0.0, 1e-70);
532        assert_almost_eq!(super::gamma_lr(500.0, 450.0), 0.0107172380912897415573958770655204965434869949241480, 1e-14);
533        assert_almost_eq!(super::gamma_lr(500.0, 500.0), 0.5059471461707603580470479574412058032802735425634263, 1e-13);
534        assert_almost_eq!(super::gamma_lr(500.0, 550.0), 0.9853855918737048059548470006900844665580616318702748, 1e-14);
535        assert_almost_eq!(super::gamma_lr(500.0, 700.0), 1.0, 1e-15);
536        assert_eq!(super::gamma_lr(1000.0, 10000.0), 1.0);
537        assert_eq!(super::gamma_lr(1e+50, 1e+48), 0.0);
538        assert_eq!(super::gamma_lr(1e+50, 1e+52), 1.0);
539    }
540
541    #[test]
542    #[should_panic]
543    fn test_gamma_lr_a_lower_bound() {
544        super::gamma_lr(-1.0, 1.0);
545    }
546
547    #[test]
548    #[should_panic]
549    fn test_gamma_lr_a_upper_bound() {
550        super::gamma_lr(f64::INFINITY, 1.0);
551    }
552
553    #[test]
554    #[should_panic]
555    fn test_gamma_lr_x_lower_bound() {
556        super::gamma_lr(1.0, -1.0);
557    }
558
559    #[test]
560    #[should_panic]
561    fn test_gamma_lr_x_upper_bound() {
562        super::gamma_lr(1.0, f64::INFINITY);
563    }
564
565    #[test]
566    fn test_checked_gamma_lr_a_lower_bound() {
567        assert!(super::checked_gamma_lr(-1.0, 1.0).is_err());
568    }
569
570    #[test]
571    fn test_checked_gamma_lr_a_upper_bound() {
572        assert!(super::checked_gamma_lr(f64::INFINITY, 1.0).is_err());
573    }
574
575    #[test]
576    fn test_checked_gamma_lr_x_lower_bound() {
577        assert!(super::checked_gamma_lr(1.0, -1.0).is_err());
578    }
579
580    #[test]
581    fn test_checked_gamma_lr_x_upper_bound() {
582        assert!(super::checked_gamma_lr(1.0, f64::INFINITY).is_err());
583    }
584
585    #[test]
586    fn test_gamma_li() {
587        assert!(super::gamma_li(f64::NAN, f64::NAN).is_nan());
588        assert_almost_eq!(super::gamma_li(0.1, 1.0), 9.2839720283798852469443229940217320532607158711056334, 1e-14);
589        assert_almost_eq!(super::gamma_li(0.1, 2.0), 9.4595297305559030536119885480983751098528458886962883, 1e-14);
590        assert_almost_eq!(super::gamma_li(0.1, 8.0), 9.5134608464704033372127589212547718314010339263844976, 1e-13);
591        assert_almost_eq!(super::gamma_li(1.5, 1.0), 0.37894469164098470380394366597039213790868855578083847, 1e-15);
592        assert_almost_eq!(super::gamma_li(1.5, 2.0), 0.65451037345177732033319477475056262302270310457635612, 1e-14);
593        assert_almost_eq!(super::gamma_li(1.5, 8.0), 0.88522195804210983776635107858848816480298923071075222, 1e-13);
594        assert_almost_eq!(super::gamma_li(2.5, 1.0), 0.20053759629003473411039172879412733941722170263949, 1e-16);
595        assert_almost_eq!(super::gamma_li(2.5, 2.0), 0.59897957413602228465664030130712917348327070206302442, 1e-15);
596        assert_almost_eq!(super::gamma_li(2.5, 8.0), 1.3202422842943799358198434659248530581833764879301293, 1e-14);
597        assert_almost_eq!(super::gamma_li(5.5, 1.0), 0.078729729026968321691794205337720556329618007004848672, 1e-16);
598        assert_almost_eq!(super::gamma_li(5.5, 2.0), 1.5746265342113649473739798668921124454837064926448459, 1e-15);
599        assert_almost_eq!(super::gamma_li(5.5, 8.0), 44.955595480196465884619737757794960132425035578313584, 1e-12);
600    }
601
602    #[test]
603    #[should_panic]
604    fn test_gamma_li_a_lower_bound() {
605        super::gamma_li(-1.0, 1.0);
606    }
607
608    #[test]
609    #[should_panic]
610    fn test_gamma_li_a_upper_bound() {
611        super::gamma_li(f64::INFINITY, 1.0);
612    }
613
614    #[test]
615    #[should_panic]
616    fn test_gamma_li_x_lower_bound() {
617        super::gamma_li(1.0, -1.0);
618    }
619
620    #[test]
621    #[should_panic]
622    fn test_gamma_li_x_upper_bound() {
623        super::gamma_li(1.0, f64::INFINITY);
624    }
625
626    #[test]
627    fn test_checked_gamma_li_a_lower_bound() {
628        assert!(super::checked_gamma_li(-1.0, 1.0).is_err());
629    }
630
631    #[test]
632    fn test_checked_gamma_li_a_upper_bound() {
633        assert!(super::checked_gamma_li(f64::INFINITY, 1.0).is_err());
634    }
635
636    #[test]
637    fn test_checked_gamma_li_x_lower_bound() {
638        assert!(super::checked_gamma_li(1.0, -1.0).is_err());
639    }
640
641    #[test]
642    fn test_checked_gamma_li_x_upper_bound() {
643        assert!(super::checked_gamma_li(1.0, f64::INFINITY).is_err());
644    }
645
646    // TODO: precision testing could be more accurate, borrowed wholesale from Math.NET
647    #[test]
648    fn test_gamma_ur() {
649        assert!(super::gamma_ur(f64::NAN, f64::NAN).is_nan());
650        assert_almost_eq!(super::gamma_ur(0.1, 1.0), 0.0241273437263277773829694356333550393309597428392044, 1e-13);
651        assert_almost_eq!(super::gamma_ur(0.1, 2.0), 0.0056738239798115280392474892193248596522723530653781, 1e-13);
652        assert_almost_eq!(super::gamma_ur(0.1, 8.0), 0.0000049248079480195131355784942177317659815795368919702, 1e-13);
653        assert_almost_eq!(super::gamma_ur(1.5, 1.0), 0.57240670447087983399904761435872810607284003197297, 1e-13);
654        assert_almost_eq!(super::gamma_ur(1.5, 2.0), 0.26146412994911062220282207597592120190281060919079, 1e-13);
655        assert_almost_eq!(super::gamma_ur(1.5, 8.0), 0.0011339842897853226567001374209697895824260180869567, 1e-13);
656        assert_almost_eq!(super::gamma_ur(2.5, 1.0), 0.84914503608460963622589311398628634965211138526581, 1e-13);
657        assert_almost_eq!(super::gamma_ur(2.5, 2.0), 0.54941595135278023260583114483306030451515482490737, 1e-13);
658        assert_almost_eq!(super::gamma_ur(2.5, 8.0), 0.0068440739224204309990606489277723868386305585456026, 1e-13);
659        assert_almost_eq!(super::gamma_ur(5.5, 1.0), 0.9984958817174161961578414788646511160282260838683, 1e-13);
660        assert_almost_eq!(super::gamma_ur(5.5, 2.0), 0.96991702387877394938482851522761264483794320341412, 1e-13);
661        assert_almost_eq!(super::gamma_ur(5.5, 8.0), 0.14113088026705815353939928144330775342264083066512, 1e-13);
662        assert_almost_eq!(super::gamma_ur(100.0, 0.5), 1.0, 1e-14);
663        assert_almost_eq!(super::gamma_ur(100.0, 1.5), 1.0, 1e-14);
664        assert_almost_eq!(super::gamma_ur(100.0, 90.0), 0.8417790108135698318950303003290894683001766542566526, 1e-12);
665        assert_almost_eq!(super::gamma_ur(100.0, 100.0), 0.4867012017208513351426857434359708365290748500720549, 1e-12);
666        assert_almost_eq!(super::gamma_ur(100.0, 110.0), 0.1582786700600870938017003790170311468066499691341777, 1e-12);
667        assert_almost_eq!(super::gamma_ur(100.0, 200.0), 0.0, 1e-14);
668        assert_almost_eq!(super::gamma_ur(500.0, 0.5), 1.0, 1e-14);
669        assert_almost_eq!(super::gamma_ur(500.0, 1.5), 1.0, 1e-14);
670        assert_almost_eq!(super::gamma_ur(500.0, 200.0), 1.0, 1e-14);
671        assert_almost_eq!(super::gamma_ur(500.0, 450.0), 0.9892827619087102584426041229344795034565130050758519, 1e-12);
672        assert_almost_eq!(super::gamma_ur(500.0, 500.0), 0.4940528538292396419529520425587941967197264574365736, 1e-12);
673        assert_almost_eq!(super::gamma_ur(500.0, 550.0), 0.0146144081262951940451529993099155334419383681297251, 1e-12);
674        assert_almost_eq!(super::gamma_ur(500.0, 700.0), 0.0, 1e-14);
675        assert_almost_eq!(super::gamma_ur(1000.0, 10000.0), 0.0, 1e-14);
676        assert_almost_eq!(super::gamma_ur(1e+50, 1e+48), 1.0, 1e-14);
677        assert_almost_eq!(super::gamma_ur(1e+50, 1e+52), 0.0, 1e-14);
678    }
679
680    #[test]
681    #[should_panic]
682    fn test_gamma_ur_a_lower_bound() {
683        super::gamma_ur(-1.0, 1.0);
684    }
685
686    #[test]
687    #[should_panic]
688    fn test_gamma_ur_a_upper_bound() {
689        super::gamma_ur(f64::INFINITY, 1.0);
690    }
691
692    #[test]
693    #[should_panic]
694    fn test_gamma_ur_x_lower_bound() {
695        super::gamma_ur(1.0, -1.0);
696    }
697
698    #[test]
699    #[should_panic]
700    fn test_gamma_ur_x_upper_bound() {
701        super::gamma_ur(1.0, f64::INFINITY);
702    }
703
704    #[test]
705    fn test_checked_gamma_ur_a_lower_bound() {
706        assert!(super::checked_gamma_ur(-1.0, 1.0).is_err());
707    }
708
709    #[test]
710    fn test_checked_gamma_ur_a_upper_bound() {
711        assert!(super::checked_gamma_ur(f64::INFINITY, 1.0).is_err());
712    }
713
714    #[test]
715    fn test_checked_gamma_ur_x_lower_bound() {
716        assert!(super::checked_gamma_ur(1.0, -1.0).is_err());
717    }
718
719    #[test]
720    fn test_checked_gamma_ur_x_upper_bound() {
721        assert!(super::checked_gamma_ur(1.0, f64::INFINITY).is_err());
722    }
723
724    #[test]
725    fn test_gamma_ui() {
726        assert!(super::gamma_ui(f64::NAN, f64::NAN).is_nan());
727        assert_almost_eq!(super::gamma_ui(0.1, 1.0), 0.2295356702888460382790772147651768201739736396141314, 1e-14);
728        assert_almost_eq!(super::gamma_ui(0.1, 2.0), 0.053977968112828232195991347726857391060870217694027, 1e-15);
729        assert_almost_eq!(super::gamma_ui(0.1, 8.0), 0.000046852198327948595220974570460669512682180005810156, 1e-19);
730        assert_almost_eq!(super::gamma_ui(1.5, 1.0), 0.50728223381177330984514007570018045349008617228036, 1e-14);
731        assert_almost_eq!(super::gamma_ui(1.5, 2.0), 0.23171655200098069331588896692000996837607162348484, 1e-15);
732        assert_almost_eq!(super::gamma_ui(1.5, 8.0), 0.0010049674106481758827326630820844265957854973504417, 1e-17);
733        assert_almost_eq!(super::gamma_ui(2.5, 1.0), 1.1288027918891022863632338837117315476809403894523, 1e-14);
734        assert_almost_eq!(super::gamma_ui(2.5, 2.0), 0.73036081404311473581698531119872971361489139002877, 1e-14);
735        assert_almost_eq!(super::gamma_ui(2.5, 8.0), 0.0090981038847570846537821465810058289147856041616617, 1e-17);
736        assert_almost_eq!(super::gamma_ui(5.5, 1.0), 52.264048055526551859457214287080473123160514369109, 1e-12);
737        assert_almost_eq!(super::gamma_ui(5.5, 2.0), 50.768151250342155233775028625526081234006425883469, 1e-12);
738        assert_almost_eq!(super::gamma_ui(5.5, 8.0), 7.3871823043570542965292707346232335470650967978006, 1e-13);
739    }
740
741    #[test]
742    #[should_panic]
743    fn test_gamma_ui_a_lower_bound() {
744        super::gamma_ui(-1.0, 1.0);
745    }
746
747    #[test]
748    #[should_panic]
749    fn test_gamma_ui_a_upper_bound() {
750        super::gamma_ui(f64::INFINITY, 1.0);
751    }
752
753    #[test]
754    #[should_panic]
755    fn test_gamma_ui_x_lower_bound() {
756        super::gamma_ui(1.0, -1.0);
757    }
758
759    #[test]
760    #[should_panic]
761    fn test_gamma_ui_x_upper_bound() {
762        super::gamma_ui(1.0, f64::INFINITY);
763    }
764
765    #[test]
766    fn test_checked_gamma_ui_a_lower_bound() {
767        assert!(super::checked_gamma_ui(-1.0, 1.0).is_err());
768    }
769
770    #[test]
771    fn test_checked_gamma_ui_a_upper_bound() {
772        assert!(super::checked_gamma_ui(f64::INFINITY, 1.0).is_err());
773    }
774
775    #[test]
776    fn test_checked_gamma_ui_x_lower_bound() {
777        assert!(super::checked_gamma_ui(1.0, -1.0).is_err());
778    }
779
780    #[test]
781    fn test_checked_gamma_ui_x_upper_bound() {
782        assert!(super::checked_gamma_ui(1.0, f64::INFINITY).is_err());
783    }
784
785    // TODO: precision testing could be more accurate
786    #[test]
787    fn test_digamma() {
788        assert!(super::digamma(f64::NAN).is_nan());
789        assert_almost_eq!(super::digamma(-1.5), 0.70315664064524318722569033366791109947350706200623256, 1e-14);
790        assert_almost_eq!(super::digamma(-0.5), 0.036489973978576520559023667001244432806840395339565891, 1e-14);
791        assert_almost_eq!(super::digamma(0.1), -10.423754940411076232100295314502760886768558023951363, 1e-14);
792        assert_almost_eq!(super::digamma(1.0), -0.57721566490153286060651209008240243104215933593992359, 1e-14);
793        assert_almost_eq!(super::digamma(1.5), 0.036489973978576520559023667001244432806840395339565888, 1e-14);
794        assert_almost_eq!(super::digamma(consts::PI / 2.0), 0.10067337642740238636795561404029690452798358068944001, 1e-14);
795        assert_almost_eq!(super::digamma(2.0), 0.42278433509846713939348790991759756895784066406007641, 1e-14);
796        assert_almost_eq!(super::digamma(2.5), 0.70315664064524318722569033366791109947350706200623255, 1e-14);
797        assert_almost_eq!(super::digamma(3.0), 0.92278433509846713939348790991759756895784066406007641, 1e-14);
798        assert_almost_eq!(super::digamma(consts::PI), 0.97721330794200673329206948640618234364083460999432603, 1e-14);
799        assert_almost_eq!(super::digamma(3.5), 1.1031566406452431872256903336679110994735070620062326, 1e-14);
800        assert_almost_eq!(super::digamma(4.0), 1.2561176684318004727268212432509309022911739973934097, 1e-14);
801        assert_almost_eq!(super::digamma(4.5), 1.3888709263595289015114046193821968137592213477205183, 1e-14);
802        assert_almost_eq!(super::digamma(5.0), 1.5061176684318004727268212432509309022911739973934097, 1e-14);
803        assert_almost_eq!(super::digamma(5.5), 1.6110931485817511237336268416044190359814435699427405, 1e-14);
804        assert_almost_eq!(super::digamma(10.1), 2.2622143570941481235561593642219403924532310597356171, 1e-14);
805    }
806
807    #[test]
808    fn test_inv_digamma() {
809        assert!(super::inv_digamma(f64::NAN).is_nan());
810        assert_eq!(super::inv_digamma(f64::NEG_INFINITY), 0.0);
811        assert_almost_eq!(super::inv_digamma(-10.423754940411076232100295314502760886768558023951363), 0.1, 1e-15);
812        assert_almost_eq!(super::inv_digamma(-0.57721566490153286060651209008240243104215933593992359), 1.0, 1e-14);
813        assert_almost_eq!(super::inv_digamma(0.036489973978576520559023667001244432806840395339565888), 1.5, 1e-14);
814        assert_almost_eq!(super::inv_digamma(0.10067337642740238636795561404029690452798358068944001), consts::PI / 2.0, 1e-14);
815        assert_almost_eq!(super::inv_digamma(0.42278433509846713939348790991759756895784066406007641), 2.0, 1e-14);
816        assert_almost_eq!(super::inv_digamma(0.70315664064524318722569033366791109947350706200623255), 2.5, 1e-14);
817        assert_almost_eq!(super::inv_digamma(0.92278433509846713939348790991759756895784066406007641), 3.0, 1e-14);
818        assert_almost_eq!(super::inv_digamma(0.97721330794200673329206948640618234364083460999432603), consts::PI, 1e-14);
819        assert_almost_eq!(super::inv_digamma(1.1031566406452431872256903336679110994735070620062326), 3.5, 1e-14);
820        assert_almost_eq!(super::inv_digamma(1.2561176684318004727268212432509309022911739973934097), 4.0, 1e-14);
821        assert_almost_eq!(super::inv_digamma(1.3888709263595289015114046193821968137592213477205183), 4.5, 1e-14);
822        assert_almost_eq!(super::inv_digamma(1.5061176684318004727268212432509309022911739973934097), 5.0, 1e-14);
823        assert_almost_eq!(super::inv_digamma(1.6110931485817511237336268416044190359814435699427405), 5.5, 1e-14);
824        assert_almost_eq!(super::inv_digamma(2.2622143570941481235561593642219403924532310597356171), 10.1, 1e-13);
825    }
826
827    #[test]
828    fn test_ln_mvlgamma() {
829        assert!(super::mvlgamma(2, f64::NAN).is_nan());
830        assert_almost_eq!(super::mvlgamma(0, 1.0), 0.0, 1e-12);
831        assert_almost_eq!(super::mvlgamma(2, 1.0), 1.1447298858494002, 1e-12);
832        assert_almost_eq!(super::mvlgamma(3, 1.5), 2.1686775340635553, 1e-12);
833        assert_almost_eq!(super::mvlgamma(3, 150.0 + 1.0e-12), 1794.2387481112528, 1e-12);
834        assert_almost_eq!(super::mvlgamma(11, 14.5), 225.2071467353132, 1e-12);
835    }
836}
statrs/function/gamma.rs

statrs/function/
gamma.rs
