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;
use crate::double_double::DoubleDouble;
use crate::err::erf::{Erf, erf_accurate, erf_fast};
use crate::exponents::{EXP_REDUCE_T0, EXP_REDUCE_T1, ldexp};
use crate::rounding::CpuRoundTiesEven;
use std::hint::black_box;

static ASYMPTOTIC_POLY: [[u64; 13]; 6] = [
    [
        0x3fe20dd750429b6d,
        0x3c61a1feb75a48a8,
        0xbfd20dd750429b6c,
        0x3fdb14c2f863e403,
        0xbff0ecf9db3af35d,
        0x400d9eb53ca6eeed,
        0xc030a945830d95c8,
        0x4056e8a963e2f1f5,
        0xc0829b7ccc8f396f,
        0x40b15e716e83c27e,
        0xc0e1cfdcfbcaf22a,
        0x4111986cc7a7e8fe,
        0xc1371f7540590a91,
    ],
    [
        0x3fe20dd750429ae7,
        0x3c863da89e801fd4,
        0xbfd20dd750400795,
        0x3fdb14c2f57c490c,
        0xbff0ecf95c8c9014,
        0x400d9e981f2321ef,
        0xc030a81482de1506,
        0x4056d662420a604b,
        0xc08233c96fff7772,
        0x40af5d62018d3e37,
        0xc0d9ae55e9554450,
        0x410052901e10d139,
        0xc1166465df1385f0,
    ],
    [
        0x3fe20dd75041e3fc,
        0xbc7c9b491c4920fc,
        0xbfd20dd74e5f1526,
        0x3fdb14c1d35a40e0,
        0xbff0ecdecd30e86b,
        0x400d9b4e7f725263,
        0xc030958b5ca8fb39,
        0x40563e3179bf609c,
        0xc0806bbd1cd2d0fd,
        0x40a7b66eb6d1d2f2,
        0xc0cce5a4b1afab75,
        0x40e8b5c6ae6f773c,
        0xc0f5475860326f86,
    ],
    [
        0x3fe20dd75025cfe9,
        0x3c55a92eef32fb20,
        0xbfd20dd71eb9d4e7,
        0x3fdb14af4c25db28,
        0xbff0ebc78a22b3d8,
        0x400d85287a0b3399,
        0xc03045f751e5ca1d,
        0x4054a0d87ddea589,
        0xc07ac6a0981d1eee,
        0x409f44822f567956,
        0xc0bcba372de71349,
        0x40d1a4a19f550ca4,
        0xc0d52a580455ed79,
    ],
    [
        0x3fe20dd74eb31d84,
        0xbc439c4054b7c090,
        0xbfd20dd561af98c4,
        0x3fdb1435165d9df1,
        0xbff0e6b60308e940,
        0x400d3ce30c140882,
        0xc02f2083e404c299,
        0x40520f113d89b42a,
        0xc0741433ebd89f19,
        0x4092f35b6a3154f6,
        0xc0ab020a4313cf3b,
        0x40b90f07e92da7ee,
        0xc0b6565e1d7665c3,
    ],
    [
        0x3fe20dd744b3517b,
        0xbc6f77ab25e01ab4,
        0xbfd20dcc62ec4024,
        0x3fdb125bfa4f66c1,
        0xbff0d80e65381970,
        0x400ca11fbcfa65b2,
        0xc02cd9eaffb88315,
        0x404e010db42e0da7,
        0xc06c5c85250ef6a3,
        0x4085e118d9c1eeaf,
        0xc098d74be13d3d30,
        0x40a211b1b2b5ac83,
        0xc09900be759fc663,
    ],
];

