statrs 0.19.0

Statistical computing library for Rust
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112

//! Provides functions related to exponential calculations

use crate::consts;
#[cfg(not(feature = "std"))]
use num_traits::Float as _;

/// Computes the generalized Exponential Integral function
/// where `x` is the argument and `n` is the integer power of the
/// denominator term.
///
/// Returns `None` if `x < 0.0` or the computation could not
/// converge after 100 iterations
///
/// # Remarks
///
/// This implementation follows the derivation in
///
/// _"Handbook of Mathematical Functions, Applied Mathematics Series, Volume
/// 55"_ - Abramowitz, M., and Stegun, I.A 1964
///
/// AND
///
/// _"Advanced mathematical methods for scientists and engineers"_ - Bender,
/// Carl M.; Steven A. Orszag (1978). page 253
///
/// The continued fraction approach is used for `x > 1.0` while the taylor
/// series expansions is used for `0.0 < x <= 1`.
// TODO: Add examples
pub fn integral(x: f64, n: u64) -> Option<f64> {
    let eps = 0.00000000000000001;
    let max_iter = 100;
    let nf64 = n as f64;
    let near_f64min = 1e-100; // needs very small value that is not quite as small as f64 min

    // special cases
    if n == 0 {
        return Some((-x).exp() / x);
    }
    if x == 0.0 {
        return Some(1.0 / (nf64 - 1.0));
    }

    if x > 1.0 {
        let mut b = x + nf64;
        let mut c = 1.0 / near_f64min;
        let mut d = 1.0 / b;
        let mut h = d;
        for i in 1..max_iter + 1 {
            let a = -(i as f64) * (nf64 - 1.0 + i as f64);
            b += 2.0;
            d = 1.0 / (a * d + b);
            c = b + a / c;
            let del = c * d;
            h *= del;
            if (del - 1.0).abs() < eps {
                return Some(h * (-x).exp());
            }
        }
        None
    } else {
        let mut factorial = 1.0;
        let mut result = if n - 1 != 0 {
            1.0 / (nf64 - 1.0)
        } else {
            -x.ln() - consts::EULER_MASCHERONI
        };
        for i in 1..max_iter + 1 {
            factorial *= -x / i as f64;
            let del = if i != n - 1 {
                -factorial / (i as f64 - nf64 + 1.0)
            } else {
                let mut psi = -consts::EULER_MASCHERONI;
                for ii in 1..n {
                    psi += 1.0 / ii as f64;
                }
                factorial * (-x.ln() + psi)
            };
            result += del;
            if del.abs() < result.abs() * eps {
                return Some(result);
            }
        }
        None
    }
}

#[rustfmt::skip]
#[cfg(test)]
mod tests {
    use crate::prec;
    use super::*;

    #[test]
    fn test_integral() {
        assert_eq!(integral(0.001, 1).unwrap(), 6.33153936413614904);
        prec::assert_abs_diff_eq!(integral(0.1, 1).unwrap(), 1.82292395841939059, epsilon = 1e-15);
        assert_eq!(integral(1.0, 1).unwrap(), 0.219383934395520286);
        prec::assert_abs_diff_eq!(integral(2.0, 1).unwrap(), 0.0489005107080611248, epsilon = 1e-15);
        prec::assert_abs_diff_eq!(integral(2.5, 1).unwrap(), 0.0249149178702697399, epsilon = 1e-15);
        prec::assert_abs_diff_eq!(integral(10.0, 1).unwrap(), 4.15696892968532464e-06, epsilon = 1e-20);
        assert_eq!(integral(0.001, 2).unwrap(), 0.992668960469238915);
        prec::assert_abs_diff_eq!(integral(0.1, 2).unwrap(), 0.722545022194020392, epsilon = 1e-15);
        prec::assert_abs_diff_eq!(integral(1.0, 2).unwrap(), 0.148495506775922048, epsilon = 1e-16);
        prec::assert_abs_diff_eq!(integral(2.0, 2).unwrap(), 0.0375342618204904527, epsilon = 1e-16);
        prec::assert_abs_diff_eq!(integral(10.0, 2).unwrap(), 3.830240465631608e-06, epsilon = 1e-20);
        assert_eq!(integral(0.001, 0).unwrap(), 999.000499833375);
        assert_eq!(integral(0.1, 0).unwrap(), 9.048374180359595);
        prec::assert_abs_diff_eq!(integral(1.0, 0).unwrap(), 0.3678794411714423, epsilon = 1e-16);
        assert_eq!(integral(2.0, 0).unwrap(), 0.06766764161830635);
        assert_eq!(integral(10.0, 0).unwrap(), 4.539992976248485e-06);
    }
}