abax 0.1.49

#[inline]
fn evaluate_polynomial(coeffs: &[f64], x: f64) -> f64 {
    let mut result = *coeffs.last().unwrap();
    for &c in coeffs[..coeffs.len() - 1].iter().rev() {
        result = result * x + c;
    }
    result
}

#[inline]
pub(crate) fn erfinv_imp(p: f64, q: f64) -> f64 {
    // Piecewise rational approximations for 64-bit and smaller floating-point types.
    if p <= 0.5 {
        // Evaluate inverse erf using: x = p(p + 10)(Y + R(p)).
        const Y: f64 = 0.089_131_474_494_934_08;
        const P: [f64; 8] = [
            -0.000_508_781_949_658_280_7,
            -0.008_368_748_197_417_368,
            0.033_480_662_540_974_46,
            -0.012_692_614_766_297_403,
            -0.036_563_797_141_176_27,
            0.021_987_868_111_116_89,
            0.008_226_878_746_769_157,
            -0.005_387_729_650_712_429,
        ];
        const Q: [f64; 10] = [
            1.0,
            -0.970_005_043_303_290_6,
            -1.565_745_582_341_758_5,
            1.562_215_583_984_230_2,
            0.662_328_840_472_003,
            -0.712_289_023_415_428_5,
            -0.052_739_638_234_009_97,
            0.079_528_368_734_157_17,
            -0.002_333_937_593_741_900_3,
            0.000_886_216_390_456_424_7,
        ];
        let g = p * (p + 10.0);
        let r = evaluate_polynomial(&P, p) / evaluate_polynomial(&Q, p);
        g * (Y + r)
    } else if q >= 0.25 {
        // Rational approximation for 0.5 > q >= 0.25:
        // x = sqrt(-2*log(q)) / (Y + R(q - 0.25)).
        const Y: f64 = 2.249_481_201_171_875;
        const P: [f64; 9] = [
            -0.202_433_508_355_938_76,
            0.105_264_680_699_391_71,
            8.370_503_283_431_2,
            17.644_729_840_837_4,
            -18.851_064_805_871_427,
            -44.638_232_444_178_7,
            17.445_385_985_570_866,
            21.129_465_544_834_053,
            -3.671_922_547_077_293_6,
        ];
        const Q: [f64; 9] = [
            1.0,
            6.242_641_248_542_475,
            3.971_343_795_334_387,
            -28.660_818_049_980_003,
            -20.143_263_468_048_52,
            48.560_921_310_873_994,
            10.826_866_735_546_016,
            -22.643_693_341_313_97,
            1.721_147_657_612_002_8,
        ];
        let g = (-2.0 * q.ln()).sqrt();
        let xs = q - 0.25;
        let r = evaluate_polynomial(&P, xs) / evaluate_polynomial(&Q, xs);
        g / (Y + r)
    } else {
        // For q < 0.25, let x = sqrt(-log(q)) and approximate the
        // result as x * (Y + R(x - B)) over progressively smaller tails.
        let x = (-q.ln()).sqrt();
        if x < 3.0 {
            const Y: f64 = 0.807_220_458_984_375;
            const P: [f64; 11] = [
                -0.131_102_781_679_951_9,
                -0.163_794_047_193_317_06,
                0.117_030_156_341_995_25,
                0.387_079_738_972_604_34,
                0.337_785_538_912_035_9,
                0.142_869_534_408_157_17,
                0.029_015_791_000_532_906,
                0.002_145_589_953_888_053,
                -0.679_465_575_181_126_4e-6,
                0.285_225_331_782_217_06e-7,
                -0.681_149_956_853_777e-9,
            ];
            const Q: [f64; 8] = [
                1.0,
                3.466_254_072_425_672_4,
                5.381_683_457_070_069,
                4.778_465_929_458_438,
                2.593_019_216_236_202_7,
                0.848_854_343_457_902,
                0.152_264_338_295_331_78,
                0.011_059_242_293_464_891,
            ];
            let xs = x - 1.125;
            x * (Y + evaluate_polynomial(&P, xs) / evaluate_polynomial(&Q, xs))
        } else if x < 6.0 {
            const Y: f64 = 0.939_955_711_364_746_1;
            const P: [f64; 9] = [
                -0.035_035_378_718_317_8,
                -0.002_224_265_292_134_479_4,
                0.018_557_330_651_423_107,
                0.009_508_047_013_259_196,
                0.001_871_234_928_195_592_3,
                0.000_157_544_617_424_960_56,
                0.460_469_890_584_318e-5,
                -0.230_404_776_911_882_6e-9,
                0.266_339_227_425_782_03e-11,
            ];
            const Q: [f64; 7] = [
                1.0,
                1.365_334_981_755_406_3,
                0.762_059_164_553_623_4,
                0.220_091_105_764_131_25,
                0.034_158_914_367_094_77,
                0.002_638_616_766_570_16,
                0.764_675_292_302_794_5e-4,
            ];
            let xs = x - 3.0;
            x * (Y + evaluate_polynomial(&P, xs) / evaluate_polynomial(&Q, xs))
        } else if x < 18.0 {
            const Y: f64 = 0.983_628_273_010_253_9;
            const P: [f64; 9] = [
                -0.016_743_100_507_663_374,
                -0.001_129_514_387_455_802_8,
                0.001_056_288_621_524_929,
                0.000_209_386_317_487_588_08,
                0.149_624_783_758_342_37e-4,
                0.449_696_789_927_706_5e-6,
                0.462_596_163_522_878_6e-8,
                -0.281_128_735_628_831_8e-13,
                0.990_557_099_733_103_3e-16,
            ];
            const Q: [f64; 7] = [
                1.0,
                0.591_429_344_886_417_5,
                0.138_151_865_749_083_32,
                0.016_074_608_709_367_65,
                0.000_964_011_807_005_165_5,
                0.275_335_474_764_726e-4,
                0.282_243_172_016_108e-6,
            ];
            let xs = x - 6.0;
            x * (Y + evaluate_polynomial(&P, xs) / evaluate_polynomial(&Q, xs))
        } else if x < 44.0 {
            const Y: f64 = 0.997_145_652_770_996_1;
            const P: [f64; 8] = [
                -0.002_497_821_279_189_813,
                -0.779_190_719_229_054e-5,
                0.254_723_037_413_027_45e-4,
                0.162_397_777_342_510_92e-5,
                0.396_341_011_304_801_17e-7,
                0.411_632_831_190_944_2e-9,
                0.145_596_286_718_675_04e-11,
                -0.116_765_012_397_184_28e-17,
            ];
            const Q: [f64; 7] = [
                1.0,
                0.207_123_112_214_422_52,
                0.016_941_083_812_097_59,
                0.000_690_538_265_622_684_6,
                0.145_007_359_818_232_64e-4,
                0.144_437_756_628_144_16e-6,
                0.509_761_276_599_778_5e-9,
            ];
            let xs = x - 18.0;
            x * (Y + evaluate_polynomial(&P, xs) / evaluate_polynomial(&Q, xs))
        } else {
            const Y: f64 = 0.999_413_490_295_410_2;
            const P: [f64; 8] = [
                -0.000_539_042_911_019_078_6,
                -0.283_987_590_047_277_2e-6,
                0.899_465_114_892_291_5e-6,
                0.229_345_859_265_920_87e-7,
                0.225_561_444_863_500_15e-9,
                0.947_846_627_503_022_7e-12,
                0.135_880_130_108_924_86e-14,
                -0.348_890_393_399_948_9e-21,
            ];
            const Q: [f64; 7] = [
                1.0,
                0.084_574_623_400_189_94,
                0.002_820_929_847_262_647,
                0.468_292_921_940_894_24e-4,
                0.399_968_812_193_862_1e-6,
                0.161_809_290_887_904_48e-8,
                0.231_558_608_310_259_6e-11,
            ];
            let xs = x - 44.0;
            x * (Y + evaluate_polynomial(&P, xs) / evaluate_polynomial(&Q, xs))
        }
    }
}