static ASYMPTOTIC_POLY_ACCURATE: [[u64; 30]; 10] = [
    [
        0x3fe20dd750429b6d,
        0x3c61ae3a912b08f0,
        0xbfd20dd750429b6d,
        0xbc51ae34c0606d68,
        0x3fdb14c2f863e924,
        0xbc796c0f4c848fc8,
        0xbff0ecf9db3e71b6,
        0x3c645d756bd288b0,
        0x400d9eb53fad4672,
        0xbcac61629de9adf2,
        0xc030a945f3d147ea,
        0x3cb8fec5ad7ece20,
        0x4056e8c02f27ca6d,
        0xc0829d1c21c363e0,
        0x40b17349b70be627,
        0xc0e28a6bb4686182,
        0x411602d1662523ca,
        0xc14ccae7625c4111,
        0x4184237d064f6e0d,
        0xc1bb1e5466ca3a2f,
        0x41e90ae06a0f6cc1,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
    ],
    [
        0x3fe20dd750429b6d,
        0x3c61adaa62435c10,
        0xbfd20dd750429b6d,
        0xbc441516126827c8,
        0x3fdb14c2f863e90b,
        0x3c7a535780ba5ed4,
        0xbff0ecf9db3e65d6,
        0xbc9089edde27ad07,
        0x400d9eb53fa52f20,
        0xbcabc9737e9464ac,
        0xc030a945f2cd7621,
        0xbcc589f28b700332,
        0x4056e8bffd7e194e,
        0xc0829d18716876e2,
        0x40b17312abe18250,
        0xc0e287e73592805c,
        0x4115ebf7394a39c1,
        0xc14c2f14d46d0cf9,
        0x4182af3d256f955e,
        0xc1b7041659ebd7aa,
        0x41e6039c232e2f71,
        0xc2070ca15c5a07cb,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
    ],
    [
        0x3fe20dd750429b6d,
        0x3c5d3c35b5d37410,
        0xbfd20dd750429b56,
        0xbc7c028415f6f81b,
        0x3fdb14c2f863c1cf,
        0x3c51bb0de6470dbc,
        0xbff0ecf9db33c363,
        0x3c80f8068459eb16,
        0x400d9eb53b9ce57b,
        0x3ca20cce33e7d84a,
        0xc030a945aa2ec4fa,
        0xbcdf6e0fcd7c6030,
        0x4056e8b824d2bfaa,
        0xc0829cc372a6d0b0,
        0x40b1703a99ddd429,
        0xc0e2749f9a267cc6,
        0x4115856a17271849,
        0xc14a8bcb4ba9753f,
        0x418035dcce882940,
        0xc1b1e5d8c5e6e043,
        0x41dfe3b4f365386e,
        0xc20398fdef2b98fe,
        0x42184234d4f4ea12,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
    ],
    [
        0x3fe20dd750429b6a,
        0x3c8ae622b765e9fd,
        0xbfd20dd750428f0e,
        0x3c703c6c67d69513,
        0x3fdb14c2f8563e8e,
        0x3c6766a6bd7aa89c,
        0xbff0ecf9d8dedd48,
        0x3c90af52e90336e3,
        0x400d9eb4aad086fe,
        0x3ca640d371d54a19,
        0xc030a93f1d01cfe0,
        0xbcc68dbd8d9c522c,
        0x4056e842e9fd5898,
        0xc08299886ef1fb80,
        0x40b15e0f0162c9a0,
        0xc0e222dbc6b04cd8,
        0x411460268db1ebdf,
        0xc1474f53ce065fb3,
        0x417961ca8553f870,
        0xc1a8788395d13798,
        0x41d35e37b25d0e81,
        0xc1f707b7457c8f5e,
        0x4211ff852df1c023,
        0xc21b75d0ec56e2cd,
        0,
        0,
        0,
        0,
        0,
        0,
    ],
    [
        0x3fe20dd750429a8f,
        0xbc766d8dda59bcea,
        0xbfd20dd7503fdbab,
        0x3c6707bdffc2b3fe,
        0x3fdb14c2f6526025,
        0xbc27fa4bb9541140,
        0xbff0ecf99c417d45,
        0xbc9748645ef7af94,
        0x400d9eaa9c712a7d,
        0x3ca79e478994ebb4,
        0xc030a8ef11fbf141,
        0x3cbb5c72d69f8954,
        0x4056e4653e0455b1,
        0xc08286909448e6cf,
        0x40b113424ce76821,
        0xc0e1346d859e76de,
        0x4111f9f6cf2293bf,
        0xc14258e6e3b337db,
        0x41714029ecd465fb,
        0xc19c530df5337a6f,
        0x41c34bc4bbccd336,
        0xc1e4a37c52641688,
        0x420019707cec2974,
        0xc21031fe736ea169,
        0x420f6b3003de3ddf,
        0,
        0,
        0,
        0,
        0,
    ],
    [
        0x3fe20dd75042756b,
        0x3c84ad9178b56910,
        0xbfd20dd74feda9e8,
        0xbc78141c70bbc8d6,
        0x3fdb14c2cb128467,
        0xbc709aebaa106821,
        0xbff0ecf603921a0b,
        0x3c97d3cb5bceaf0b,
        0x400d9e3e1751ca59,
        0x3c76622ae5642670,
        0xc030a686af57f547,
        0x3cc083b320aff6b6,
        0x4056cf0b6c027326,
        0xc0823afcb69443d3,
        0x40b03ab450d9f1b9,
        0xc0de74cdb76bcab4,
        0x410c671b60e607f1,
        0xc138f1376d324ce4,
        0x4163b64276234676,
        0xc18aff0ce13c5a8e,
        0x41aef20247251e87,
        0xc1cc9f5662f721f6,
        0x41e4687858e185e1,
        0xc1f4fa507be073c2,
        0x41fb99ac35ee4acc,
        0xc1f16cb585ee3fa9,
        0,
        0,
        0,
        0,
    ],
    [
        0x3fe20dd7503e730d,
        0x3c84e524a098a467,
        0xbfd20dd7498fa6b2,
        0x3c260a4e27751c80,
        0x3fdb14c061bd2a0c,
        0x3c695a8f847d2fc2,
        0xbff0ecd0f11b8c7d,
        0xbc94126deea76061,
        0x400d9b1344463548,
        0x3cafe09a4eca9b0e,
        0xc030996ea52a87ed,
        0xbca924f920db26c0,
        0x40567a2264b556b0,
        0xc0815dfc2c86b6b5,
        0x40accc291b62efe4,
        0xc0d81375a78e746a,
        0x41033a6f15546329,
        0xc12c1e9dc1216010,
        0x4152397ea3d43fda,
        0xc174661e5b2ea512,
        0x4193412367ca5d45,
        0xc1ade56b9d7f37c4,
        0x41c2851d9722146d,
        0xc1d19027baf0c3fe,
        0x41d7e7b8b6ab58ac,
        0xc1d4c446d56aaf22,
        0x41c1492190400505,
        0,
        0,
        0,
    ],
    [
        0x3fe20dd74ff10852,
        0x3c8a32f26deff875,
        0xbfd20dd6f06c491c,
        0x3c770c16e1793358,
        0x3fdb14a7d5e7fd4a,
        0x3c7479998b54db5b,
        0xbff0ebbdb3889c5f,
        0xbc759b853e11369c,
        0x400d89dd249d7ef8,
        0xbc84b5edf0c8c314,
        0xc0306526fb386114,
        0xbc840d04eed7c7e0,
        0x40557ff657e429ce,
        0xc07ef63e90d38630,
        0x40a6d4f34c4ea3da,
        0xc0d04542b9e36a54,
        0x40f577bf19097738,
        0xc119702fe47c736d,
        0x413a7ae12b54fdc6,
        0xc157ca3f0f7c4fa9,
        0x417225d983963cbf,
        0xc1871a6eac612f9e,
        0x4198086324225e1e,
        0xc1a3de68670a7716,
        0x41a91674de4dcbe9,
        0xc1a6b44cc15b76c2,
        0x419a36dae0f30d80,
        0xc17cffc1747ea3dc,
        0,
        0,
    ],
    [
        0x3fe20dd74ba8f300,
        0xbc59dd256871d210,
        0xbfd20dd3593675bc,
        0x3c7ec0e7ffa91ad9,
        0x3fdb13eef86a077a,
        0xbc74fb5d78d411b8,
        0xbff0e5cf52a11f3a,
        0xbc851f36c779dc8c,
        0x400d4417a08b39d5,
        0x3c91be9fb5956638,
        0xc02f91b9f6ce80c3,
        0xbccc9c99dd42829c,
        0x405356439f45bb43,
        0xc078c0ca12819b48,
        0x409efcad2ecd6671,
        0xc0c21b0af6fc1039,
        0x40e327d215ee30c9,
        0xc101fabda96167b0,
        0x411d82e4373b315d,
        0xc134ed9e2ff591e9,
        0x41495c85dcd8eab5,
        0xc159f016f0a3d62a,
        0x41660e89d918b96f,
        0xc16e97be202cba64,
        0x4170d8a081619793,
        0xc16c5422b4fcfc65,
        0x4161131a9dc6aed1,
        0xc14a457d9dced257,
        0x4123605e980e8b86,
        0,
    ],
    [
        0x3fe20dd7319d4d25,
        0x3c82b02992c3b7ab,
        0xbfd20dc29c13ab1b,
        0xbc7d78d79b4ad767,
        0x3fdb115a57b5ab13,
        0xbc6aa8c45be0aa2e,
        0xbff0d58ec437efd7,
        0xbc5994f00a15e850,
        0x400cb1742e229f23,
        0xbca8000471d54399,
        0xc02d99a5edf7b946,
        0xbcbaf76ed7e35cde,
        0x4050a8b71058eb28,
        0xc072d88289da5bfc,
        0x40943ddf24168edb,
        0xc0b3e9dfc38b6d1a,
        0x40d18d4df97ab3df,
        0xc0eb550fc62dcab5,
        0x41029cb71f116ed1,
        0xc115fc9cc4e854e3,
        0x41265915fd0567b1,
        0xc1335eb5fca0e46d,
        0x413c5261ecc0d789,
        0xc14138932dc4eafc,
        0x414117d4eb18facd,
        0xc13af96163e35eca,
        0x4130454a3a63c766,
        0xc11c2ebc1d39b44a,
        0x40ff3327698e0e6b,
        0xc0d094febc3dff35,
    ],
];

