abax 0.1.20

A lightweight Rust library providing high-precision mathematical constants and special functions, including Bernoulli numbers, Riemann Zeta values, and robust incomplete gamma functions.
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
113
114

use crate::erfc::erfc;

/// Calculates the scaled complementary error function `erfcx(x) = exp(x^2) * erfc(x)`.
///
/// The scaled complementary error function is useful for large values of `x` where
/// `erfc(x)` would underflow to zero and `exp(x^2)` would overflow to infinity.
/// As `x` increases, `erfcx(x)` decays slowly toward zero as `1 / (x * sqrt(pi))`.
///
/// # Numerical stability
/// - For `x < 0`, the function grows very rapidly and is calculated as
///   `erfcx(x) = 2 * exp(x^2) - erfcx(|x|)`.
/// - For small `x >= 0`, it is calculated using the product `exp(x^2) * erfc(x)`.
/// - For large `x`, it is calculated using a continued fraction expansion to
///   maintain precision and avoid overflow/underflow.
///
/// # Examples
/// ```
/// use abax::erfcx;
///
/// assert!((erfcx(0.0) - 1.0).abs() < 1e-15);
/// assert!((erfcx(1.0) - 0.4275835761558071).abs() < 1e-15);
/// ```
pub fn erfcx(x: f64) -> f64 {
    if x.is_nan() {
        return f64::NAN;
    }
    if x.is_infinite() {
        return if x.is_sign_positive() {
            0.0
        } else {
            f64::INFINITY
        };
    }

    if x < 0.0 {
        let x2 = x * x;
        // exp(x^2) overflows for x^2 > log(f64::MAX) ~= 709.78
        if x2 > 709.78 {
            return f64::INFINITY;
        }
        // Reflection formula: erfcx(x) = exp(x^2) * (2 - erfc(|x|)) = 2*exp(x^2) - erfcx(|x|)
        return 2.0 * x2.exp() - erfcx(-x);
    }

    // For positive x, erfcx is monotonically decreasing from 1 to 0.
    if x < 2.5 {
        // Standard definition is safe for small x.
        return (x * x).exp() * erfc(x);
    }

    // For large x, use Laplace's continued fraction expansion:
    // erfcx(x) = (1/sqrt(pi)) * [1 / (x + (1/2 / (x + (1 / (x + (3/2 / (x + ...)))))))]
    erfcx_cf(x)
}

/// Evaluates the continued fraction for erfcx(x) using the modified Lentz's method.
fn erfcx_cf(x: f64) -> f64 {
    const INV_SQRT_PI: f64 = 0.5641895835477563;
    const TINY: f64 = 1.0e-100;

    let mut f = TINY;
    let mut c = f;
    let mut d = 0.0;

    for i in 1..200 {
        // Coefficients for the continued fraction:
        // a_1 = 1, b_1 = x
        // a_n = (n-1)/2, b_n = x  for n > 1
        let a = if i == 1 { 1.0 } else { (i as f64 - 1.0) * 0.5 };
        let b = x;

        d = b + a * d;
        if d == 0.0 { d = TINY; }
        c = b + a / c;
        if c == 0.0 { c = TINY; }

        d = 1.0 / d;
        let delta = c * d;
        f *= delta;

        if (delta - 1.0).abs() < 1e-16 {
            break;
        }
    }

    f * INV_SQRT_PI
}

#[cfg(test)]
mod tests {
    use super::*;

    fn assert_close(actual: f64, expected: f64, tol: f64) {
        let diff = (actual - expected).abs();
        assert!(diff <= tol, "actual={actual:?}, expected={expected:?}, diff={diff:?}");
    }

    #[test]
    fn test_special_cases() {
        assert!(erfcx(f64::NAN).is_nan());
        assert_eq!(erfcx(f64::INFINITY), 0.0);
        assert_eq!(erfcx(f64::NEG_INFINITY), f64::INFINITY);
        assert_close(erfcx(0.0), 1.0, 1e-15);
    }

    #[test]
    fn test_known_values() {
        assert_close(erfcx(0.5), 0.615697223530000, 1e-14);
        assert_close(erfcx(1.0), 0.4275835761558071, 1e-15);
        assert_close(erfcx(5.0), 0.1116123285081272, 1e-15);
        assert_close(erfcx(10.0), 0.0561459483563720, 1e-15);
        assert_close(erfcx(-1.0), 5.008977599118573, 1e-14);
    }
}