pxfm 0.1.29

/*
 * // Copyright (c) Radzivon Bartoshyk 7/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
use crate::common::{dd_fmla, dyad_fmla, f_fmla};
use crate::double_double::DoubleDouble;
use crate::err::erf_poly::{ERF_POLY, ERF_POLY_C2};
use crate::rounding::CpuFloor;
/* double-double approximation of 2/sqrt(pi) to nearest */
const TWO_OVER_SQRT_PI: DoubleDouble = DoubleDouble::new(
    f64::from_bits(0x3c71ae3a914fed80),
    f64::from_bits(0x3ff20dd750429b6d),
);

pub(crate) struct Erf {
    pub(crate) result: DoubleDouble,
    pub(crate) err: f64,
}

/* for |z| < 1/8, assuming z >= 2^-61, thus no underflow can occur */
#[cold]
fn cr_erf_accurate_tiny(x: f64) -> DoubleDouble {
    static P: [u64; 15] = [
        0x3ff20dd750429b6d,
        0x3c71ae3a914fed80,
        0xbfd812746b0379e7,
        0x3c6ee12e49ca96ba,
        0x3fbce2f21a042be2,
        0xbc52871bc0a0a0d0,
        0xbf9b82ce31288b51,
        0x3c21003accf1355c,
        0x3f7565bcd0e6a53f,
        0xbf4c02db40040cc3,
        0x3f1f9a326fa3cf50,
        0xbeef4d25e3c73ce9,
        0x3ebb9eb332b31646,
        0xbe864a4bd5eca4d7,
        0x3e6c0acc2502e94e,
    ];
    let z2 = x * x;
    let mut h = f64::from_bits(P[21 / 2 + 4]); /* degree 21 */
    for a in (12..=19).rev().step_by(2) {
        h = dd_fmla(h, z2, f64::from_bits(P[(a / 2 + 4) as usize]))
    }
    let mut l = 0.;
    for a in (8..=11).rev().step_by(2) {
        let mut t = DoubleDouble::from_exact_mult(h, x);
        t.lo = dd_fmla(l, x, t.lo);
        let mut k = DoubleDouble::from_exact_mult(t.hi, x);
        k.lo = dd_fmla(t.lo, x, k.lo);
        let p = DoubleDouble::from_exact_add(f64::from_bits(P[(a / 2 + 4) as usize]), k.hi);
        l = k.lo + p.lo;
        h = p.hi;
    }
    for a in (1..=7).rev().step_by(2) {
        let mut t = DoubleDouble::from_exact_mult(h, x);
        t.lo = dd_fmla(l, x, t.lo);
        let mut k = DoubleDouble::from_exact_mult(t.hi, x);
        k.lo = dd_fmla(t.lo, x, k.lo);
        let p = DoubleDouble::from_exact_add(f64::from_bits(P[a - 1]), k.hi);
        l = k.lo + p.lo + f64::from_bits(P[a]);
        h = p.hi;
    }
    /* multiply by z */
    let p = DoubleDouble::from_exact_mult(h, x);
    l = dd_fmla(l, x, p.lo);
    DoubleDouble::new(l, p.hi)
}

/* Assuming 0 <= z <= 0x1.7afb48dc96626p+2, put in h+l an accurate
approximation of erf(z).
Assumes z >= 2^-61, thus no underflow can occur. */
#[cold]
#[inline(never)]
pub(crate) fn erf_accurate(x: f64) -> DoubleDouble {
    if x < 0.125
    /* z < 1/8 */
    {
        return cr_erf_accurate_tiny(x);
    }

    let v = (8.0 * x).cpu_floor();
    let i: u32 = (8.0 * x) as u32;
    let z = (x - 0.0625) - 0.125 * v;
    /* now |z| <= 1/16 */
    let p = ERF_POLY_C2[(i - 1) as usize];
    let mut h = f64::from_bits(p[26]); /* degree-18 */
    for a in (11..=17).rev() {
        h = dd_fmla(h, z, f64::from_bits(p[(8 + a) as usize])); /* degree j */
    }
    let mut l: f64 = 0.;
    for a in (8..=10).rev() {
        let mut t = DoubleDouble::from_exact_mult(h, z);
        t.lo = dd_fmla(l, z, t.lo);
        let p = DoubleDouble::from_exact_add(f64::from_bits(p[(8 + a) as usize]), t.hi);
        h = p.hi;
        l = p.lo + t.lo;
    }
    for a in (0..=7).rev() {
        let mut t = DoubleDouble::from_exact_mult(h, z);
        t.lo = dd_fmla(l, z, t.lo);
        /* add p[2*j] + p[2*j+1] to th + tl: we use two_sum() instead of
        fast_two_sum because for example for i=3, the coefficient of
        degree 7 is tiny (0x1.060b78c935b8ep-13) with respect to that
        of degree 8 (0x1.678b51a9c4b0ap-7) */
        let v = DoubleDouble::from_exact_add(f64::from_bits(p[(2 * a) as usize]), t.hi);
        h = v.hi;
        l = v.lo + t.lo + f64::from_bits(p[(2 * a + 1) as usize]);
    }
    DoubleDouble::new(l, h)
}