// Approximation for the fast path of exp(z) for z=zh+zl,
// with |z| < 0.000130273 < 2^-12.88 and |zl| < 2^-42.6
// (assuming x^y does not overflow or underflow)
#[inline]
fn q_1(z_dd: DoubleDouble) -> DoubleDouble {
    const C: [u64; 5] = [
        0x3ff0000000000000,
        0x3ff0000000000000,
        0x3fe0000000000000,
        0x3fc5555555995d37,
        0x3fa55555558489dc,
    ];
    let z = z_dd.to_f64();
    let mut q = dd_fmla(f64::from_bits(C[4]), z_dd.hi, f64::from_bits(C[3]));

    q = dd_fmla(q, z, f64::from_bits(C[2]));

    let mut v = DoubleDouble::from_exact_add(f64::from_bits(C[1]), q * z);
    v = DoubleDouble::quick_mult(z_dd, v);
    DoubleDouble::f64_add(f64::from_bits(C[0]), v)
}

#[inline]
fn exp_1(x: DoubleDouble) -> DoubleDouble {
    const INVLOG2: f64 = f64::from_bits(0x40b71547652b82fe); /* |INVLOG2-2^12/log(2)| < 2^-43.4 */
    let k = (x.hi * INVLOG2).cpu_round_ties_even();

    const LOG2_DD: DoubleDouble = DoubleDouble::new(
        f64::from_bits(0x3bbabc9e3b39803f),
        f64::from_bits(0x3f262e42fefa39ef),
    );
    let k_dd = DoubleDouble::quick_f64_mult(k, LOG2_DD);
    let mut y_dd = DoubleDouble::from_exact_add(x.hi - k_dd.hi, x.lo);
    y_dd.lo -= k_dd.lo;

    let ki: i64 = k as i64; /* Note: k is an integer, this is just a conversion. */
    let mi = (ki >> 12).wrapping_add(0x3ff);
    let i2: i64 = (ki >> 6) & 0x3f;
    let i1: i64 = ki & 0x3f;
    let t1 = DoubleDouble::new(
        f64::from_bits(EXP_REDUCE_T0[i2 as usize].0),
        f64::from_bits(EXP_REDUCE_T0[i2 as usize].1),
    );
    let t2 = DoubleDouble::new(
        f64::from_bits(EXP_REDUCE_T1[i1 as usize].0),
        f64::from_bits(EXP_REDUCE_T1[i1 as usize].1),
    );
    let mut v = DoubleDouble::quick_mult(t2, t1);
    let q = q_1(y_dd);
    v = DoubleDouble::quick_mult(v, q);

    let scale = f64::from_bits((mi as u64) << 52);
    v.hi *= scale;
    v.lo *= scale;
    v
}

