sophy 0.2.0

//! # Special Mathematical Functions Module
//!
//! This module provides implementations of special mathematical functions that appear
//! frequently in advanced mathematics, physics, and engineering applications.
//!
//! ## Available Functions
//!
//! ### Special Functions
//! - [`gamma`]: Gamma function Γ(x)
//! - [`zeta`]: Riemann zeta function ζ(s)
//! - [`erf`]: Error function erf(x)
//! - [`eta`]: Dirichlet eta function η(s)
//! - [`sigma`]: Sum of divisors function σ(n)
//!
//! ## Usage Examples
//!
//! ```rust
//! use sophy::specials::{gamma, zeta, erf};
//!
//! // Gamma function: Γ(5) = 4! = 24
//! let factorial_4 = gamma(5.0);
//!
//! // Riemann zeta function: ζ(2) = π²/6
//! let zeta_2 = zeta(2.0);
//!
//! // Error function
//! let error_func = erf(1.0);
//! ```
//!
//! ## Mathematical Background
//!
//! These functions are fundamental in:
//! - **Probability theory**: Error function in normal distributions
//! - **Number theory**: Zeta and eta functions in prime number studies
//! - **Statistical mechanics**: Gamma function in probability distributions
//! - **Combinatorics**: Gamma function as generalization of factorial

use crate::base::numbers::PI;

/// Gamma function Γ(x) - generalization of factorial to real numbers
///
/// The gamma function extends the factorial function to real and complex numbers.
/// For positive integers n, Γ(n) = (n-1)!
///
/// ## Mathematical Definition
///
/// Γ(x) = ∫₀^∞ t^(x-1) * e^(-t) dt
///
/// ## Properties
///
/// - Γ(n) = (n-1)! for positive integers n
/// - Γ(x+1) = x * Γ(x) (recurrence relation)
/// - Γ(1/2) = √π
///
/// ## Arguments
///
/// * `x` - Input value (must be positive)
///
/// ## Returns
///
/// The value of Γ(x)
///
/// ## Examples
///
/// ```rust
/// use sophy::specials::gamma;
///
/// // Γ(5) = 4! = 24
/// assert!((gamma(5.0) - 24.0).abs() < 1e-10);
///
/// // Γ(1) = 0! = 1
/// assert!((gamma(1.0) - 1.0).abs() < 1e-10);
///
/// // Γ(1/2) = √π
/// assert!((gamma(0.5) - std::f64::consts::PI.sqrt()).abs() < 1e-10);
/// ```
///
/// ## Implementation Notes
///
/// This implementation uses Lanczos approximation for computational efficiency
/// while maintaining high precision for x > 0.5.
pub fn gamma(x: f64) -> f64 {
    // Handle special cases
    if x <= 0.0 {
        panic!("Gamma function undefined for non-positive values");
    }

    // Use recurrence relation to shift x into range [1, 2)
    if x < 1.0 {
        return gamma(x + 1.0) / x;
    }

    // For x >= 2, use recurrence relation Γ(x) = (x-1) * Γ(x-1)
    if x >= 2.0 {
        return (x - 1.0) * gamma(x - 1.0);
    }

    // Lanczos approximation for x in [1, 2)
    const G: f64 = 7.0;
    const COEFFICIENTS: [f64; 9] = [
        0.999_999_999_999_809_9,
        676.520_368_121_885_1,
        -1_259.139_216_722_402_8,
        771.323_428_777_653_1,
        -176.615_029_162_140_6,
        12.507_343_278_686_905,
        -0.138_571_095_265_720_12,
        9.984_369_578_019_572e-6,
        1.505_632_735_149_311_6e-7,
    ];

    let z = x - 1.0;
    let mut a = COEFFICIENTS[0];

    for (i, &coeff) in COEFFICIENTS.iter().enumerate().skip(1) {
        a += coeff / (z + i as f64);
    }

    let t = z + G + 0.5;
    (2.0 * PI).sqrt() * t.powf(z + 0.5) * (-t).exp() * a
}

