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
113
114
115
116
117

//! Provides functions related to exponential calculations

use crate::{consts, Result, StatsError};

/// Computes the generalized Exponential Integral function
/// where `x` is the argument and `n` is the integer power of the
/// denominator term.
///
/// # Errors
///
/// Returns an error if `x < 0.0` or the computation could not
/// converge after 100 iterations
///
/// # Remarks
///
/// This implementation follows the derivation in
/// <br />
/// <div>
/// <i>"Handbook of Mathematical Functions, Applied Mathematics Series, Volume
/// 55"</i> - Abramowitz, M., and Stegun, I.A 1964
/// </div>
/// AND
/// <br />
/// <div>
/// <i>"Advanced mathematical methods for scientists and engineers" - Bender,
/// Carl M.; Steven A. Orszag (1978). page 253
/// </div>
/// <br />
/// The continued fraction approac is used for `x > 1.0` while the taylor
/// series expansions
/// is used for `0.0 < x <= 1`
///
/// # Examples
///
/// ```
/// ```
pub fn integral(x: f64, n: u64) -> Result<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 Ok((-1.0 * x).exp() / x);
    }
    if x == 0.0 {
        return Ok(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 = -1.0 * 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 Ok(h * (-x).exp());
            }
        }
        Err(StatsError::ComputationFailedToConverge)
    } else {
        let mut factorial = 1.0;
        let mut result = if n - 1 != 0 {
            1.0 / (nf64 - 1.0)
        } else {
            -1.0 * x.ln() - consts::EULER_MASCHERONI
        };
        for i in 1..max_iter + 1 {
            factorial *= -1.0 * x / i as f64;
            let del = if i != n - 1 {
                -factorial / (i as f64 - nf64 + 1.0)
            } else {
                let mut psi = -1.0 * consts::EULER_MASCHERONI;
                for ii in 1..n {
                    psi += 1.0 / ii as f64;
                }
                factorial * (-1.0 * x.ln() + psi)
            };
            result += del;
            if del.abs() < result.abs() * eps {
                return Ok(result);
            }
        }
        Err(StatsError::ComputationFailedToConverge)
    }
}

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