statrs/function/
exponential.rs

1//! Provides functions related to exponential calculations
2
3use crate::consts;
4
5/// Computes the generalized Exponential Integral function
6/// where `x` is the argument and `n` is the integer power of the
7/// denominator term.
8///
9/// Returns `None` if `x < 0.0` or the computation could not
10/// converge after 100 iterations
11///
12/// # Remarks
13///
14/// This implementation follows the derivation in
15///
16/// _"Handbook of Mathematical Functions, Applied Mathematics Series, Volume
17/// 55"_ - Abramowitz, M., and Stegun, I.A 1964
18///
19/// AND
20///
21/// _"Advanced mathematical methods for scientists and engineers"_ - Bender,
22/// Carl M.; Steven A. Orszag (1978). page 253
23///
24/// The continued fraction approach is used for `x > 1.0` while the taylor
25/// series expansions is used for `0.0 < x <= 1`.
26// TODO: Add examples
27pub fn integral(x: f64, n: u64) -> Option<f64> {
28    let eps = 0.00000000000000001;
29    let max_iter = 100;
30    let nf64 = n as f64;
31    let near_f64min = 1e-100; // needs very small value that is not quite as small as f64 min
32
33    // special cases
34    if n == 0 {
35        return Some((-1.0 * x).exp() / x);
36    }
37    if x == 0.0 {
38        return Some(1.0 / (nf64 - 1.0));
39    }
40
41    if x > 1.0 {
42        let mut b = x + nf64;
43        let mut c = 1.0 / near_f64min;
44        let mut d = 1.0 / b;
45        let mut h = d;
46        for i in 1..max_iter + 1 {
47            let a = -1.0 * i as f64 * (nf64 - 1.0 + i as f64);
48            b += 2.0;
49            d = 1.0 / (a * d + b);
50            c = b + a / c;
51            let del = c * d;
52            h *= del;
53            if (del - 1.0).abs() < eps {
54                return Some(h * (-x).exp());
55            }
56        }
57        None
58    } else {
59        let mut factorial = 1.0;
60        let mut result = if n - 1 != 0 {
61            1.0 / (nf64 - 1.0)
62        } else {
63            -1.0 * x.ln() - consts::EULER_MASCHERONI
64        };
65        for i in 1..max_iter + 1 {
66            factorial *= -1.0 * x / i as f64;
67            let del = if i != n - 1 {
68                -factorial / (i as f64 - nf64 + 1.0)
69            } else {
70                let mut psi = -1.0 * consts::EULER_MASCHERONI;
71                for ii in 1..n {
72                    psi += 1.0 / ii as f64;
73                }
74                factorial * (-1.0 * x.ln() + psi)
75            };
76            result += del;
77            if del.abs() < result.abs() * eps {
78                return Some(result);
79            }
80        }
81        None
82    }
83}
84
85#[rustfmt::skip]
86#[cfg(test)]
87mod tests {
88    #[test]
89    fn test_integral() {
90        assert_eq!(super::integral(0.001, 1).unwrap(), 6.33153936413614904);
91        assert_almost_eq!(super::integral(0.1, 1).unwrap(), 1.82292395841939059, 1e-15);
92        assert_eq!(super::integral(1.0, 1).unwrap(), 0.219383934395520286);
93        assert_almost_eq!(super::integral(2.0, 1).unwrap(), 0.0489005107080611248, 1e-15);
94        assert_almost_eq!(super::integral(2.5, 1).unwrap(), 0.0249149178702697399, 1e-15);
95        assert_almost_eq!(super::integral(10.0, 1).unwrap(), 4.15696892968532464e-06, 1e-20);
96        assert_eq!(super::integral(0.001, 2).unwrap(), 0.992668960469238915);
97        assert_almost_eq!(super::integral(0.1, 2).unwrap(), 0.722545022194020392, 1e-15);
98        assert_almost_eq!(super::integral(1.0, 2).unwrap(), 0.148495506775922048, 1e-16);
99        assert_almost_eq!(super::integral(2.0, 2).unwrap(), 0.0375342618204904527, 1e-16);
100        assert_almost_eq!(super::integral(10.0, 2).unwrap(), 3.830240465631608e-06, 1e-20);
101        assert_eq!(super::integral(0.001, 0).unwrap(), 999.000499833375);
102        assert_eq!(super::integral(0.1, 0).unwrap(), 9.048374180359595);
103        assert_almost_eq!(super::integral(1.0, 0).unwrap(), 0.3678794411714423, 1e-16);
104        assert_eq!(super::integral(2.0, 0).unwrap(), 0.06766764161830635);
105        assert_eq!(super::integral(10.0, 0).unwrap(), 4.539992976248485e-06);
106    }
107}
statrs/function/exponential.rs

statrs/function/
exponential.rs
