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

statrs/function/
gamma.rs
