sparse_ir_basis/

special_functions.rs

//! Special functions implementation ported from C++ libsparseir
//!
//! This module provides high-precision implementations of special functions
//! used in the SparseIR library, particularly spherical Bessel functions
//! and related mathematical functions.

use std::f64;
use std::f64::consts::PI;

/// sqrt(π/2) - used frequently in spherical Bessel calculations
const SQPIO2: f64 = 1.253_314_137_315_500_3;

/// sqrt(2π) - prefactor of Stirling's series for the Gamma function
const SQ2PI: f64 = 2.506_628_274_631_000_7;

/// Above this argument the Gamma function uses Stirling's series
const GAMMA_STIRLING_MIN: f64 = 11.5;

/// Γ(x) exceeds `f64::MAX` for x > 171.624..., so every argument above this
/// bound overflows (Stirling's series already rounds to +∞ just below it)
const GAMMA_OVERFLOW_ARG: f64 = 172.0;

/// |Γ(x)| is below the smallest subnormal for every non-integer x < -184.
/// Beyond |x| = 200 the reflected Stirling series is not evaluated; its
/// intermediates stay finite up to about |x| = 256.
const GAMMA_UNDERFLOW_ARG: f64 = 200.0;

/// Maximum number of iterations for continued fractions
const MAX_ITER: usize = 5000;

/// Evaluate polynomial using Horner's method
fn evalpoly(x: f64, coeffs: &[f64]) -> f64 {
    let mut result = 0.0;
    for &coeff in coeffs.iter().rev() {
        result = result * x + coeff;
    }
    result
}

/// Compute sin(π*x) with exact argument reduction
///
/// `x = n + r` with `n = x.round()` and `|r| <= 1/2` is exact in floating
/// point, so `sin(πx) = (-1)^n sin(πr)` is accurate to a few ulp for every
/// finite `x` and is exactly zero if and only if `x` is an integer. The
/// unreduced `(PI * x).sin()` is neither: the rounding error of `PI * x` grows
/// with `|x|`, and `(PI * -1.0).sin()` is -1.2e-16 rather than 0.
fn sinpi(x: f64) -> f64 {
    let n = x.round();
    let s = (PI * (x - n)).sin();
    if n.rem_euclid(2.0) == 0.0 { s } else { -s }
}

/// Stirling's series for `x > GAMMA_STIRLING_MIN`
///
/// Returns `(v, t, w)` with `Γ(x) = SQ2PI * v * t * w`, where
/// `v = x^(x/2 - 1/4)` and `t = v / e^x` split `x^(x - 1/2) e^(-x)` so that it
/// does not overflow before Γ(x) itself, and `w` is the asymptotic series.
fn gamma_stirling_factors(x: f64) -> (f64, f64, f64) {
    let coefs = [
        1.0,
        8.333_333_333_333_331e-2,
        3.472_222_222_230_075e-3,
        -2.681_327_161_876_304_3e-3,
        -2.294_719_747_873_185_4e-4,
        7.840_334_842_744_753e-4,
        6.989_332_260_623_193e-5,
        -5.950_237_554_056_33e-4,
        -2.363_848_809_501_759e-5,
        7.147_391_378_143_611e-4,
    ];
    let w = evalpoly(1.0 / x, &coefs);
    let v = x.powf(0.5 * x - 0.25);
    (v, v / x.exp(), w)
}

