statrs/function/
exponential.rs

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

statrs/function/
exponential.rs
