expr_solver/
utils.rs

1const TAYLOR_COEFFICIENTS: [f64; 29] = [
2    -0.00000000000000000023,
3    0.00000000000000000141,
4    0.00000000000000000119,
5    -0.00000000000000011813,
6    0.00000000000000122678,
7    -0.00000000000000534812,
8    -0.00000000000002058326,
9    0.00000000000051003703,
10    -0.00000000000369680562,
11    0.00000000000778226344,
12    0.00000000010434267117,
13    -0.00000000118127457049,
14    0.00000000500200764447,
15    0.00000000611609510448,
16    -0.00000020563384169776,
17    0.00000113302723198170,
18    -0.00000125049348214267,
19    -0.00002013485478078824,
20    0.00012805028238811619,
21    -0.00021524167411495097,
22    -0.00116516759185906511,
23    0.00721894324666309954,
24    -0.00962197152787697356,
25    -0.04219773455554433675,
26    0.16653861138229148950,
27    -0.04200263503409523553,
28    -0.65587807152025388108,
29    0.57721566490153286061,
30    1.00000000000000000000,
31];
32
33const INITIAL_SUM: f64 = 0.00000000000000000002;
34
35// Calculates gamma of a f64.
36// https://en.wikipedia.org/wiki/Gamma_function
37fn gamma(x: f64) -> f64 {
38    TAYLOR_COEFFICIENTS
39        .iter()
40        .fold(INITIAL_SUM, |sum, coefficient| {
41            sum * (x - 1.0) + coefficient
42        })
43        .recip()
44}
45
46// Calculates simple factorial of a f64.
47fn simple_factorial(x: f64) -> f64 {
48    if x <= 1.0 {
49        return 1.0;
50    }
51
52    simple_factorial(x - 1.0) * x
53}
54
55/// Wrapper function to calculate factorial of a function.
56/// Returns simple factorial of x if the number can be converted to integer without loss.
57/// Returns gamma of x if the number cannot be converted to integer without loss.
58/// # Arguments
59/// * x - the number to find factorial of.
60pub fn factorial(x: f64) -> f64 {
61    log::trace!("getting factorial of {}", x);
62
63    // for integers.
64    if x.fract() == 0.0 {
65        return simple_factorial(x);
66    }
67
68    // for non integers.
69    gamma(x + 1_f64)
70}