statrs/function/
factorial.rs

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

statrs/function/
factorial.rs