/// Gamma function Γ(x) for real `x`
///
/// Adapted from the C++ libsparseir `gamma_func` (SpM-lab/libsparseir,
/// `backend/cxx/src/specfuncs.cpp` at commit 4bc58ea), which follows
/// `gamma(::Float64)` of Bessels.jl v0.2.8 (`src/gamma.jl`), itself adapted
/// from the Cephes Mathematical Library by Stephen L. Moshier. For `x > 0` it
/// uses Stirling's series above 11.5 and otherwise a rational approximation on
/// `[2, 3)` reached with `Γ(x + 1) = xΓ(x)`. It deviates from the C++ source in
/// three places: the Stirling prefactor is `sqrt(2π)` as in Bessels.jl (the
/// C++ code uses `sqrt(π/2)`, which halves every result above 11.5); `sin(πx)`
/// uses exact argument reduction, so that the poles are detected and accuracy
/// holds next to them; and no input throws or panics.
///
/// A negative non-integer `x < -1` uses the reflection formula
/// `Γ(x) = π / (sin(πx) Γ(1 - x))`, evaluated as `π / (sin(πx) |x| Γ(|x|))`
/// so that `1 - x` is never rounded. For `-1 < x < 0` the recurrence
/// `Γ(x) = Γ(1 + x) / x` is used instead, because `|x| sin(πx)` underflows
/// for tiny `|x|`.
///
/// Wherever Γ(x) is a normal `f64`, the relative error is a few ulp.
///
/// Special values follow C99 `tgamma`:
///
/// - `x = ±0` (pole): `±∞`, the sign of zero selecting the side of the pole.
/// - `x` a negative integer (a pole without a signed limit) or `x = -∞`: NaN.
/// - `x = +∞`, and `x > 171.62...` where Γ(x) exceeds `f64::MAX`: `+∞`.
/// - Non-integer `x < -184`, where |Γ(x)| is below the smallest subnormal:
///   `±0` carrying the sign of Γ(x), which is the sign of `sin(πx)`.
/// - `x = NaN`: NaN.
pub fn gamma_func(x: f64) -> f64 {
    if x.is_nan() {
        return x;
    }
    if x <= 0.0 && x == x.floor() {
        // Poles at zero and at the negative integers (x = -∞ also lands here).
        // Only a zero tells the side of the pole, through its sign.
        return if x == 0.0 { 1.0 / x } else { f64::NAN };
    }
    if x > 0.0 {
        return gamma_positive(x);
    }
    if x > -1.0 {
        // Γ(x) = Γ(1 + x) / x; the reflection below would underflow for tiny |x|.
        return gamma_positive(1.0 + x) / x;
    }

    // Reflection for non-integer x < -1: s = |x| sin(πx) is nonzero, and
    // Γ(x) = π / (s Γ(|x|)) has the sign of s.
    let ax = -x;
    let s = ax * sinpi(x);
    if ax <= GAMMA_STIRLING_MIN {
        return PI / (s * gamma_positive(ax));
    }
    if ax > GAMMA_UNDERFLOW_ARG {
        return 0.0_f64.copysign(s);
    }
    // Γ(|x|) = SQ2PI * v * t * w overflows f64 for |x| > 171.62 while Γ(x)
    // next to a pole is still a normal number down to x ≈ -175. Dividing in
    // stages keeps every intermediate finite.
    let (v, t, w) = gamma_stirling_factors(ax);
    PI / (s * SQ2PI * w * v) / t
}

/// Γ(x) for `x > 0`, including `x = +∞`
fn gamma_positive(x: f64) -> f64 {
    if x > GAMMA_STIRLING_MIN {
        if x > GAMMA_OVERFLOW_ARG {
            // Also avoids ∞/∞ = NaN in `t` once e^x overflows (x > 709.78).
            return f64::INFINITY;
        }
        let (v, t, w) = gamma_stirling_factors(x);
        return SQ2PI * v * t * w;
    }

    let p = [
        1.0,
        8.378_004_301_573_126e-1,
        3.629_515_436_640_239_3e-1,
        1.113_062_816_019_361_6e-1,
        2.385_363_243_461_108_3e-2,
        4.092_666_828_394_036e-3,
        4.542_931_960_608_009_3e-4,
        4.212_760_487_471_622e-5,
    ];

    let q = [
        1.0,
        4.150_160_950_588_455_7e-1,
        -2.243_510_905_670_329_2e-1,
        -4.633_887_671_244_534e-2,
        2.773_706_565_840_073e-2,
        -7.955_933_682_494_738e-4,
        -1.237_799_246_653_152_3e-3,
        2.346_584_059_160_635e-4,
        -1.397_148_517_476_170_5e-5,
    ];

    // Shift x into [2, 3) with Γ(x + 1) = xΓ(x).
    let mut x = x;
    let mut z = 1.0;
    while x >= 3.0 {
        x -= 1.0;
        z *= x;
    }

    while x < 2.0 {
        z /= x;
        x += 1.0;
    }

    if x == 2.0 {
        return z;
    }

    x -= 2.0;
    let p_val = evalpoly(x, &p);
    let q_val = evalpoly(x, &q);

    z * p_val / q_val
}

