injective_math/fp_decimal/
factorial.rs1use 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}