mathhook_core/functions/special/
error_functions.rs

1//! Error function implementations
2
3use crate::core::{Expression, MathConstant, Number};
4
5/// Error function erf(x)
6///
7/// # Mathematical Definition
8///
9/// erf(x) = (2/√π) ∫₀ˣ e^(-t²) dt
10///
11/// # Arguments
12///
13/// * `arg` - Expression to compute error function of
14///
15/// # Returns
16///
17/// Error function expression
18///
19/// # Examples
20///
21/// ```
22/// use mathhook_core::functions::special::error_functions::erf;
23/// use mathhook_core::expr;
24///
25/// let result = erf(&expr!(0));
26/// assert_eq!(result, expr!(0));
27/// ```
28pub fn erf(arg: &Expression) -> Expression {
29    match arg {
30        Expression::Number(Number::Integer(0)) => Expression::integer(0),
31        Expression::Constant(MathConstant::Infinity) => Expression::integer(1),
32        Expression::Constant(MathConstant::NegativeInfinity) => Expression::integer(-1),
33        Expression::Number(Number::Float(f)) => Expression::float(erf_approx(*f)),
34        Expression::Number(Number::Integer(i)) => Expression::float(erf_approx(*i as f64)),
35        _ => Expression::function("erf", vec![arg.clone()]),
36    }
37}
38
39/// Complementary error function erfc(x)
40///
41/// # Mathematical Definition
42///
43/// erfc(x) = 1 - erf(x) = (2/√π) ∫ₓ^∞ e^(-t²) dt
44///
45/// # Arguments
46///
47/// * `arg` - Expression to compute complementary error function of
48///
49/// # Returns
50///
51/// Complementary error function expression
52///
53/// # Examples
54///
55/// ```
56/// use mathhook_core::functions::special::error_functions::erfc;
57/// use mathhook_core::expr;
58///
59/// let result = erfc(&expr!(0));
60/// assert_eq!(result, expr!(1));
61/// ```
62pub fn erfc(arg: &Expression) -> Expression {
63    match arg {
64        Expression::Number(Number::Integer(0)) => Expression::integer(1),
65        Expression::Constant(MathConstant::Infinity) => Expression::integer(0),
66        Expression::Constant(MathConstant::NegativeInfinity) => Expression::integer(2),
67        Expression::Number(Number::Float(f)) => Expression::float(erfc_approx(*f)),
68        Expression::Number(Number::Integer(i)) => Expression::float(erfc_approx(*i as f64)),
69        _ => Expression::function("erfc", vec![arg.clone()]),
70    }
71}
72
73fn erf_approx(x: f64) -> f64 {
74    if x.is_infinite() {
75        return if x.is_sign_positive() { 1.0 } else { -1.0 };
76    }
77
78    let a1 = 0.254829592;
79    let a2 = -0.284496736;
80    let a3 = 1.421413741;
81    let a4 = -1.453152027;
82    let a5 = 1.061405429;
83    let p = 0.3275911;
84
85    let sign = if x < 0.0 { -1.0 } else { 1.0 };
86    let x = x.abs();
87
88    let t = 1.0 / (1.0 + p * x);
89    let y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * (-x * x).exp();
90
91    sign * y
92}
93
94fn erfc_approx(x: f64) -> f64 {
95    1.0 - erf_approx(x)
96}
97
98#[cfg(test)]
99mod tests {
100    use super::*;
101
102    #[test]
103    fn test_erf_zero() {
104        assert_eq!(erf(&Expression::integer(0)), Expression::integer(0));
105    }
106
107    #[test]
108    fn test_erfc_zero() {
109        assert_eq!(erfc(&Expression::integer(0)), Expression::integer(1));
110    }
111
112    #[test]
113    fn test_erf_infinity() {
114        assert_eq!(
115            erf(&Expression::constant(MathConstant::Infinity)),
116            Expression::integer(1)
117        );
118    }
119
120    #[test]
121    fn test_erf_negative_infinity() {
122        assert_eq!(
123            erf(&Expression::constant(MathConstant::NegativeInfinity)),
124            Expression::integer(-1)
125        );
126    }
127
128    #[test]
129    fn test_erfc_infinity() {
130        assert_eq!(
131            erfc(&Expression::constant(MathConstant::Infinity)),
132            Expression::integer(0)
133        );
134    }
135}