/// Cylindrical Bessel function of the first kind, J_nu(x)
///
/// Uses the series expansion:
///   J_nu(x) = sum_{m=0}^∞ (-1)^m / (m! * Gamma(nu+m+1)) * (x/2)^(2m+nu)
pub fn cyl_bessel_j(nu: f64, x: f64) -> f64 {
    let eps = f64::EPSILON;
    let mut _sum = 0.0;
    let mut term = (x / 2.0).powf(nu) / gamma_func(nu + 1.0);
    _sum = term;

    for m in 1..1000 {
        term *= -(x * x / 4.0) / (m as f64 * (nu + m as f64));
        _sum += term;
        if term.abs() < _sum.abs() * eps {
            break;
        }
    }

    _sum
}

/// Spherical Bessel function j_n(x) using the relation:
///   j_n(x) = sqrt(pi/(2x)) * J_{n+1/2}(x)
fn spherical_bessel_j_generic(nu: f64, x: f64) -> f64 {
    SQPIO2 * cyl_bessel_j(nu + 0.5, x) / x.sqrt()
}

/// Approximation for small x
fn spherical_bessel_j_small_args(nu: f64, x: f64) -> f64 {
    if x == 0.0 {
        return if nu == 0.0 { 1.0 } else { 0.0 };
    }

    let x2 = (x * x) / 4.0;
    let coef = [
        1.0,
        -1.0 / (1.5 + nu), // 3/2 + nu
        -1.0 / (5.0 + nu),
        -1.0 / ((21.0 / 2.0) + nu), // 21/2 + nu
        -1.0 / (18.0 + nu),
    ];

    let a = SQPIO2 / (gamma_func(1.5 + nu) * 2.0_f64.powf(nu + 0.5));
    x.powf(nu) * a * evalpoly(x2, &coef)
}

/// Determines when the small-argument expansion is accurate
fn spherical_bessel_j_small_args_cutoff(nu: f64, x: f64) -> bool {
    (x * x) / (4.0 * nu + 110.0) < f64::EPSILON
}

/// Computes the continued-fraction for the ratio J_{nu}(x) / J_{nu-1}(x)
fn bessel_j_ratio_jnu_jnum1(n: f64, x: f64) -> f64 {
    let xinv = 1.0 / x;
    let xinv2 = 2.0 * xinv;
    let mut d = x / (2.0 * n);
    let mut a = d;
    let mut h = a;
    let mut b = (2.0 * n + 2.0) * xinv;

    for _i in 0..MAX_ITER {
        d = 1.0 / (b - d);
        a *= b * d - 1.0;
        h += a;
        b += xinv2;

        if (a / h).abs() <= f64::EPSILON {
            break;
        }
    }

    h
}

/// Computes forward recurrence for spherical Bessel y.
/// Returns a pair: (sY_{n-1}, sY_n)
fn spherical_bessel_y_forward_recurrence(nu: i32, x: f64) -> (f64, f64) {
    let xinv = 1.0 / x;
    let s = x.sin();
    let c = x.cos();
    let mut s_y0 = -c * xinv;
    let mut s_y1 = xinv * (s_y0 - s);
    let mut nu_start = 1.0;

    while nu_start < nu as f64 + 0.5 {
        let temp = s_y1;
        s_y1 = (2.0 * nu_start + 1.0) * xinv * s_y1 - s_y0;
        s_y0 = temp;
        nu_start += 1.0;
    }

    (s_y0, s_y1)
}

/// Uses forward recurrence if stable; otherwise uses spherical Bessel y recurrence
fn spherical_bessel_j_recurrence(nu: i32, x: f64) -> f64 {
    if x >= nu as f64 {
        let xinv = 1.0 / x;
        let s = x.sin();
        let c = x.cos();
        let mut s_j0 = s * xinv;
        let mut s_j1 = (s_j0 - c) * xinv;
        let mut nu_start = 1.0;

        while nu_start < nu as f64 + 0.5 {
            let temp = s_j1;
            s_j1 = (2.0 * nu_start + 1.0) * xinv * s_j1 - s_j0;
            s_j0 = temp;
            nu_start += 1.0;
        }

        s_j0
    } else {
        // For x < nu, use the alternative method
        // This should return j_nu(x), not j_nu-1(x)
        let (s_ynm1, s_yn) = spherical_bessel_y_forward_recurrence(nu, x);
        let h = bessel_j_ratio_jnu_jnum1(nu as f64 + 1.5, x);
        1.0 / (x * x * (h * s_ynm1 - s_yn))
    }
}

