injective_math/fp_decimal/
factorial.rs

1use crate::fp_decimal::FPDecimal;
2use std::str::FromStr;
3impl FPDecimal {
4    fn _factorial(x: FPDecimal) -> FPDecimal {
5        if x.is_zero() {
6            FPDecimal::ONE
7        } else if x < FPDecimal::ZERO {
8            x * FPDecimal::_factorial(x + FPDecimal::ONE)
9        } else {
10            x * FPDecimal::_factorial(x - FPDecimal::ONE)
11        }
12    }
13    pub fn factorial(&self) -> FPDecimal {
14        FPDecimal::_factorial(*self)
15    }
16
17    pub fn gamma(x: FPDecimal) -> FPDecimal {
18        let table: [FPDecimal; 30] = [
19            FPDecimal::from_str("1.00000000000000000000").unwrap(),
20            FPDecimal::from_str("0.57721566490153286061").unwrap(),
21            FPDecimal::from_str("-0.65587807152025388108").unwrap(),
22            FPDecimal::from_str("-0.04200263503409523553").unwrap(),
23            FPDecimal::from_str("0.16653861138229148950").unwrap(),
24            FPDecimal::from_str("-0.04219773455554433675").unwrap(),
25            FPDecimal::from_str("-0.00962197152787697356").unwrap(),
26            FPDecimal::from_str("0.00721894324666309954").unwrap(),
27            FPDecimal::from_str("-0.00116516759185906511").unwrap(),
28            FPDecimal::from_str("-0.00021524167411495097").unwrap(),
29            FPDecimal::from_str("0.00012805028238811619").unwrap(),
30            FPDecimal::from_str("-0.00002013485478078824").unwrap(),
31            FPDecimal::from_str("-0.00000125049348214267").unwrap(),
32            FPDecimal::from_str("0.00000113302723198170").unwrap(),
33            FPDecimal::from_str("-0.00000020563384169776").unwrap(),
34            FPDecimal::from_str("0.00000000611609510448").unwrap(),
35            FPDecimal::from_str("0.00000000500200764447").unwrap(),
36            FPDecimal::from_str("-0.00000000118127457049").unwrap(),
37            FPDecimal::from_str("0.00000000010434267117").unwrap(),
38            FPDecimal::from_str("0.00000000000778226344").unwrap(),
39            FPDecimal::from_str("-0.00000000000369680562").unwrap(),
40            FPDecimal::from_str("0.00000000000051003703").unwrap(),
41            FPDecimal::from_str("-0.00000000000002058326").unwrap(),
42            FPDecimal::from_str("-0.00000000000000534812").unwrap(),
43            FPDecimal::from_str("0.00000000000000122678").unwrap(),
44            FPDecimal::from_str("-0.00000000000000011813").unwrap(),
45            FPDecimal::from_str("0.00000000000000000119").unwrap(),
46            FPDecimal::from_str("0.00000000000000000141").unwrap(),
47            FPDecimal::from_str("-0.00000000000000000023").unwrap(),
48            FPDecimal::from_str("0.00000000000000000002").unwrap(),
49        ];
50        let y = x - FPDecimal::ONE;
51
52        let mut sm = table[table.len() - 1];
53        for i in (0..table.len() - 2).rev() {
54            sm = sm * y + table[i];
55        }
56        FPDecimal::ONE / sm
57    }
58}
59#[cfg(test)]
60mod tests {
61
62    use crate::FPDecimal;
63    use std::str::FromStr;
64
65    #[test]
66    fn test_factorial_nine() {
67        let nine = FPDecimal::NINE;
68        assert_eq!(nine, FPDecimal::from_str("9").unwrap());
69        assert_eq!(FPDecimal::from_str("362880").unwrap(), nine.factorial());
70    }
71
72    #[test]
73    fn test_factorial_negative_nine() {
74        let negative_nine = -FPDecimal::NINE;
75        assert_eq!(negative_nine, FPDecimal::from_str("-9").unwrap());
76        assert_eq!(FPDecimal::from_str("-362880").unwrap(), negative_nine.factorial());
77    }
78}