struct Exp {
    e: i32,
    result: DoubleDouble,
}

fn exp_accurate(x_dd: DoubleDouble) -> Exp {
    static E2: [u64; 28] = [
        0x3ff0000000000000,
        0xb960000000000000,
        0x3ff0000000000000,
        0xb9be200000000000,
        0x3fe0000000000000,
        0x3a03c00000000000,
        0x3fc5555555555555,
        0x3c655555555c78d9,
        0x3fa5555555555555,
        0x3c455555545616e2,
        0x3f81111111111111,
        0x3c011110121fc314,
        0x3f56c16c16c16c17,
        0xbbef49e06ee3a56e,
        0x3f2a01a01a01a01a,
        0x3b6b053e1eeab9c0,
        0x3efa01a01a01a01a,
        0x3ec71de3a556c733,
        0x3e927e4fb7789f66,
        0x3e5ae64567f54abe,
        0x3e21eed8eff8958b,
        0x3de6124613837216,
        0x3da93974aaf26a57,
        0x3d6ae7f4fd6d0bd9,
        0x3d2ae7e982620b25,
        0x3ce94e4ca59460d8,
        0x3ca69a2a4b7ef36d,
        0x3c6abfe1602308c9,
    ];
    const LOG2INV: f64 = f64::from_bits(0x3ff71547652b82fe);
    let k: i32 = unsafe {
        (x_dd.hi * LOG2INV)
            .cpu_round_ties_even()
            .to_int_unchecked::<i32>()
    };

    const LOG2_H: f64 = f64::from_bits(0x3fe62e42fefa39ef);
    /* we approximate LOG2Lacc ~ log(2) - LOG2H with 38 bits, so that
    k*LOG2Lacc is exact (k has at most 11 bits) */
    const LOG2_L: f64 = f64::from_bits(0x3c7abc9e3b398000);
    const LOG2_TINY: f64 = f64::from_bits(0x398f97b57a079a19);
    let yh = dd_fmla(-k as f64, LOG2_H, x_dd.hi);
    /* since |xh+xl| >= 2.92 we have |k| >= 4;
    (|k|-1/2)*log(2) <= |x| <= (|k|+1/2)*log(2) thus
    1-1/(2|k|) <= |x/(k*log(2))| <= 1+1/(2|k|) thus by Sterbenz theorem
    yh is exact too */
    let mut t = DoubleDouble::from_full_exact_add(-k as f64 * LOG2_L, x_dd.lo);
    let mut y_dd = DoubleDouble::from_exact_add(yh, t.hi);
    y_dd.lo = dd_fmla(-k as f64, LOG2_TINY, y_dd.lo + t.lo);
    /* now yh+yl approximates xh + xl - k*log(2), and we approximate p(yh+yl)
    in h + l */
    /* Since |xh| <= 742, we assume |xl| <= ulp(742) = 2^-43. Then since
       |k| <= round(742/log(2)) = 1070, |yl| <= 1070*LOG2L + 2^-42 < 2^-42.7.
       Since |yh| <= log(2)/2, the contribution of yl is negligible as long
       as |i*p[i]*yh^(i-1)*yl| < 2^-104, which holds for i >= 16.
       Thus for coefficients of degree 16 or more, we don't take yl into account.
    */
    let mut h = f64::from_bits(E2[19 + 8]); // degree 19
    for a in (16..=18).rev() {
        h = dd_fmla(h, y_dd.hi, f64::from_bits(E2[a + 8])); // degree i
    }
    /* degree 15: h*(yh+yl)+E2[15 + 8] */
    t = DoubleDouble::from_exact_mult(h, y_dd.hi);
    t.lo = dd_fmla(h, y_dd.lo, t.lo);
    let mut v = DoubleDouble::from_exact_add(f64::from_bits(E2[15 + 8]), t.hi);
    v.lo += t.lo;
    for a in (8..=14).rev() {
        /* degree i: (h+l)*(yh+yl)+E2[i+8] */
        t = DoubleDouble::quick_mult(v, y_dd);
        v = DoubleDouble::from_exact_add(f64::from_bits(E2[a + 8]), t.hi);
        v.lo += t.lo;
    }
    for a in (0..=7).rev() {
        /* degree i: (h+l)*(yh+yl)+E2[2i]+E2[2i+1] */
        t = DoubleDouble::quick_mult(v, y_dd);
        v = DoubleDouble::from_exact_add(f64::from_bits(E2[2 * a]), t.hi);
        v.lo += t.lo + f64::from_bits(E2[2 * a + 1]);
    }

    Exp { e: k, result: v }
}