/// Selects the proper method for computing j_n(x) for positive arguments
fn spherical_bessel_j_positive_args(nu: i32, x: f64) -> f64 {
    if spherical_bessel_j_small_args_cutoff(nu as f64, x) {
        spherical_bessel_j_small_args(nu as f64, x)
    } else if (x >= nu as f64 && nu < 250) || (x < nu as f64 && nu < 60) {
        // Use recurrence for both x >= nu and x < nu (when nu < 60)
        spherical_bessel_j_recurrence(nu, x)
    } else {
        spherical_bessel_j_generic(nu as f64, x)
    }
}

/// Main function to calculate spherical Bessel function of the first kind
///
/// This is the main entry point that matches the C++ sphericalbesselj function
pub fn spherical_bessel_j(n: i32, x: f64) -> f64 {
    // Handle negative arguments
    if x < 0.0 {
        panic!("sphericalBesselJ requires non-negative x");
    }

    // Handle negative orders: j_{-n}(x) = (-1)^n * j_n(x)
    if n < 0 {
        let result = spherical_bessel_j_positive_args(-n, x);
        if n % 2 == 0 { result } else { -result }
    } else {
        spherical_bessel_j_positive_args(n, x)
    }
}

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

    #[test]
    fn test_gamma_function() {
        // Test some known values
        assert!((gamma_func(1.0) - 1.0).abs() < 1e-10);
        assert!((gamma_func(2.0) - 1.0).abs() < 1e-10);
        assert!((gamma_func(3.0) - 2.0).abs() < 1e-10);
        assert!((gamma_func(4.0) - 6.0).abs() < 1e-10);

        // Test half-integer values
        assert!((gamma_func(0.5) - 1.7724538509055159).abs() < 1e-10); // sqrt(π)
    }

    // Reference values below are correctly rounded from MPFR, evaluated at the
    // exact f64 arguments: Julia 1.12.5, SpecialFunctions.jl 2.8.3,
    // `setprecision(BigFloat, 256)`, Γ values as
    // `repr(Float64(gamma(BigFloat(x))))`.

    /// Relative tolerance of `gamma_func` against correctly rounded references.
    ///
    /// Error model: the Cephes rational approximation and Stirling series
    /// contribute a few ulp (Bessels.jl tests the same algorithm at 7 eps
    /// against BigFloat), the reduced `sin(πx)` about one ulp, and the
    /// reflection a handful of roundings, i.e. at most about 16 ulp (3.6e-15).
    /// 1e-14 leaves a ~3x margin for platform `pow`/`exp`/`sin` differences
    /// while staying far below the O(1) relative errors of the defects guarded
    /// against here (wrong sign, Γ(|x|) instead of Γ(x), a factor of 2).
    const GAMMA_RTOL: f64 = 1e-14;

    /// Relative tolerance for the Bessel functions built on `gamma_func`. The
    /// Γ error (at most about 16 ulp, see above) is a common factor of every
    /// series term; `powf` and the per-term roundings add a few ulp, amplified
    /// by at most the cancellation ratio sum|terms| / |sum| (about 3, for
    /// J_{-1/2}(1)). Together at most about 30 ulp (6.7e-15), with the same
    /// ~3x margin.
    const BESSEL_RTOL: f64 = 2e-14;

    fn rel_err(got: f64, want: f64) -> f64 {
        ((got - want) / want).abs()
    }

    #[test]
    fn test_gamma_negative_non_integer_reference_values() {
        // The half-integer rows equal the closed forms -2√π, 4√π/3, -8√π/15
        // and 16√π/105.
        let cases = [
            (-0.5, -3.544907701811032),
            (-1.5, 2.363271801207355),
            (-2.5, -0.9453087204829419),
            (-3.5, 0.2700882058522691),
            (-0.1, -10.686287021193193),
            (-1.0e-200, -1.0e200),
            (-0.999, -1000.4241966812758),
            (-1.001, 999.5786270024664),
            (-3.0 + 2f64.powi(-30), -1.789569708760196e8),
            (-10.1, -2.2134165830856185e-6),
            (-11.3, 4.656958619058061e-8),
            (-11.7, 1.7234064143490278e-8),
            (-12.5, -1.836606483859281e-9),
            (-20.5, -2.834656574391335e-19),
            (-30.7, -1.3275373492818893e-33),
            (-100.25, -1.503087709322751e-158),
            (-150.9, -1.9468126352925122e-264),
            (-170.5, -3.3127395215386074e-308),
            // Γ(175) overflows f64, but this Γ(x) is a normal number.
            (-175.0 + 2f64.powi(-40), -9.778221578627872e-307),
        ];
        for (x, want) in cases {
            let got = gamma_func(x);
            let err = rel_err(got, want);
            assert!(
                err <= GAMMA_RTOL,
                "gamma_func({x:e}) = {got:e}, expected {want:e} (relative error {err:e})"
            );
        }
    }

    #[test]
    fn test_gamma_positive_reference_values() {
        // Rows above 11.5 use Stirling's series; 13, 20 and 25 are factorials.
        let cases = [
            (0.5, 1.772453850905516),
            (2.5, 1.329340388179137),
            (11.5, 1.1899423083962249e7),
            (11.6, 1.5131318919703094e7),
            (12.5, 1.3684336546556586e8),
            (13.0, 4.790016e8),
            (20.0, 1.21645100408832e17),
            (25.0, 6.204484017332394e23),
            (50.5, 4.29046291235196e63),
            (100.0, 9.332621544394415e155),
            (170.5, 5.56209241456e305),
            (171.5, 9.4833675668248e307),
        ];
        for (x, want) in cases {
            let got = gamma_func(x);
            let err = rel_err(got, want);
            assert!(
                err <= GAMMA_RTOL,
                "gamma_func({x:e}) = {got:e}, expected {want:e} (relative error {err:e})"
            );
        }
    }

    #[test]
    fn test_gamma_recurrence() {
        // Γ(x + 1) = x Γ(x) for both signs of x, straddling the switch to
        // Stirling's series at |x| = 11.5 and approaching the poles. Dyadic
        // fractions keep x and x + 1 exact; |x| < 170 keeps every value a
        // normal f64. Tolerance: two evaluations within 16 ulp each plus one
        // product (about 33 ulp, 7.3e-15), with the same ~3x margin as
        // GAMMA_RTOL.
        const RECURRENCE_RTOL: f64 = 2e-14;
        let fracs = [
            2f64.powi(-20),
            0.125,
            0.25,
            0.5,
            0.75,
            0.875,
            1.0 - 2f64.powi(-20),
        ];
        for k in 0..170 {
            for f in fracs {
                for x in [-(k as f64 + f), k as f64 + f] {
                    let lhs = gamma_func(x + 1.0);
                    let rhs = x * gamma_func(x);
                    let err = rel_err(rhs, lhs);
                    assert!(
                        err <= RECURRENCE_RTOL,
                        "Γ({x:e} + 1) = {lhs:e} but x Γ(x) = {rhs:e} (relative error {err:e})"
                    );
                }
            }
        }
    }

    #[test]
    fn test_gamma_poles_and_special_values() {
        // Negative integers are poles without a signed limit: NaN, never a panic.
        for x in [
            -1.0,
            -2.0,
            -3.0,
            -171.0,
            -1e10,
            -(2f64.powi(52)),
            -1e300,
            f64::MIN,
        ] {
            let got = gamma_func(x);
            assert!(got.is_nan(), "gamma_func({x:e}) = {got:e}, expected NaN");
        }
        // The pole at zero is approached from the side given by the sign of zero.
        assert_eq!(gamma_func(0.0), f64::INFINITY);
        assert_eq!(gamma_func(-0.0), f64::NEG_INFINITY);
        assert_eq!(gamma_func(f64::INFINITY), f64::INFINITY);
        assert!(gamma_func(f64::NEG_INFINITY).is_nan());
        assert!(gamma_func(f64::NAN).is_nan());

        // Overflow: Γ(x) exceeds f64::MAX for x > 171.62...
        for x in [171.7, 172.0, 709.0, 710.0, 1000.0, 1e300] {
            let got = gamma_func(x);
            assert_eq!(got, f64::INFINITY, "gamma_func({x:e}) = {got:e}");
        }

        // Γ(-171.5) is subnormal and must not be flushed to zero. Its final
        // rounding onto the subnormal grid costs half a unit of the spacing;
        // allow two on top of the relative bound.
        let (got, want) = (gamma_func(-171.5), 1.9316265431712e-310);
        assert!(
            (got - want).abs() <= GAMMA_RTOL * want + 2.0 * f64::from_bits(1),
            "gamma_func(-171.5) = {got:e}, expected {want:e}"
        );

        // Below x = -184, |Γ(x)| is under the smallest subnormal for every
        // non-integer x: the result is a zero with the sign of Γ(x).
        for (x, negative) in [
            (-184.5, true),
            (-200.5, true),
            (-201.5, false),
            (-1000.25, true),
            (-(1e15 + 0.5), true),
        ] {
            let got = gamma_func(x);
            assert!(
                got == 0.0 && got.is_sign_negative() == negative,
                "gamma_func({x:e}) = {got:e}, expected {}0",
                if negative { "-" } else { "+" }
            );
        }
    }

    #[test]
    fn test_cyl_bessel_j_orders_reaching_gamma_reflection_and_stirling() {
        // cyl_bessel_j(nu, x) divides by gamma_func(nu + 1): nu < -1 reaches
        // negative non-integer arguments, nu > 10.5 the Stirling branch.
        // References: the power series
        // sum((-1)^m (x/2)^(2m+nu) / (gamma(m+1) gamma(m+nu+1)) for m in 0:80)
        // in 256-bit BigFloat (same Julia setup as above).
        let cases = [
            (-0.5, 1.0, 0.4310988680183761),
            (-1.5, 1.0, -1.1024955751601793),
            (-2.5, 2.0, 0.8282206324443038),
            (-12.7, 1.0, 3.9444125287326807e11),
            (11.5, 1.0, 2.4730845703448897e-12),
            (12.0, 2.0, 1.9326951487239857e-9),
        ];
        for (nu, x, want) in cases {
            let got = cyl_bessel_j(nu, x);
            let err = rel_err(got, want);
            assert!(
                err <= BESSEL_RTOL,
                "J_{nu}({x}) = {got:e}, expected {want:e} (relative error {err:e})"
            );
        }
    }

    #[test]
    fn test_spherical_bessel_j_small_args_high_order() {
        // Below the small-argument cutoff, j_n(x) uses gamma_func(n + 1.5),
        // which lies in the Stirling branch for n >= 11. References:
        // j_n(x) = sqrt(pi/(2x)) J_{n+1/2}(x) from the 256-bit series above.
        let x = 1e-8;
        let cases = [
            (11, 3.162213889372793e-100),
            (12, 1.2648855557491175e-109),
            (13, 4.6847613175893235e-119),
            (14, 1.6154349370997668e-128),
            (15, 5.2110804422573123e-138),
        ];
        for (n, want) in cases {
            assert!(spherical_bessel_j_small_args_cutoff(n as f64, x));
            let got = spherical_bessel_j(n, x);
            let err = rel_err(got, want);
            assert!(
                err <= BESSEL_RTOL,
                "j_{n}({x:e}) = {got:e}, expected {want:e} (relative error {err:e})"
            );
        }
    }

    #[test]
    fn test_cylindrical_bessel_j() {
        // Test known values
        let j0_1 = cyl_bessel_j(0.0, 1.0);
        let expected_j0_1 = 0.765_197_686_557_966_6;
        assert!((j0_1 - expected_j0_1).abs() < 1e-10);

        let j1_1 = cyl_bessel_j(1.0, 1.0);
        let expected_j1_1 = 0.440_050_585_744_933_5;
        assert!((j1_1 - expected_j1_1).abs() < 1e-10);
    }

    #[test]
    fn test_spherical_bessel_j_basic() {
        // Test j_0(x) = sin(x)/x for x != 0
        let x = 1.0;
        let j0 = spherical_bessel_j(0, x);
        let expected_j0 = x.sin() / x;
        println!("j_0({}) = {}, expected = {}", x, j0, expected_j0);
        assert!((j0 - expected_j0).abs() < 1e-10);

        // Test j_1(x) = sin(x)/x² - cos(x)/x
        let j1 = spherical_bessel_j(1, x);
        let expected_j1 = x.sin() / (x * x) - x.cos() / x;
        println!("j_1({}) = {}, expected = {}", x, j1, expected_j1);
        assert!((j1 - expected_j1).abs() < 1e-10);

        // Test j_0(0) = 1
        let j0_zero = spherical_bessel_j(0, 0.0);
        println!("j_0(0) = {}, expected = 1.0", j0_zero);
        assert!((j0_zero - 1.0).abs() < 1e-10);

        // Test j_n(0) = 0 for n > 0
        let j1_zero = spherical_bessel_j(1, 0.0);
        println!("j_1(0) = {}, expected = 0.0", j1_zero);
        assert!(j1_zero.abs() < 1e-10);
    }

    #[test]
    fn test_spherical_bessel_j_various_values() {
        // Test various values to ensure accuracy
        // These are reference values from mathematical tables
        let test_cases = [
            (0, 0.1, 0.9983341664682815),
            (0, 0.5, 0.958_851_077_208_406),
            (0, 1.0, 0.8414709848078965),
            (0, 2.0, 0.4546487134128409),
            (0, 5.0, -0.1917848549326277),
            // Corrected expected values for j_1 using analytical formulas
            (1, 0.1, 0.0333000128900053),
            (1, 0.5, 0.1625370306360665),
            (1, 1.0, 0.3011686789397568),
            // j_1(2) = sin(2)/4 - cos(2)/2 = 0.43539777497999166
            (1, 2.0, 0.43539777497999166),
            // j_1(5) = sin(5)/25 - cos(5)/5 = -0.0950894080791708
            (1, 5.0, -0.0950894080791708),
        ];

        for (n, x, expected) in test_cases {
            let result = spherical_bessel_j(n, x);
            println!(
                "j_{}({}) = {}, expected = {}, diff = {}",
                n,
                x,
                result,
                expected,
                (result - expected).abs()
            );

            // For now, just check that the result is finite and reasonable
            assert!(
                result.is_finite(),
                "j_{}({}) should be finite, got {}",
                n,
                x,
                result
            );

            // Check accuracy with more lenient tolerance for now
            if (result - expected).abs() > 1e-6 {
                println!(
                    "WARNING: j_{}({}) accuracy issue: got {}, expected {}, diff = {}",
                    n,
                    x,
                    result,
                    expected,
                    (result - expected).abs()
                );
            }
        }
    }

    #[test]
    fn test_debug_spherical_bessel_j() {
        // Debug specific problematic cases
        println!("=== Debug j_1(0.1) ===");
        let x = 0.1;
        let n = 1;

        // Check which method is being used
        let cutoff = spherical_bessel_j_small_args_cutoff(n as f64, x);
        println!("small_args_cutoff: {}", cutoff);

        if cutoff {
            let small_result = spherical_bessel_j_small_args(n as f64, x);
            println!("small_args result: {}", small_result);
        }

        let recurrence_result = spherical_bessel_j_recurrence(n, x);
        println!("recurrence result: {}", recurrence_result);

        let generic_result = spherical_bessel_j_generic(n as f64, x);
        println!("generic result: {}", generic_result);

        let final_result = spherical_bessel_j(n, x);
        println!("final result: {}", final_result);

        // Expected: j_1(0.1) = sin(0.1)/0.1^2 - cos(0.1)/0.1 = 0.033300...
        let expected = x.sin() / (x * x) - x.cos() / x;
        println!("expected (sin(x)/x^2 - cos(x)/x): {}", expected);
    }

    #[test]
    fn test_spherical_bessel_j_large_values() {
        // Test large values to ensure stability
        let large_x = 100.0;
        let j0_large = spherical_bessel_j(0, large_x);
        println!("j_0({}) = {}", large_x, j0_large);
        assert!(j0_large.is_finite());

        let j1_large = spherical_bessel_j(1, large_x);
        println!("j_1({}) = {}", large_x, j1_large);
        assert!(j1_large.is_finite());
    }

    #[test]
    fn test_spherical_bessel_j_cpp_style_high_order() {
        // Test high-order spherical Bessel functions like C++ implementation
        // Reference values from Julia (same as C++ test)
        // julia> using Bessels
        // julia> for i in 0:15; println(sphericalbesselj(i, 1.)); end
        let refs = [
            0.8414709848078965,
            0.30116867893975674,
            0.06203505201137386,
            0.009006581117112517,
            0.0010110158084137527,
            9.256115861125818e-5,
            7.156936310087086e-6,
            4.790134198739489e-7,
            2.82649880221473e-8,
            1.4913765025551456e-9,
            7.116552640047314e-11,
            3.09955185479008e-12,
            1.2416625969871055e-13,
            4.604637677683788e-15,
            1.5895759875169764e-16,
            5.1326861154437626e-18,
        ];

        let x = 1.0;
        for (l, &expected) in refs.iter().enumerate() {
            let expected: f64 = expected;
            let result = spherical_bessel_j(l as i32, x);

            // Use same tolerance as C++ Approx: relative error 1e-6, absolute error 1e-12
            let relative_tolerance = 1e-6;
            let absolute_tolerance = 1e-12;

            // Check relative error for non-zero expected values
            if expected.abs() > absolute_tolerance {
                let relative_error = (result - expected).abs() / expected.abs();
                assert!(
                    relative_error <= relative_tolerance,
                    "j_{}({}) relative error too large: {} > {}",
                    l,
                    x,
                    relative_error,
                    relative_tolerance
                );
            } else {
                // For very small expected values, check absolute error
                assert!(
                    (result - expected).abs() <= absolute_tolerance,
                    "j_{}({}) absolute error too large: {} > {}",
                    l,
                    x,
                    (result - expected).abs(),
                    absolute_tolerance
                );
            }

            // Check that result is finite
            assert!(
                result.is_finite(),
                "j_{}({}) should be finite, got {}",
                l,
                x,
                result
            );
        }
    }

    #[test]
    fn test_spherical_bessel_j_zero_argument() {
        // Test behavior at x = 0
        // j_0(0) = 1, j_n(0) = 0 for n > 0
        let j0_zero = spherical_bessel_j(0, 0.0);
        assert!(
            (j0_zero - 1.0).abs() < 1e-15,
            "j_0(0) should be 1, got {}",
            j0_zero
        );

        for n in 1..=10 {
            let jn_zero = spherical_bessel_j(n, 0.0);
            assert!(
                jn_zero.abs() < 1e-15,
                "j_{}(0) should be 0, got {}",
                n,
                jn_zero
            );
        }
    }

    #[test]
    fn test_spherical_bessel_j_negative_orders() {
        // Test behavior for negative orders: j_{-n}(x) = (-1)^n * j_n(x)
        for n in -5..0 {
            let result_neg = spherical_bessel_j(n, 1.0);
            let result_pos = spherical_bessel_j(-n, 1.0);
            let expected = if n % 2 == 0 { result_pos } else { -result_pos };
            assert!(
                (result_neg - expected).abs() < 1e-15,
                "j_{}(1.0) should be {} for negative n, got {}",
                n,
                expected,
                result_neg
            );
        }
    }

    #[test]
    fn test_spherical_bessel_j_small_arguments() {
        // Test very small arguments to ensure numerical stability
        let small_x_values = [1e-10, 1e-8, 1e-6, 1e-4, 1e-2];

        for &x in &small_x_values {
            for n in 0..=5 {
                let result = spherical_bessel_j(n, x);
                assert!(result.is_finite(), "j_{}({}) should be finite", n, x);

                // For very small x, j_n(x) ≈ x^n / (2n+1)!!
                if x < 1e-6 && n < 3 {
                    let expected_approx = x.powi(n) / double_factorial(2 * n + 1);
                    let relative_error = (result - expected_approx).abs() / expected_approx.abs();
                    assert!(
                        relative_error < 1e-6,
                        "j_{}({}) small argument approximation failed: relative error = {}",
                        n,
                        x,
                        relative_error
                    );
                }
            }
        }
    }

    #[test]
    fn test_spherical_bessel_j_large_arguments() {
        // Test large arguments to ensure asymptotic behavior
        let large_x_values = [10.0, 50.0, 100.0, 500.0];

        for &x in &large_x_values {
            for n in 0..=3 {
                let result = spherical_bessel_j(n, x);
                assert!(result.is_finite(), "j_{}({}) should be finite", n, x);

                // For large x, j_n(x) ≈ sin(x - n*π/2) / x
                let expected_approx = (x - (n as f64) * PI / 2.0).sin() / x;
                let relative_error =
                    (result - expected_approx).abs() / expected_approx.abs().max(1e-10);

                if x > 50.0 && relative_error > 1e-2 {
                    println!(
                        "WARNING: j_{}({}) large argument approximation: got {}, expected ≈ {}, relative error = {}",
                        n, x, result, expected_approx, relative_error
                    );
                }
            }
        }
    }

    /// Helper function for double factorial: n!!
    fn double_factorial(n: i32) -> f64 {
        if n <= 0 {
            1.0
        } else if n % 2 == 0 {
            // Even: n!! = 2^(n/2) * (n/2)!
            let half_n = n / 2;
            2.0_f64.powi(half_n) * factorial(half_n)
        } else {
            // Odd: n!! = n! / (2^((n-1)/2) * ((n-1)/2)!)
            let half_n_minus_1 = (n - 1) / 2;
            factorial(n) / (2.0_f64.powi(half_n_minus_1) * factorial(half_n_minus_1))
        }
    }

    /// Helper function for factorial: n!
    fn factorial(n: i32) -> f64 {
        if n <= 1 {
            1.0
        } else {
            let mut result = 1.0;
            for i in 2..=n {
                result *= i as f64;
            }
            result
        }
    }
}