/* Assuming 0 <= z <= 5.9215871957945065, put in h+l an approximation
of erf(z). Return err the maximal relative error:
|(h + l)/erf(z) - 1| < err*|h+l| */
#[inline]
pub(crate) fn erf_fast(x: f64) -> Erf {
    /* we split [0,5.9215871957945065] into intervals i/16 <= z < (i+1)/16,
       and for each interval, we use a minimax polynomial:
       * for i=0 (0 <= z < 1/16) we use a polynomial evaluated at zero,
         since if we evaluate in the middle 1/32, we will get bad accuracy
         for tiny z, and moreover z-1/32 might not be exact
       * for 1 <= i <= 94, we use a polynomial evaluated in the middle of
         the interval, namely i/16+1/32
    */
    if x < 0.0625
    /* z < 1/16 */
    {
        /* the following is a degree-11 minimax polynomial for erf(x) on [0,1/16]
        generated by Sollya, with double-double coefficients for degree 1 and 3,
        and double coefficients for degrees 5 to 11 (file erf0.sollya).
        The maximal relative error is 2^-68.935. */
        let z2 = DoubleDouble::from_exact_mult(x, x);
        const C: [u64; 8] = [
            0x3ff20dd750429b6d,
            0x3c71ae3a7862d9c4,
            0xbfd812746b0379e7,
            0x3c6f1a64d72722a2,
            0x3fbce2f21a042b7f,
            0xbf9b82ce31189904,
            0x3f7565bbf8a0fe0b,
            0xbf4bf9f8d2c202e4,
        ];
        let z4 = z2.hi * z2.hi;
        let c9 = dd_fmla(f64::from_bits(C[7]), z2.hi, f64::from_bits(C[6]));
        let mut c5 = dd_fmla(f64::from_bits(C[5]), z2.hi, f64::from_bits(C[4]));
        c5 = dd_fmla(c9, z4, c5);
        /* compute c0[2] + c0[3] + z2h*c5 */
        let mut t = DoubleDouble::from_exact_mult(z2.hi, c5);
        let mut v = DoubleDouble::from_exact_add(f64::from_bits(C[2]), t.hi);
        v.lo += t.lo + f64::from_bits(C[3]);
        /* compute c0[0] + c0[1] + (z2h + z2l)*(h + l) */
        t = DoubleDouble::from_exact_mult(z2.hi, v.hi);
        let h_c = v.hi;
        t.lo += dd_fmla(z2.hi, v.lo, f64::from_bits(C[1]));
        v = DoubleDouble::from_exact_add(f64::from_bits(C[0]), t.hi);
        v.lo += dd_fmla(z2.lo, h_c, t.lo);
        v = DoubleDouble::quick_mult_f64(v, x);
        return Erf {
            result: v,
            err: f64::from_bits(0x3ba7800000000000),
        }; /* err < 2.48658249618372e-21, cf Analyze0() */
    }

    let v = (16.0 * x).cpu_floor();
    let i: u32 = (16.0 * x) as u32;
    /* i/16 <= z < (i+1)/16 */
    /* For 0.0625 0 <= z <= 0x1.7afb48dc96626p+2, z - 0.03125 is exact:
    (1) either z - 0.03125 is in the same binade as z, then 0.03125 is
        an integer multiple of ulp(z), so is z - 0.03125
    (2) if z - 0.03125 is in a smaller binade, both z and 0.03125 are
        integer multiple of the ulp() of that smaller binade.
    Also, subtracting 0.0625 * v is exact. */
    let z = (x - 0.03125) - 0.0625 * v;
    /* now |z| <= 1/32 */
    let c = ERF_POLY[(i - 1) as usize];
    let z2 = z * z;
    let z4 = z2 * z2;
    /* the degree-10 coefficient is c[12] */
    let c9 = dd_fmla(f64::from_bits(c[12]), z, f64::from_bits(c[11]));
    let mut c7 = dd_fmla(f64::from_bits(c[10]), z, f64::from_bits(c[9]));
    let c5 = dd_fmla(f64::from_bits(c[8]), z, f64::from_bits(c[7]));
    /* c3h, c3l <- c[5] + z*c[6] */
    let mut c3 = DoubleDouble::from_exact_add(f64::from_bits(c[5]), z * f64::from_bits(c[6]));
    c7 = dd_fmla(c9, z2, c7);
    /* c3h, c3l <- c3h, c3l + c5*z2 */
    let p = DoubleDouble::from_exact_add(c3.hi, c5 * z2);
    c3.hi = p.hi;
    c3.lo += p.lo;
    /* c3h, c3l <- c3h, c3l + c7*z4 */
    let p = DoubleDouble::from_exact_add(c3.hi, c7 * z4);
    c3.hi = p.hi;
    c3.lo += p.lo;
    /* c2h, c2l <- c[4] + z*(c3h + c3l) */
    let mut t = DoubleDouble::from_exact_mult(z, c3.hi);
    let mut c2 = DoubleDouble::from_exact_add(f64::from_bits(c[4]), t.hi);
    c2.lo += dd_fmla(z, c3.lo, t.lo);
    /* compute c[2] + c[3] + z*(c2h + c2l) */
    t = DoubleDouble::from_exact_mult(z, c2.hi);
    let mut v = DoubleDouble::from_exact_add(f64::from_bits(c[2]), t.hi);
    v.lo += t.lo + dd_fmla(z, c2.lo, f64::from_bits(c[3]));
    /* compute c[0] + c[1] + z*(h + l) */
    t = DoubleDouble::from_exact_mult(z, v.hi);
    t.lo = dd_fmla(z, v.lo, t.lo);
    v = DoubleDouble::from_exact_add(f64::from_bits(c[0]), t.hi);
    v.lo += t.lo + f64::from_bits(c[1]);
    Erf {
        result: v,
        err: f64::from_bits(0x3ba1100000000000),
    } /* err < 1.80414390200020e-21, cf analyze_p(1)
    (larger values of i yield smaller error bounds) */
}

