statrs/function/
factorial.rs

1//! Provides functions related to factorial calculations (e.g. binomial
2//! coefficient, factorial, multinomial)
3
4use crate::error::StatsError;
5use crate::function::gamma;
6use crate::Result;
7use core::f64::INFINITY as INF;
8
9/// The maximum factorial representable
10/// by a 64-bit floating point without
11/// overflowing
12pub const MAX_FACTORIAL: usize = 170;
13
14/// Computes the factorial function `x -> x!` for
15/// `170 >= x >= 0`. All factorials larger than `170!`
16/// will overflow an `f64`.
17///
18/// # Remarks
19///
20/// Returns `f64::INFINITY` if `x > 170`
21pub fn factorial(x: u64) -> f64 {
22    let x = x as usize;
23    FCACHE.get(x).map_or(INF, |&fac| fac)
24}
25
26/// Computes the logarithmic factorial function `x -> ln(x!)`
27/// for `x >= 0`.
28///
29/// # Remarks
30///
31/// Returns `0.0` if `x <= 1`
32pub fn ln_factorial(x: u64) -> f64 {
33    let x = x as usize;
34    FCACHE
35        .get(x)
36        .map_or_else(|| gamma::ln_gamma(x as f64 + 1.0), |&fac| fac.ln())
37}
38
39/// Computes the binomial coefficient `n choose k`
40/// where `k` and `n` are non-negative values.
41///
42/// # Remarks
43///
44/// Returns `0.0` if `k > n`
45pub fn binomial(n: u64, k: u64) -> f64 {
46    if k > n {
47        0.0
48    } else {
49        (0.5 + (ln_factorial(n) - ln_factorial(k) - ln_factorial(n - k)).exp()).floor()
50    }
51}
52
53/// Computes the natural logarithm of the binomial coefficient
54/// `ln(n choose k)` where `k` and `n` are non-negative values
55///
56/// # Remarks
57///
58/// Returns `f64::NEG_INFINITY` if `k > n`
59pub fn ln_binomial(n: u64, k: u64) -> f64 {
60    if k > n {
61        f64::NEG_INFINITY
62    } else {
63        ln_factorial(n) - ln_factorial(k) - ln_factorial(n - k)
64    }
65}
66
67/// Computes the multinomial coefficient: `n choose n1, n2, n3, ...`
68///
69/// # Panics
70///
71/// If the elements in `ni` do not sum to `n`
72pub fn multinomial(n: u64, ni: &[u64]) -> f64 {
73    checked_multinomial(n, ni).unwrap()
74}
75
76/// Computes the multinomial coefficient: `n choose n1, n2, n3, ...`
77///
78/// # Errors
79///
80/// If the elements in `ni` do not sum to `n`
81pub fn checked_multinomial(n: u64, ni: &[u64]) -> Result<f64> {
82    let (sum, ret) = ni.iter().fold((0, ln_factorial(n)), |acc, &x| {
83        (acc.0 + x, acc.1 - ln_factorial(x))
84    });
85    if sum != n {
86        Err(StatsError::ContainerExpectedSumVar("ni", "n"))
87    } else {
88        Ok((0.5 + ret.exp()).floor())
89    }
90}
91
92// Initialization for pre-computed cache of 171 factorial
93// values 0!...170!
94lazy_static! {
95    static ref FCACHE: [f64; MAX_FACTORIAL + 1] = {
96        let mut fcache = [1.0; MAX_FACTORIAL + 1];
97        fcache
98            .iter_mut()
99            .enumerate()
100            .skip(1)
101            .fold(1.0, |acc, (i, elt)| {
102                let fac = acc * i as f64;
103                *elt = fac;
104                fac
105            });
106        fcache
107    };
108}
109
110#[cfg(test)]
111mod tests {
112    use super::*;
113
114    #[test]
115    fn test_factorial_and_ln_factorial() {
116        let mut fac = 1.0;
117        assert_eq!(factorial(0), fac);
118        for i in 1..171 {
119            fac *= i as f64;
120            assert_eq!(factorial(i), fac);
121            assert_eq!(ln_factorial(i), fac.ln());
122        }
123    }
124
125    #[test]
126    fn test_factorial_overflow() {
127        assert_eq!(factorial(172), INF);
128        assert_eq!(factorial(u64::MAX), INF);
129    }
130
131    #[test]
132    fn test_ln_factorial_does_not_overflow() {
133        assert_eq!(ln_factorial(1 << 10), 6078.2118847500501140);
134        assert_almost_eq!(ln_factorial(1 << 12), 29978.648060844048236, 1e-11);
135        assert_eq!(ln_factorial(1 << 15), 307933.81973375485425);
136        assert_eq!(ln_factorial(1 << 17), 1413421.9939462073242);
137    }
138
139    #[test]
140    fn test_binomial() {
141        assert_eq!(binomial(1, 1), 1.0);
142        assert_eq!(binomial(5, 2), 10.0);
143        assert_eq!(binomial(7, 3), 35.0);
144        assert_eq!(binomial(1, 0), 1.0);
145        assert_eq!(binomial(0, 1), 0.0);
146        assert_eq!(binomial(5, 7), 0.0);
147    }
148
149    #[test]
150    fn test_ln_binomial() {
151        assert_eq!(ln_binomial(1, 1), 1f64.ln());
152        assert_almost_eq!(ln_binomial(5, 2), 10f64.ln(), 1e-14);
153        assert_almost_eq!(ln_binomial(7, 3), 35f64.ln(), 1e-14);
154        assert_eq!(ln_binomial(1, 0), 1f64.ln());
155        assert_eq!(ln_binomial(0, 1), 0f64.ln());
156        assert_eq!(ln_binomial(5, 7), 0f64.ln());
157    }
158
159    #[test]
160    fn test_multinomial() {
161        assert_eq!(1.0, multinomial(1, &[1, 0]));
162        assert_eq!(10.0, multinomial(5, &[3, 2]));
163        assert_eq!(10.0, multinomial(5, &[2, 3]));
164        assert_eq!(35.0, multinomial(7, &[3, 4]));
165    }
166
167    #[test]
168    #[should_panic]
169    fn test_multinomial_bad_ni() {
170        multinomial(1, &[1, 1]);
171    }
172
173    #[test]
174    fn test_checked_multinomial_bad_ni() {
175        assert!(checked_multinomial(1, &[1, 1]).is_err());
176    }
177}
statrs/function/factorial.rs

statrs/function/
factorial.rs