#[cold]
fn erfc_asympt_accurate(x: f64) -> f64 {
    /* subnormal exceptions */
    if x == f64::from_bits(0x403a8f7bfbd15495) {
        return dd_fmla(
            f64::from_bits(0x0000000000000001),
            -0.25,
            f64::from_bits(0x000667bd620fd95b),
        );
    }
    let u_dd = DoubleDouble::from_exact_mult(x, x);
    let exp_result = exp_accurate(DoubleDouble::new(-u_dd.lo, -u_dd.hi));

    /* compute 1/x as double-double */
    let yh = 1.0 / x;
    /* Newton's iteration for 1/x is y -> y + y*(1-x*y) */
    let yl = yh * dd_fmla(-x, yh, 1.0);
    // yh+yl approximates 1/x
    static THRESHOLD: [u64; 10] = [
        0x3fb4500000000000,
        0x3fbe000000000000,
        0x3fc3f00000000000,
        0x3fc9500000000000,
        0x3fcf500000000000,
        0x3fd3100000000000,
        0x3fd7100000000000,
        0x3fdbc00000000000,
        0x3fe0b00000000000,
        0x3fe3000000000000,
    ];
    let mut i = 0usize;
    while i < THRESHOLD.len() && yh > f64::from_bits(THRESHOLD[i]) {
        i += 1;
    }
    let p = ASYMPTOTIC_POLY_ACCURATE[i];
    let mut u_dd = DoubleDouble::from_exact_mult(yh, yh);
    u_dd.lo = dd_fmla(2.0 * yh, yl, u_dd.lo);
    /* the polynomial p has degree 29+2i, and its coefficient of largest
    degree is p[14+6+i] */
    let mut z_dd = DoubleDouble::new(0., f64::from_bits(p[14 + 6 + i]));
    for a in (13..=27 + 2 * i).rev().step_by(2) {
        /* degree j: (zh+zl)*(uh+ul)+p[(j-1)/2+6]] */
        let v = DoubleDouble::quick_mult(z_dd, u_dd);
        z_dd = DoubleDouble::from_full_exact_add(f64::from_bits(p[(a - 1) / 2 + 6]), v.hi);
        z_dd.lo += v.lo;
    }
    for a in (1..=11).rev().step_by(2) {
        let v = DoubleDouble::quick_mult(z_dd, u_dd);
        z_dd = DoubleDouble::from_full_exact_add(f64::from_bits(p[a - 1]), v.hi);
        z_dd.lo += v.lo + f64::from_bits(p[a]);
    }

    /* multiply by yh+yl */
    u_dd = DoubleDouble::quick_mult(z_dd, DoubleDouble::new(yl, yh));
    /* now uh+ul approximates p(1/x), i.e., erfc(x)*exp(x^2) */
    /* now multiply (uh+ul)*(eh+el), after normalizing uh+ul to reduce the
    number of exceptional cases */
    u_dd = DoubleDouble::from_exact_add(u_dd.hi, u_dd.lo);
    let v = DoubleDouble::quick_mult(u_dd, exp_result.result);
    /* multiply by 2^e */
    /* multiply by 2^e */
    let mut res = ldexp(v.to_f64(), exp_result.e);
    if res < f64::from_bits(0x0010000000000000) {
        /* for erfc(x) in the subnormal range, we have to perform a special
        rounding */
        let mut corr = v.hi - ldexp(res, -exp_result.e);
        corr += v.lo;
        /* add corr*2^e */
        res += ldexp(corr, exp_result.e);
    }
    res
}