/// Riemann zeta function ζ(s)
///
/// The Riemann zeta function is one of the most important functions in number theory,
/// defined as the infinite series ζ(s) = Σ(1/n^s) for n from 1 to infinity.
///
/// ## Mathematical Definition
///
/// ζ(s) = Σ_{n=1}^∞ 1/n^s for Re(s) > 1
///
/// ## Special Values
///
/// - ζ(2) = π²/6 ≈ 1.6449340668
/// - ζ(4) = π⁴/90 ≈ 1.0823232337
/// - ζ(0) = -1/2
///
/// ## Arguments
///
/// * `s` - Input value (must be > 1 for convergence)
///
/// ## Returns
///
/// The value of ζ(s)
///
/// ## Examples
///
/// ```rust
/// use sophy::specials::zeta;
///
/// // ζ(2) = π²/6
/// let zeta_2 = zeta(2.0);
/// let expected = std::f64::consts::PI.powi(2) / 6.0;
/// assert!((zeta_2 - expected).abs() < 1e-10);
///
/// // ζ(4) = π⁴/90
/// let zeta_4 = zeta(4.0);
/// let expected_4 = std::f64::consts::PI.powi(4) / 90.0;
/// assert!((zeta_4 - expected_4).abs() < 1e-8);
/// ```
pub fn zeta(s: f64) -> f64 {
    if s <= 1.0 {
        panic!("Zeta function implementation requires s > 1");
    }

    // For known values, return exact results
    if (s - 2.0).abs() < 1e-15 {
        return PI * PI / 6.0;
    }
    if (s - 4.0).abs() < 1e-15 {
        return PI.powi(4) / 90.0;
    }

    // Series approximation: ζ(s) = Σ(1/n^s)
    let mut sum = 0.0;
    let mut n = 1.0f64;
    let tolerance = 1e-15;

    loop {
        let term = 1.0 / n.powf(s);
        if term < tolerance {
            break;
        }
        sum += term;
        n += 1.0;

        // Prevent infinite loops
        if n > 1_000_000.0 {
            break;
        }
    }

    sum
}

/// Error function erf(x)
///
/// The error function is a special function that appears in probability theory,
/// statistics, and partial differential equations. It's closely related to the
/// normal distribution.
///
/// ## Mathematical Definition
///
/// erf(x) = (2/√π) * ∫₀^x e^(-t²) dt
///
/// ## Properties
///
/// - erf(0) = 0
/// - erf(∞) = 1
/// - erf(-x) = -erf(x) (odd function)
/// - erf'(x) = (2/√π) * e^(-x²)
///
/// ## Arguments
///
/// * `x` - Input value
///
/// ## Returns
///
/// The value of erf(x)
///
/// ## Examples
///
/// ```rust
/// use sophy::specials::erf;
///
/// // erf(0) = 0
/// assert!(erf(0.0).abs() < 1e-15);
///
/// // erf is an odd function
/// let x = 1.5;
/// assert!((erf(-x) + erf(x)).abs() < 1e-15);
///
/// // erf(∞) approaches 1
/// assert!((erf(5.0) - 1.0).abs() < 1e-10);
/// ```
pub fn erf(x: f64) -> f64 {
    // Handle special cases
    if x == 0.0 {
        return 0.0;
    }

    // Use odd function property: erf(-x) = -erf(x)
    if x < 0.0 {
        return -erf(-x);
    }

    // For large x, erf(x) approaches 1
    if x > 5.0 {
        return 1.0;
    }

    // Abramowitz and Stegun approximation
    // Maximum error: 1.5 × 10^(-7)
    const A1: f64 = 0.254829592;
    const A2: f64 = -0.284496736;
    const A3: f64 = 1.421413741;
    const A4: f64 = -1.453152027;
    const A5: f64 = 1.061405429;
    const P: f64 = 0.3275911;

    let t = 1.0 / (1.0 + P * x);
    let poly = A1 * t + A2 * t * t + A3 * t * t * t + A4 * t.powi(4) + A5 * t.powi(5);

    1.0 - poly * (-x * x).exp()
}

