pxfm 0.1.29

/*
 * // Copyright (c) Radzivon Bartoshyk 9/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::double_double::DoubleDouble;
use crate::pow_exec::exp_dd_fast;

#[inline]
fn core_erfcx(x: f64) -> DoubleDouble {
    if x <= 8. {
        // Rational approximant for erfcx generated by Wolfram Mathematica:
        // <<FunctionApproximations`
        // ClearAll["Global`*"]
        // f[x_]:=Exp[x^2]Erfc[x]
        // {err0,approx,err1}=MiniMaxApproximation[f[z],{z,{1, 8},11,11},WorkingPrecision->75,MaxIterations->100]
        // num=Numerator[approx];
        // den=Denominator[approx];
        // coeffs=CoefficientList[num,z];
        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
        // coeffs=CoefficientList[den,z];
        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
        const P: [(u64, u64); 12] = [
            (0xbc836faeb9a312bb, 0x3ff000000000ee8e),
            (0x3c91842f891bec6a, 0x4002ca20a78aaf8f),
            (0x3c7916e8a1c30681, 0x4005e955f70aed5b),
            (0x3cabad150d828d82, 0x4000646f5807ad07),
            (0xbc6f482680d66d9c, 0x3ff1449e03ed381c),
            (0xbc7188796156ae19, 0x3fdaa7e997e3b034),
            (0xbc5c8af0642761e3, 0x3fbe836282058d4a),
            (0xbc372829be2d072f, 0x3f99a2b2adc2ec05),
            (0x3c020cc8b96000ab, 0x3f6e6cc3d120a955),
            (0x3bdd138e6c136806, 0x3f3743d6735eaf13),
            (0xbb9fbd22f0675122, 0x3ef1c1d36ebe29a2),
            (0xb89093cc981c934c, 0xbc43c18bc6385c74),
        ];

        let x2 = DoubleDouble::from_exact_mult(x, x);
        let x4 = x2 * x2;
        let x8 = x4 * x4;

        let e0 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(P[1]),
            x,
            DoubleDouble::from_bit_pair(P[0]),
        );
        let e1 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(P[3]),
            x,
            DoubleDouble::from_bit_pair(P[2]),
        );
        let e2 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(P[5]),
            x,
            DoubleDouble::from_bit_pair(P[4]),
        );
        let e3 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(P[7]),
            x,
            DoubleDouble::from_bit_pair(P[6]),
        );
        let e4 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(P[9]),
            x,
            DoubleDouble::from_bit_pair(P[8]),
        );
        let e5 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(P[11]),
            x,
            DoubleDouble::from_bit_pair(P[10]),
        );

        let f0 = DoubleDouble::mul_add(x2, e1, e0);
        let f1 = DoubleDouble::mul_add(x2, e3, e2);
        let f2 = DoubleDouble::mul_add(x2, e5, e4);

        let g0 = DoubleDouble::mul_add(x4, f1, f0);

        let p_num = DoubleDouble::mul_add(x8, f2, g0);

        const Q: [(u64, u64); 12] = [
            (0x0000000000000000, 0x3ff0000000000000),
            (0xbc95d65be031374e, 0x400bd10c4fb1dbe5),
            (0x3cb2d8f661db08a0, 0x4016a649ff973199),
            (0x3ca32cbcfdc0ea93, 0x4016daab399c1ffc),
            (0xbca2982868536578, 0x400fd61ab892d14c),
            (0xbca2e29199e17fd9, 0x40001f56c4d495a3),
            (0x3c412ce623a1790a, 0x3fe852b582135164),
            (0x3c61152eaf4b0dc5, 0x3fcb760564da7cde),
            (0xbc1b57ff91d81959, 0x3fa6e146988df835),
            (0x3c17183d8445f19a, 0x3f7b06599b5e912f),
            (0xbbd0ada61b85ff98, 0x3f449e39467b73d0),
            (0xbb658d84fc735e67, 0x3eff794442532b51),
        ];

        let e0 = DoubleDouble::mul_f64_add_f64(
            DoubleDouble::from_bit_pair(Q[1]),
            x,
            f64::from_bits(0x3ff0000000000000),
        );
        let e1 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(Q[3]),
            x,
            DoubleDouble::from_bit_pair(Q[2]),
        );
        let e2 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(Q[5]),
            x,
            DoubleDouble::from_bit_pair(Q[4]),
        );
        let e3 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(Q[7]),
            x,
            DoubleDouble::from_bit_pair(Q[6]),
        );
        let e4 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(Q[9]),
            x,
            DoubleDouble::from_bit_pair(Q[8]),
        );
        let e5 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(Q[11]),
            x,
            DoubleDouble::from_bit_pair(Q[10]),
        );

        let f0 = DoubleDouble::mul_add(x2, e1, e0);
        let f1 = DoubleDouble::mul_add(x2, e3, e2);
        let f2 = DoubleDouble::mul_add(x2, e5, e4);

        let g0 = DoubleDouble::mul_add(x4, f1, f0);

        let p_den = DoubleDouble::mul_add(x8, f2, g0);
        return DoubleDouble::div(p_num, p_den);
    }

    // for large x erfcx(x) ~ 1/sqrt(pi) / x * R(1/x)
    const ONE_OVER_SQRT_PI: DoubleDouble =
        DoubleDouble::from_bit_pair((0x3c61ae3a914fed80, 0x3fe20dd750429b6d));
    let r = DoubleDouble::from_quick_recip(x);
    // Rational approximant generated by Wolfram:
    // <<FunctionApproximations`
    // ClearAll["Global`*"]
    // f[x_]:=Exp[1/x^2]Erfc[1/x]/x*Sqrt[Pi]
    // {err0,approx}=MiniMaxApproximation[f[z],{z,{2^-23,1/8},8,8},WorkingPrecision->75,MaxIterations->100]
    // num=Numerator[approx][[1]];
    // den=Denominator[approx][[1]];
    // coeffs=CoefficientList[num,z];
    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
    // coeffs=CoefficientList[den,z];
    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
    const P: [(u64, u64); 9] = [
        (0xbb1d2ee37e46a4cd, 0x3ff0000000000000),
        (0x3ca2e575a4ce3d30, 0x4001303ab00c8bac),
        (0xbccf38381e5ee521, 0x4030a97aeed54c9f),
        (0xbcc3a2842df0dd3d, 0x4036f7733c9fd2f9),
        (0xbcfeaf46506f16ed, 0x4051c5f382750553),
        (0x3ccbb9f5e11d176a, 0x404ac0081e0749e0),
        (0xbcf374f8966ae2a5, 0x4052082526d99a5c),
        (0x3cbb5530b924f224, 0x402feabbf6571c29),
        (0xbcbcdd50a3ca4776, 0x40118726e1f2d204),
    ];
    const Q: [(u64, u64); 9] = [
        (0x0000000000000000, 0x3ff0000000000000),
        (0x3ca2e4613c9e0017, 0x4001303ab00c8bac),
        (0xbcce5f17cf14e51d, 0x4031297aeed54c9f),
        (0xbcdf7e0fed176f92, 0x40380a76e7a09bb2),
        (0x3cfc57b67a2797af, 0x4053bb22e04faf3e),
        (0xbcd3e63b7410b46b, 0x404ff46317ae9483),
        (0xbce122c15db2653f, 0x405925ef8a428c36),
        (0x3ce174ebe3e52c8e, 0x4040f49acfe692e1),
        (0xbcda0e267ce6e2e6, 0x40351a07878bfbd3),
    ];
    let mut p_num = DoubleDouble::mul_add(
        DoubleDouble::from_bit_pair(P[8]),
        r,
        DoubleDouble::from_bit_pair(P[7]),
    );
    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[6]));
    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[5]));
    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[4]));
    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[3]));
    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[2]));
    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[1]));
    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[0]));

    let mut p_den = DoubleDouble::mul_add(
        DoubleDouble::from_bit_pair(Q[8]),
        r,
        DoubleDouble::from_bit_pair(Q[7]),
    );
    p_den = DoubleDouble::mul_add(p_den, r, DoubleDouble::from_bit_pair(Q[6]));
    p_den = DoubleDouble::mul_add(p_den, r, DoubleDouble::from_bit_pair(Q[5]));
    p_den = DoubleDouble::mul_add(p_den, r, DoubleDouble::from_bit_pair(Q[4]));
    p_den = DoubleDouble::mul_add(p_den, r, DoubleDouble::from_bit_pair(Q[3]));
    p_den = DoubleDouble::mul_add(p_den, r, DoubleDouble::from_bit_pair(Q[2]));
    p_den = DoubleDouble::mul_add(p_den, r, DoubleDouble::from_bit_pair(Q[1]));
    p_den = DoubleDouble::mul_add_f64(p_den, r, f64::from_bits(0x3ff0000000000000));

    let v0 = DoubleDouble::quick_mult(ONE_OVER_SQRT_PI, r);
    let v1 = DoubleDouble::div(p_num, p_den);
    DoubleDouble::quick_mult(v0, v1)
}

/// Scaled complementary error function (exp(x^2)*erfc(x))
pub fn f_erfcx(x: f64) -> f64 {
    let ux = x.to_bits().wrapping_shl(1);

    if ux >= 0x7ffu64 << 53 || ux <= 0x7960000000000000u64 {
        // x == NaN, x == inf, x == 0, |x| <= f64::EPSILON
        if x.is_nan() {
            return f64::NAN;
        }
        if x.to_bits().wrapping_shl(1) == 0 {
            return 1.;
        }
        if x.is_infinite() {
            return if x.is_sign_positive() {
                0.
            } else {
                f64::INFINITY
            };
        }

        if ux <= 0x7888f5c28f5c28f6u64 {
            // |x| <= 2.2204460492503131e-18
            return 1.;
        }

        // |x| <= f64::EPSILON
        use crate::common::f_fmla;
        const M_TWO_OVER_SQRT_PI: DoubleDouble =
            DoubleDouble::from_bit_pair((0xbc71ae3a914fed80, 0xbff20dd750429b6d));
        return f_fmla(
            M_TWO_OVER_SQRT_PI.lo,
            x,
            f_fmla(M_TWO_OVER_SQRT_PI.hi, x, 1.),
        );
    }

    if x.to_bits() >= 0xc03aa449ebc84dd6 {
        // x <= -sqrt(709.783) ~ -26.6417
        return f64::INFINITY;
    }

    let ax = x.to_bits() & 0x7fff_ffff_ffff_ffffu64;

    if ax <= 0x3ff0000000000000u64 {
        // |x| <= 1
        // Rational approximant generated by Wolfram Mathematica:
        // <<FunctionApproximations`
        // ClearAll["Global`*"]
        // f[x_]:=Exp[x^2]Erfc[x]
        // {err0,approx}=MiniMaxApproximation[f[z],{z,{-1, 1},10,10},WorkingPrecision->75,MaxIterations->100]
        // num=Numerator[approx][[1]];
        // den=Denominator[approx][[1]];
        // coeffs=CoefficientList[num,z];
        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
        // coeffs=CoefficientList[den,z];
        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
        const P: [(u64, u64); 11] = [
            (0xbb488611350b1950, 0x3ff0000000000000),
            (0xbc86ae482c7f2342, 0x3ff9c5d39e89602f),
            (0x3c6702d70b807254, 0x3ff5a4c406d6468b),
            (0x3c7fe41fc43cfed5, 0x3fe708e7f401bd0c),
            (0x3c73a4a355172c6d, 0x3fd0d9a0c1a7126c),
            (0x3c5f4c372faa270f, 0x3fb154722e30762e),
            (0xbc04c0227976379e, 0x3f88ecebb62ce646),
            (0xbbdc9ea151b9eb33, 0x3f580ea84143877b),
            (0xbb6dae7001a91491, 0x3f1c3c5f95579b0a),
            (0x3b6aca5e82c52897, 0x3ecea4db51968d9e),
            (0x3a41c4edd175d2af, 0x3dbc0dccea7fc8ed),
        ];

        let x2 = DoubleDouble::from_exact_mult(x, x);
        let x4 = x2 * x2;
        let x8 = x4 * x4;

        let q0 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(P[1]),
            x,
            DoubleDouble::from_bit_pair(P[0]),
        );
        let q1 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(P[3]),
            x,
            DoubleDouble::from_bit_pair(P[2]),
        );
        let q2 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(P[5]),
            x,
            DoubleDouble::from_bit_pair(P[4]),
        );
        let q3 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(P[7]),
            x,
            DoubleDouble::from_bit_pair(P[6]),
        );
        let q4 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(P[9]),
            x,
            DoubleDouble::from_bit_pair(P[8]),
        );

        let r0 = DoubleDouble::mul_add(x2, q1, q0);
        let r1 = DoubleDouble::mul_add(x2, q3, q2);

        let s0 = DoubleDouble::mul_add(x4, r1, r0);
        let s1 = DoubleDouble::mul_add(x2, DoubleDouble::from_bit_pair(P[10]), q4);
        let p_num = DoubleDouble::mul_add(x8, s1, s0);

        const Q: [(u64, u64); 11] = [
            (0x0000000000000000, 0x3ff0000000000000),
            (0xbc7bae414cad99c8, 0x4005e9d57765fdce),
            (0x3c8fa553bed15758, 0x400b8c670b3fbcda),
            (0x3ca6c7ad610f1019, 0x4004f2ca59958153),
            (0x3c87787f336cc4e6, 0x3ff55c267090315a),
            (0xbc6ef55d4b2c4150, 0x3fde8b84b64b6f4e),
            (0x3c570d63c94be3a3, 0x3fbf0d5e36017482),
            (0x3c1882a745ef572e, 0x3f962f73633506c1),
            (0xbc0850bb6fc82764, 0x3f65593e0dc46acd),
            (0xbbb9dc0097d7d776, 0x3f290545603e2f94),
            (0xbb776e5781e3889d, 0x3edb29c49d18cf89),
        ];

        let q0 = DoubleDouble::mul_f64_add_f64(
            DoubleDouble::from_bit_pair(Q[1]),
            x,
            f64::from_bits(0x3ff0000000000000),
        );
        let q1 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(Q[3]),
            x,
            DoubleDouble::from_bit_pair(Q[2]),
        );
        let q2 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(Q[5]),
            x,
            DoubleDouble::from_bit_pair(Q[4]),
        );
        let q3 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(Q[7]),
            x,
            DoubleDouble::from_bit_pair(Q[6]),
        );
        let q4 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(Q[9]),
            x,
            DoubleDouble::from_bit_pair(Q[8]),
        );

        let r0 = DoubleDouble::mul_add(x2, q1, q0);
        let r1 = DoubleDouble::mul_add(x2, q3, q2);

        let s0 = DoubleDouble::mul_add(x4, r1, r0);
        let s1 = DoubleDouble::mul_add(x2, DoubleDouble::from_bit_pair(Q[10]), q4);
        let p_den = DoubleDouble::mul_add(x8, s1, s0);

        let v = DoubleDouble::div(p_num, p_den);
        return v.to_f64();
    }

    let mut erfcx_abs_x = core_erfcx(f64::from_bits(ax));
    if x < 0. {
        // exp(x^2)erfc(-x) = 2*exp(x^2) - erfcx(|x|)
        erfcx_abs_x = DoubleDouble::from_exact_add(erfcx_abs_x.hi, erfcx_abs_x.lo);
        let d2x = DoubleDouble::from_exact_mult(x, x);
        let expd2x = exp_dd_fast(d2x);
        return DoubleDouble::mul_f64_add(expd2x, 2., -erfcx_abs_x).to_f64();
    }
    erfcx_abs_x.to_f64()
}

#[cfg(test)]
mod tests {
    use crate::f_erfcx;

    #[test]
    fn test_erfcx() {
        assert_eq!(f_erfcx(2.2204460492503131e-18), 1.0);
        assert_eq!(f_erfcx(-2.2204460492503131e-18), 1.0);
        assert_eq!(f_erfcx(-f64::EPSILON), 1.0000000000000002);
        assert_eq!(f_erfcx(f64::EPSILON), 0.9999999999999998);
        assert_eq!(f_erfcx(-173.), f64::INFINITY);
        assert_eq!(f_erfcx(-9.4324165432), 8.718049147018359e38);
        assert_eq!(f_erfcx(9.4324165432), 0.059483265496416374);
        assert_eq!(f_erfcx(-1.32432512125), 11.200579112797806);
        assert_eq!(f_erfcx(1.32432512125), 0.3528722004785406);
        assert_eq!(f_erfcx(-0.532431235), 2.0560589406595384);
        assert_eq!(f_erfcx(0.532431235), 0.5994337293294584);
        assert_eq!(f_erfcx(1e-26), 1.0);
        assert_eq!(f_erfcx(-0.500000000023073), 1.952360489253639);
        assert_eq!(f_erfcx(-175.), f64::INFINITY);
        assert_eq!(f_erfcx(f64::INFINITY), 0.);
        assert_eq!(f_erfcx(f64::NEG_INFINITY), f64::INFINITY);
        assert!(f_erfcx(f64::NAN).is_nan());
    }
}