#[cold]
fn erfc_accurate(x: f64) -> f64 {
    if x < 0. {
        let mut v_dd = erf_accurate(-x);
        let t = DoubleDouble::from_exact_add(1.0, v_dd.hi);
        v_dd.hi = t.hi;
        v_dd.lo += t.lo;
        return v_dd.to_f64();
    } else if x <= f64::from_bits(0x3ffb59ffb450828c) {
        // erfc(x) >= 2^-6
        let mut v_dd = erf_accurate(x);
        let t = DoubleDouble::from_exact_add(1.0, -v_dd.hi);
        v_dd.hi = t.hi;
        v_dd.lo = t.lo - v_dd.lo;
        return v_dd.to_f64();
    }
    // now 0x1.b59ffb450828cp+0 < x < 0x1.b39dc41e48bfdp+4
    erfc_asympt_accurate(x)
}

/* Fast path for 0x1.713786d9c7c09p+1 < x < 0x1.b39dc41e48bfdp+4,
using the asymptotic formula erfc(x) = exp(-x^2) * p(1/x)*/
fn erfc_asympt_fast(x: f64) -> Erf {
    /* for x >= 0x1.9db1bb14e15cap+4, erfc(x) < 2^-970, and we might encounter
    underflow issues in the computation of l, thus we delegate this case
    to the accurate path */
    if x >= f64::from_bits(0x4039db1bb14e15ca) {
        return Erf {
            err: 1.0,
            result: DoubleDouble::default(),
        };
    }

    let mut u = DoubleDouble::from_exact_mult(x, x);
    let e_dd = exp_1(DoubleDouble::new(-u.lo, -u.hi));

    /* the assumptions from exp_1 are satisfied:
    * a_mul ensures |ul| <= ulp(uh), thus |ul/uh| <= 2^-52
    * since |x| < 0x1.9db1bb14e15cap+4 we have
      |ul| < ulp(0x1.9db1bb14e15cap+4^2) = 2^-43 */
    /* eh+el approximates exp(-x^2) with maximal relative error 2^-74.139 */

    /* compute 1/x as double-double */
    let yh = 1.0 / x;
    /* Assume 1 <= x < 2, then 0.5 <= yh <= 1,
    and yh = 1/x + eps with |eps| <= 2^-53. */
    /* Newton's iteration for 1/x is y -> y + y*(1-x*y) */
    let yl = yh * dd_fmla(-x, yh, 1.0);
    /* x*yh-1 = x*(1/x+eps)-1 = x*eps
       with |x*eps| <= 2^-52, thus the error on the FMA is bounded by
       ulp(2^-52.1) = 2^-105.
       Now |yl| <= |yh| * 2^-52 <= 2^-52, thus the rounding error on
       yh * __builtin_fma (-x, yh, 1.0) is bounded also by ulp(2^-52.1) = 2^-105.
       From [6], Lemma 3.7, if yl was computed exactly, then yh+yl would differ
       from 1/x by at most yh^2/theta^3*(1/x-yh)^2 for some theta in [yh,1/x]
       or [1/x,yh].
       Since yh, 1/x <= 1, this gives eps^2 <= 2^-106.
       Adding the rounding errors, we have:
       |yh + yl - 1/x| <= 2^-105 + 2^-105 + 2^-106 < 2^-103.67.
       For the relative error, since |yh| >= 1/2, this gives:
       |yh + yl - 1/x| < 2^-102.67 * |yh+yl|
    */
    const THRESHOLD: [u64; 6] = [
        0x3fbd500000000000,
        0x3fc59da6ca291ba6,
        0x3fcbc00000000000,
        0x3fd0c00000000000,
        0x3fd3800000000000,
        0x3fd6300000000000,
    ];
    let mut i = 0usize;
    while i < THRESHOLD.len() && yh > f64::from_bits(THRESHOLD[i]) {
        i += 1;
    }
    let p = ASYMPTOTIC_POLY[i];
    u = DoubleDouble::from_exact_mult(yh, yh);
    /* Since |yh| <= 1, we have |uh| <= 1 and |ul| <= 2^-53. */
    u.lo = dd_fmla(2.0 * yh, yl, u.lo);
    /* uh+ul approximates (yh+yl)^2, with absolute error bounded by
       ulp(ul) + yl^2, where ulp(ul) is the maximal rounding error in
       the FMA, and yl^2 is the neglected term.
       Since |ul| <= 2^-53, ulp(ul) <= 2^-105, and since |yl| <= 2^-52,
       this yields |uh + ul - yh^2| <= 2^-105 + 2^-104 < 2^-103.41.
       For the relative error, since |(yh+yl)^2| >= 1/4:
       |uh + ul - yh^2| < 2^-101.41 * |uh+ul|.
       And relatively to 1/x^2:
       yh + yl = 1/x * (1 + eps1)       with |eps1| < 2^-102.67
       uh + ul = (yh+yl)^2 * (1 + eps2) with |eps2| < 2^-101.41
       This yields:
       |uh + ul - 1/x| < 2^-100.90 * |uh+ul|.
    */

    /* evaluate p(uh+ul) */
    let mut zh: f64 = f64::from_bits(p[12]); // degree 23
    zh = dd_fmla(zh, u.hi, f64::from_bits(p[11])); // degree 21
    zh = dd_fmla(zh, u.hi, f64::from_bits(p[10])); // degree 19

    /* degree 17: zh*(uh+ul)+p[i] */
    let mut v = DoubleDouble::quick_f64_mult(zh, u);
    let mut z_dd = DoubleDouble::from_exact_add(f64::from_bits(p[9]), v.hi);
    z_dd.lo += v.lo;

    for a in (3..=15).rev().step_by(2) {
        v = DoubleDouble::quick_mult(z_dd, u);
        z_dd = DoubleDouble::from_exact_add(f64::from_bits(p[((a + 1) / 2) as usize]), v.hi);
        z_dd.lo += v.lo;
    }

    /* degree 1: (zh+zl)*(uh+ul)+p[0]+p[1] */
    v = DoubleDouble::quick_mult(z_dd, u);
    z_dd = DoubleDouble::from_exact_add(f64::from_bits(p[0]), v.hi);
    z_dd.lo += v.lo + f64::from_bits(p[1]);
    /* multiply by yh+yl */
    u = DoubleDouble::quick_mult(z_dd, DoubleDouble::new(yl, yh));
    /* now uh+ul approximates p(1/x) */
    /* now multiply (uh+ul)*(eh+el) */
    v = DoubleDouble::quick_mult(u, e_dd);

    /* Write y = 1/x.  We have the following errors:
       * the maximal mathematical error is:
         |erfc(x)*exp(x^2) - p(y)| < 2^-71.804 * |p(y)| (for i=3) thus
         |erfc(x) - exp(-x^2)*p(y)| < 2^-71.804 * |exp(-x^2)*p(y)|
       * the error in approximating exp(-x^2) by eh+el:
         |eh + el - exp(-x^2)| < 2^-74.139 * |eh + el|
       * the fact that we evaluate p on yh+yl instead of 1/x
         this error is bounded by |p'| * |yh+yl - 1/x|, where
         |yh+yl - 1/x| < 2^-102.67 * |yh+yl|, and the relative
         error is bounded by |p'/p| * |yh+yl - 1/x|.
         Since the maximal value of |p'/p| is bounded by 27.2 (for i=0),
         this yields 27.2 * 2^-102.67 < 2^-97.9
       * the rounding errors when evaluating p on yh+yl: this error is bounded
         (relatively) by 2^-67.184 (for i=5), see analyze_erfc_asympt_fast()
         in erfc.sage
       * the rounding error in (uh+ul)*(eh+el): we assume this error is bounded
         by 2^-80 (relatively)
       This yields a global relative bound of:
       (1+2^-71.804)*(1+2^-74.139)*(1+2^-97.9)*(1+2^-67.184)*(1+2^-80)-1
       < 2^-67.115
    */
    if v.hi >= f64::from_bits(0x044151b9a3fdd5c9) {
        Erf {
            err: f64::from_bits(0x3bbd900000000000) * v.hi,
            result: v,
        } /* 2^-67.115 < 0x1.d9p-68 */
    } else {
        Erf {
            result: v,
            err: f64::from_bits(0x0010000000000000),
        } // this overestimates 0x1.d9p-68 * h
    }
}