/// Dirichlet eta function η(s)
///
/// The Dirichlet eta function is related to the Riemann zeta function and is
/// defined as the alternating zeta function.
///
/// ## Mathematical Definition
///
/// η(s) = Σ_{n=1}^∞ (-1)^(n+1) / n^s = (1 - 2^(1-s)) * ζ(s)
///
/// ## Properties
///
/// - η(1) = ln(2) ≈ 0.693147
/// - η(s) = (1 - 2^(1-s)) * ζ(s)
/// - Converges for Re(s) > 0
///
/// ## Arguments
///
/// * `s` - Input value (must be > 0)
///
/// ## Returns
///
/// The value of η(s)
///
/// ## Examples
///
/// ```rust
/// use sophy::specials::eta;
///
/// // η(1) = ln(2)
/// let eta_1 = eta(1.0);
/// assert!((eta_1 - 2.0_f64.ln()).abs() < 1e-10);
///
/// // η(2) = π²/12
/// let eta_2 = eta(2.0);
/// let expected = std::f64::consts::PI.powi(2) / 12.0;
/// assert!((eta_2 - expected).abs() < 1e-10);
/// ```
pub fn eta(s: f64) -> f64 {
    if s <= 0.0 {
        panic!("Eta function implementation requires s > 0");
    }

    // Special case: η(1) = ln(2)
    if (s - 1.0).abs() < 1e-15 {
        return 2.0_f64.ln();
    }

    // For s != 1, use relation: η(s) = (1 - 2^(1-s)) * ζ(s)
    if s > 1.0 {
        let factor = 1.0 - 2.0_f64.powf(1.0 - s);
        return factor * zeta(s);
    }

    // Direct series calculation for 0 < s <= 1
    let mut sum = 0.0;
    let mut sign = 1.0;
    let tolerance = 1e-15;

    for n in 1..=1000000 {
        let term = sign / (n as f64).powf(s);
        if term.abs() < tolerance {
            break;
        }
        sum += term;
        sign *= -1.0;
    }

    sum
}

/// Sum of divisors function σ(n)
///
/// The sum of divisors function returns the sum of all positive divisors of n,
/// including 1 and n itself.
///
/// ## Mathematical Definition
///
/// σ(n) = Σ d where d divides n
///
/// ## Properties
///
/// - σ(1) = 1
/// - σ(p) = p + 1 for prime p
/// - σ is multiplicative: σ(mn) = σ(m)σ(n) if gcd(m,n) = 1
///
/// ## Arguments
///
/// * `n` - Input positive integer
///
/// ## Returns
///
/// The sum of all divisors of n
///
/// ## Examples
///
/// ```rust
/// use sophy::specials::sigma;
///
/// // σ(1) = 1
/// assert_eq!(sigma(1), 1);
///
/// // σ(6) = 1 + 2 + 3 + 6 = 12
/// assert_eq!(sigma(6), 12);
///
/// // σ(12) = 1 + 2 + 3 + 4 + 6 + 12 = 28
/// assert_eq!(sigma(12), 28);
/// ```
pub fn sigma(n: u64) -> u64 {
    if n == 0 {
        panic!("Sum of divisors undefined for n = 0");
    }

    if n == 1 {
        return 1;
    }

    let mut sum = 0;
    let sqrt_n = (n as f64).sqrt() as u64;

    for i in 1..=sqrt_n {
        if n.is_multiple_of(i) {
            sum += i;
            // Add the corresponding divisor n/i, but avoid double counting
            if i != n / i {
                sum += n / i;
            }
        }
    }

    sum
}