rs-stats 3.0.0

//! # Error Function (erf)
//!
//! This module implements the error function, a special function that occurs in probability,
//! statistics, and partial differential equations.
//!
//! ## Mathematical Definition
//! The error function is defined as:
//!
//! erf(x) = (2/√π) ∫₀ˣ e^(-t²) dt
//!
//! ## Key Properties
//! - erf(-x) = -erf(x) (odd function)
//! - erf(0) = 0
//! - erf(∞) = 1
//! - erf(-∞) = -1
//!
//! ## Implementation Details
//! Computed via the regularised lower incomplete gamma function:
//!   erf(x) = sign(x) · P(1/2, x²)
//! This delegates to [`crate::utils::special_functions::regularized_incomplete_gamma`]
//! and is accurate to ~1e-12 across the real line — substantially better
//! than the previous Abramowitz-Stegun 7.1.26 approximation (max error 1.5e-7).

use crate::error::{StatsError, StatsResult};
use crate::utils::special_functions::regularized_incomplete_gamma;
use num_traits::ToPrimitive;

/// Calculate the error function (erf) of a value
///
/// The error function is related to the normal distribution and is used
/// in probability calculations.
///
/// # Arguments
/// * `x` - The value at which to evaluate the error function
///
/// # Returns
/// The value of the error function at x
///
/// # Examples
/// ```
/// use rs_stats::prob::erf;
///
/// // Calculate erf(1.0)
/// let x: f64 = 1.0;
/// let result = erf(x).unwrap();
/// assert!((result - 0.842_700_792_949_714_9).abs() < 1e-12);
///
/// // Verify symmetry property
/// assert!((erf(x).unwrap() + erf(-x).unwrap()).abs() < 1e-12);
/// ```
#[inline]
pub fn erf<T>(x: T) -> StatsResult<f64>
where
    T: ToPrimitive,
{
    let x = x.to_f64().ok_or_else(|| StatsError::ConversionError {
        message: "prob::erf: Failed to convert x to f64".to_string(),
    })?;
    if x == 0.0 {
        return Ok(0.0);
    }
    if x.is_nan() {
        return Ok(f64::NAN);
    }
    if x.is_infinite() {
        return Ok(if x > 0.0 { 1.0 } else { -1.0 });
    }
    let sign = x.signum();
    // erf(x) = sign(x) · P(1/2, x²) via the canonical incomplete gamma.
    Ok(sign * regularized_incomplete_gamma(0.5, x * x))
}

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

    #[test]
    fn test_erf_special_cases() {
        assert!((erf(f64::INFINITY).unwrap() - 1.0).abs() < 1e-10);
        assert!((erf(f64::NEG_INFINITY).unwrap() + 1.0).abs() < 1e-10);
        assert!(erf(f64::NAN).unwrap().is_nan());
    }

    #[test]
    fn test_erf_against_known_values() {
        let test_cases = vec![
            (-3.0, -0.999977909503),
            (-2.0, -0.995322265019),
            (-1.0, -0.842700792950),
            (0.0, 0.0),
            (0.5, 0.520499877813),
            (1.0, 0.842700792950),
            (2.0, 0.995322265019),
            (3.0, 0.999977909503),
        ];

        for (x, expected) in test_cases {
            let actual = erf(x).unwrap();
            assert!(
                (actual - expected).abs() < 1e-6,
                "For x = {}, expected {}, but got {}",
                x,
                expected,
                actual
            );
        }
    }

    #[test]
    fn test_erf_symmetry() {
        let x = 0.7;
        let actual = erf(x).unwrap() + erf(-x).unwrap();
        assert!(
            actual.abs() < 1e-10,
            "erf(x) + erf(-x) should be 0.0, but got {}",
            actual
        );
    }

    #[test]
    fn test_erf_limits() {
        // Test erf approaching its limits
        assert!((erf(10.0).unwrap() - 1.0).abs() < 1e-15); // erf(x) -> 1 as x -> +inf
        assert!((erf(-10.0).unwrap() + 1.0).abs() < 1e-15); // erf(x) -> -1 as x -> -inf
    }

    #[test]
    fn test_erf_large_negative() {
        let x = -8.0;
        let actual = erf(x).unwrap();
        assert!(
            (actual + 1.0).abs() < 1e-10,
            "For large negative x, erf(x) should be close to -1.0, but got {}",
            actual
        );
    }
}