#[inline]
fn erfc_fast(x: f64) -> Erf {
    if x < 0.
    // erfc(x) = 1 - erf(x) = 1 + erf(-x)
    {
        let res = erf_fast(-x);
        /* h+l approximates erf(-x), with relative error bounded by err,
        where err <= 0x1.78p-69 */
        let err = res.err * res.result.hi; /* convert into absolute error */
        let mut t = DoubleDouble::from_exact_add(1.0, res.result.hi);
        t.lo += res.result.lo;
        // since h <= 2, the fast_two_sum() error is bounded by 2^-105*h <= 2^-104
        /* After the fast_two_sum() call, we have |t| <= ulp(h) <= ulp(2) = 2^-51
        thus assuming |l| <= 2^-51 after the cr_erf_fast() call,
        we have |t| <= 2^-50 here, thus the rounding
        error on t -= *l is bounded by ulp(2^-50) = 2^-102.
        The absolute error is thus bounded by err + 2^-104 + 2^-102
        = err + 0x1.4p-102.
        The maximal value of err here is for |x| < 0.0625, where cr_erf_fast()
        returns 0x1.78p-69, and h=1/2, yielding err = 0x1.78p-70 here.
        Adding 0x1.4p-102 is thus exact. */
        return Erf {
            err: err + f64::from_bits(0x3994000000000000),
            result: t,
        };
    } else if x <= f64::from_bits(0x400713786d9c7c09) {
        let res = erf_fast(x);
        /* h+l approximates erf(x), with relative error bounded by err,
        where err <= 0x1.78p-69 */
        let err = res.err * res.result.hi; /* convert into absolute error */
        let mut t = DoubleDouble::from_exact_add(1.0, -res.result.hi);
        t.lo -= res.result.lo;
        /* for x >= 0x1.e861fbb24c00ap-2, erf(x) >= 1/2, thus 1-h is exact
        by Sterbenz theorem, thus t = 0 in fast_two_sum(), and we have t = -l
        here, thus the absolute error is err */
        if x >= f64::from_bits(0x3fde861fbb24c00a) {
            return Erf { err, result: t };
        }
        /* for x < 0x1.e861fbb24c00ap-2, the error in fast_two_sum() is bounded
        by 2^-105*h, and since h <= 1/2, this yields 2^-106.
        After the fast_two_sum() call, we have |t| <= ulp(h) <= ulp(1/2) = 2^-53
        thus assuming |l| <= 2^-53 after the cr_erf_fast() call,
        we have |t| <= 2^-52 here, thus the rounding
        error on t -= *l is bounded by ulp(2^-52) = 2^-104.
        The absolute error is thus bounded by err + 2^-106 + 2^-104
        The maximal value of err here is for x < 0.0625, where cr_erf_fast()
        returns 0x1.78p-69, and h=1/2, yielding err = 0x1.78p-70 here.
        Adding 0x1.4p-104 is thus exact. */
        return Erf {
            err: err + f64::from_bits(0x3974000000000000),
            result: t,
        };
    }
    /* Now THRESHOLD1 < x < 0x1.b39dc41e48bfdp+4 thus erfc(x) < 0.000046. */
    /* on a i7-8700 with gcc 12.2.0, for x in [THRESHOLD1,+5.0],
    the average reciprocal throughput is about 111 cycles
    (among which 20 cycles for exp_1) */
    erfc_asympt_fast(x)
}