/// Error function
///
/// Max ULP 0.5
pub fn f_erf(x: f64) -> f64 {
    let z = f64::from_bits(x.to_bits() & 0x7fff_ffff_ffff_ffff);
    let mut t = z.to_bits();
    let ux = t;
    /* erf(x) rounds to +/-1 for RNDN for |x| > 0x4017afb48dc96626 */
    if ux > 0x4017afb48dc96626
    // |x| > 0x4017afb48dc96626
    {
        let os = f64::copysign(1.0, x);
        const MASK: u64 = 0x7ff0000000000000u64;
        if ux > MASK {
            return x + x; /* NaN */
        }
        if ux == MASK {
            return os; /* +/-Inf */
        }
        return f_fmla(-f64::from_bits(0x3c90000000000000), os, os);
    }

    /* now |x| <= 0x4017afb48dc96626 */
    if z < f64::from_bits(0x3c20000000000000) {
        /* for x=-0 the code below returns +0 which is wrong */
        if x == 0. {
            return x;
        }
        /* tiny x: erf(x) ~ 2/sqrt(pi) * x + O(x^3), where the ratio of the O(x^3)
        term to the main term is in x^2/3, thus less than 2^-123 */
        let y = TWO_OVER_SQRT_PI.hi * x; /* tentative result */
        /* scale x by 2^106 to get out the subnormal range */
        let sx = x * f64::from_bits(0x4690000000000000);
        let mut p = DoubleDouble::quick_mult_f64(TWO_OVER_SQRT_PI, sx);
        /* now compute the residual h + l - y */
        p.lo += f_fmla(-y, f64::from_bits(0x4690000000000000), p.hi); /* h-y*2^106 is exact since h and y are very close */
        let res = dyad_fmla(p.lo, f64::from_bits(0x3950000000000000), y);
        return res;
    }

    let result = erf_fast(z);
    let mut u = result.result.hi.to_bits();
    let mut v = result.result.lo.to_bits();
    t = x.to_bits();

    const SIGN_MASK: u64 = 0x8000000000000000u64;
    u ^= t & SIGN_MASK;
    v ^= t & SIGN_MASK;

    let left = f64::from_bits(u) + f_fmla(result.err, -f64::from_bits(u), f64::from_bits(v));
    let right = f64::from_bits(u) + f_fmla(result.err, f64::from_bits(u), f64::from_bits(v));

    if left == right {
        return left;
    }

    let a_results = erf_accurate(z);

    if x >= 0. {
        a_results.to_f64()
    } else {
        (-a_results.hi) + (-a_results.lo)
    }
}

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

    #[test]
    fn test_erf() {
        assert_eq!(f_erf(0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000009456563898732),
                   0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010670589695636709);
        assert_eq!(f_erf(0.), 0.);
        assert_eq!(f_erf(1.), 0.8427007929497149);
        assert_eq!(f_erf(0.49866735123), 0.5193279892991808);
        assert_eq!(f_erf(-0.49866735123), -0.5193279892991808);
        assert!(f_erf(f64::NAN).is_nan());
        assert_eq!(f_erf(f64::INFINITY), 1.0);
        assert_eq!(f_erf(f64::NEG_INFINITY), -1.0);
    }
}