/// Calculates the inverse error function `erf⁻¹(x)`.
///
/// The inverse error function returns the value `y` such that `erf(y) = x`.
/// Valid finite inputs are in the closed interval `[-1, 1]`; values outside
/// that interval return `NaN`, while the endpoints return signed infinity.
///
/// This implementation uses 64-bit piecewise rational approximations
/// for inverse `erf`/`erfc`.
///
/// # Examples
/// ```
/// use abax::{erf, erfinv};
///
/// assert_eq!(erfinv(0.0), 0.0);
/// assert!((erf(erfinv(0.5)) - 0.5).abs() < 1e-15);
/// ```
pub fn erfinv(x: f64) -> f64 {
    if x.is_nan() || !(-1.0..=1.0).contains(&x) {
        return f64::NAN;
    }
    if x == 1.0 {
        return f64::INFINITY;
    }
    if x == -1.0 {
        return f64::NEG_INFINITY;
    }
    if x == 0.0 {
        return x;
    }

    let (p, q, sign) = if x < 0.0 {
        let p = -x;
        (p, 1.0 - p, -1.0)
    } else {
        (x, 1.0 - x, 1.0)
    };

    sign * erfinv_imp(p, q)
}

#[cfg(test)]
mod tests {
    use super::erfinv;
    use crate::erf;

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

    #[test]
    fn handles_special_values() {
        assert!(erfinv(f64::NAN).is_nan());
        assert!(erfinv(-1.000_000_000_000_000_2).is_nan());
        assert!(erfinv(1.000_000_000_000_000_2).is_nan());
        assert_eq!(erfinv(1.0), f64::INFINITY);
        assert_eq!(erfinv(-1.0), f64::NEG_INFINITY);
        assert_eq!(erfinv(0.0), 0.0);
        assert_eq!(erfinv(-0.0).to_bits(), (-0.0f64).to_bits());
    }

    #[test]
    fn matches_reference_values() {
        let cases = [
            (-0.999_999, -3.458_910_737_275_498, 5e-15),
            (-0.99, -1.821_386_367_718_449_6, 3e-15),
            (-0.9, -1.163_087_153_676_674_3, 2e-15),
            (-0.5, -0.476_936_276_204_469_9, 1e-15),
            (-0.1, -0.088_855_990_494_257_69, 1e-16),
            (0.1, 0.088_855_990_494_257_69, 1e-16),
            (0.5, 0.476_936_276_204_469_9, 1e-15),
            (0.9, 1.163_087_153_676_674_3, 2e-15),
            (0.99, 1.821_386_367_718_449_6, 3e-15),
            (0.999_999, 3.458_910_737_275_498, 5e-15),
        ];

        for (x, expected, tol) in cases {
            assert_close(erfinv(x), expected, tol);
        }
    }

    #[test]
    fn round_trips_through_erf() {
        let cases = [-0.99, -0.9, -0.5, -0.125, 0.125, 0.5, 0.9, 0.99, 0.999_999];

        for x in cases {
            assert_close(erf(erfinv(x)), x, 3e-15);
        }
    }
}