/*

 * // 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());

    }

}