/// Complementary error function
///
/// Max ulp 0.5
pub fn f_erfc(x: f64) -> f64 {
    let t: u64 = x.to_bits();
    let at: u64 = t & 0x7fff_ffff_ffff_ffff;

    if t >= 0x8000000000000000u64
    // x = -NaN or x <= 0 (excluding +0)
    {
        // for x <= -0x1.7744f8f74e94bp2, erfc(x) rounds to 2 (to nearest)
        if t >= 0xc017744f8f74e94bu64
        // x = NaN or x <= -0x1.7744f8f74e94bp2
        {
            if t >= 0xfff0000000000000u64 {
                // -Inf or NaN
                if t == 0xfff0000000000000u64 {
                    return 2.0;
                } // -Inf
                return x + x; // NaN
            }
            return black_box(2.0) - black_box(f64::from_bits(0x3c90000000000000)); // rounds to 2 or below(2)
        }

        // for -9.8390953768041405e-17 <= x <= 0, erfc(x) rounds to 1 (to nearest)
        if f64::from_bits(0xbc9c5bf891b4ef6a) <= x {
            return dd_fmla(-x, f64::from_bits(0x3c90000000000000), 1.0);
        }
    } else
    // x = +NaN or x >= 0 (excluding -0)
    {
        // for x >= 0x1.b39dc41e48bfdp+4, erfc(x) < 2^-1075: rounds to 0 or 2^-1074
        if at >= 0x403b39dc41e48bfdu64
        // x = NaN or x >= 0x1.b39dc41e48bfdp+4
        {
            if at >= 0x7ff0000000000000u64 {
                // +Inf or NaN
                if at == 0x7ff0000000000000u64 {
                    return 0.0;
                } // +Inf
                return x + x; // NaN
            }
            return black_box(f64::from_bits(0x0000000000000001)) * black_box(0.25); // 0 or 2^-1074 wrt rounding
        }

        // for 0 <= x <= 0x1.c5bf891b4ef6ap-55, erfc(x) rounds to 1 (to nearest)
        if x <= f64::from_bits(0x3c8c5bf891b4ef6a) {
            return dd_fmla(-x, f64::from_bits(0x3c90000000000000), 1.0);
        }
    }

    /* now -0x1.7744f8f74e94bp+2 < x < -0x1.c5bf891b4ef6ap-54
    or 0x1.c5bf891b4ef6ap-55 < x < 0x1.b39dc41e48bfdp+4 */

    let result = erfc_fast(x);

    let left = result.result.hi + (result.result.lo - result.err);
    let right = result.result.hi + (result.result.lo + result.err);
    if left == right {
        return left;
    }
    erfc_accurate(x)
}

#[cfg(test)]
mod tests {
    use super::*;
    #[test]
    fn test_erfc() {
        assert_eq!(f_erfc(1.0), 0.15729920705028513);
        assert_eq!(f_erfc(0.5), 0.4795001221869535);
        assert_eq!(f_erfc(0.000000005), 0.9999999943581042);
        assert_eq!(f_erfc(-0.00000000000065465465423305), 1.0000000000007387);
        assert!(f_erfc(f64::NAN).is_nan());
        assert_eq!(f_erfc(f64::INFINITY), 0.0);
        assert_eq!(f_erfc(f64::NEG_INFINITY), 2.0);
    }
}