/*

 * // 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::bessel::bessel_exp::i0_exp_accurate;

use crate::bessel::i0::{

    bessel_rsqrt_hard, eval_small_hard_3p6_to_7p5, i0_0_to_3p6_dd, i0_0_to_3p6_hard,

    i0_3p6_to_7p5_dd,

};

use crate::bessel::i0_exp;

use crate::common::f_fmla;

use crate::double_double::DoubleDouble;

use crate::dyadic_float::{DyadicFloat128, DyadicSign};

/// Modified exponentially scaled Bessel of the first kind of order 0

///

/// Computes exp(-|x|)*I0(x)

pub fn f_i0e(x: f64) -> f64 {

    let ux = x.to_bits().wrapping_shl(1);

    if ux >= 0x7ffu64 << 53 || ux <= 0x7960000000000000u64 {

        // |x| <= f64::EPSILON, |x| == inf, x == NaN

        if ux <= 0x760af31dc4611874u64 {

            // |x| <= 2.2204460492503131e-24f64

            return 1.;

        }

        if ux <= 0x7960000000000000u64 {

            // |x| <= f64::EPSILON

            // Power series of I0(x)Exp[-x] ~ 1 - x + O(x^2)

            return 1. - x;

        }

        if x.is_infinite() {

            return 0.;

        }

        return x + f64::NAN; // x = NaN

    }

    let xb = x.to_bits() & 0x7fff_ffff_ffff_ffff;

    if xb <= 0x4023000000000000u64 {

        // |x| <= 9.5

        if xb <= 0x400ccccccccccccdu64 {

            // |x| <= 3.6

            return i0e_0_to_3p6_exec(f64::from_bits(xb));

        } else if xb <= 0x401e000000000000u64 {

            // |x| <= 7.5

            return i0e3p6_to_7p5(f64::from_bits(xb));

        }

        return i0e_7p5_to_9p5(f64::from_bits(xb));

    }

    i0e_asympt(f64::from_bits(xb))

}

/**

Computes I0 on interval [-7.5; -3.6], [3.6; 7.5]

**/

#[inline]

fn i0e3p6_to_7p5(x: f64) -> f64 {

    let mut r = i0_3p6_to_7p5_dd(x);

    let v_exp = i0_exp(-x);

    r = DoubleDouble::quick_mult(r, v_exp);

    let err = f_fmla(

        r.hi,

        f64::from_bits(0x3c56a09e667f3bcd), // 2^-57.5

        f64::from_bits(0x3c00000000000000), // 2^-63

    );

    let ub = r.hi + (r.lo + err);

    let lb = r.hi + (r.lo - err);

    if ub == lb {

        return ub;

    }

    let v = eval_small_hard_3p6_to_7p5(x);

    let v_exp_accurate = i0_exp_accurate(-x);

    DoubleDouble::quick_mult(v, v_exp_accurate).to_f64()

}

#[inline]

fn i0e_0_to_3p6_exec(x: f64) -> f64 {

    let mut r = i0_0_to_3p6_dd(x);

    let v_exp = i0_exp(-x);

    r = DoubleDouble::quick_mult(r, v_exp);

    let err = f_fmla(

        r.hi,

        f64::from_bits(0x3c40000000000000), // 2^-59

        f64::from_bits(0x3be0000000000000), // 2^-66

    );

    let ub = r.hi + (r.lo + err);

    let lb = r.hi + (r.lo - err);

    if ub == lb {

        return ub;

    }

    let v = i0_0_to_3p6_hard(x);

    let v_exp_accurate = i0_exp_accurate(-x);

    DoubleDouble::quick_mult(v, v_exp_accurate).to_f64()

}

/**

Mid-interval [7.5;9.5] generated by Wolfram:

I0(x)=R(1/x)/sqrt(x)*Exp(x)

```text

<<FunctionApproximations`

ClearAll["Global`*"]

f[x_]:=Sqrt[x] Exp[-x] BesselI[0,x]

g[z_]:=f[1/z]

{err, approx}=MiniMaxApproximation[g[z],{z,{1/9.5,1/7.5},11,11},WorkingPrecision->120]

num=Numerator[approx][[1]];

den=Denominator[approx][[1]];

poly=den;

coeffs=CoefficientList[poly,z];

TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]

```

**/

#[inline]

fn i0e_7p5_to_9p5(x: f64) -> f64 {

    let dx = x;

    let recip = DoubleDouble::from_quick_recip(x);

    const P: [(u64, u64); 12] = [

        (0x3c778e3de1f76f48, 0x3fd988450531281b),

        (0xbcb572f6149f389e, 0xc01a786676fb4d3a),

        (0x3cf2f373365347ed, 0x405c0e8405fdb642),

        (0x3d276a94c8f1e627, 0xc0885e4718dfb761),

        (0x3d569f8a993434e2, 0x40b756d52d5fa90c),

        (0xbd6f953f7dd1a223, 0xc0c8818365c47790),

        (0xbd74247967fbf7b2, 0x40e8cf89daf87353),

        (0x3db449add7abb056, 0x41145d3c2d96d159),

        (0xbdc5cc822b71f891, 0xc123694c58fd039b),

        (0x3da2047ac1a6fba8, 0x415462e630bf3e7e),

        (0xbdc2f2c06eda6a95, 0xc14c6984ebdd6792),

        (0xbdf51fa85dafeca5, 0x4166a437c202d27b),

    ];

    let x2 = DoubleDouble::quick_mult(recip, recip);

    let x4 = DoubleDouble::quick_mult(x2, x2);

    let x8 = DoubleDouble::quick_mult(x4, x4);

    let e0 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(P[1]),

        DoubleDouble::from_bit_pair(P[0]),

    );

    let e1 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(P[3]),

        DoubleDouble::from_bit_pair(P[2]),

    );

    let e2 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(P[5]),

        DoubleDouble::from_bit_pair(P[4]),

    );

    let e3 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(P[7]),

        DoubleDouble::from_bit_pair(P[6]),

    );

    let e4 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(P[9]),

        DoubleDouble::from_bit_pair(P[8]),

    );

    let e5 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(P[11]),

        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),

        (0x3cde08e4cbf324d1, 0xc030b67bd69af0ca),

        (0x3cec5e4ee7e77024, 0x4071b54c0f58409c),

        (0xbd340e1131739e2f, 0xc09f140a738b14b3),

        (0x3d607673189d6644, 0x40cdb44bd822add2),

        (0xbd7793a4f1dd74d1, 0xc0e03fe2689b802d),

        (0xbd8415501228a87e, 0x410009beafea72cc),

        (0x3dcecdac2702661f, 0x4128c2073da9a447),

        (0xbdd8386404f3dec5, 0xc1389ec7d7186bf4),

        (0xbe06eb53a3e86436, 0x4168b7a2dc85ed0b),

        (0x3e098e2cefaf8299, 0xc1604f8cf34af02c),

        (0x3e1a5e496b547001, 0x41776b1e0153d1e9),

    ];

    let e0 = DoubleDouble::mul_add_f64(

        recip,

        DoubleDouble::from_bit_pair(Q[1]),

        f64::from_bits(0x3ff0000000000000),

    );

    let e1 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(Q[3]),

        DoubleDouble::from_bit_pair(Q[2]),

    );

    let e2 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(Q[5]),

        DoubleDouble::from_bit_pair(Q[4]),

    );

    let e3 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(Q[7]),

        DoubleDouble::from_bit_pair(Q[6]),

    );

    let e4 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(Q[9]),

        DoubleDouble::from_bit_pair(Q[8]),

    );

    let e5 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(Q[11]),

        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);

    let z = DoubleDouble::div(p_num, p_den);

    let r_sqrt = DoubleDouble::from_rsqrt_fast(dx);

    let r = z * r_sqrt;

    let err = f_fmla(

        r.hi,

        f64::from_bits(0x3bc0000000000000),

        f64::from_bits(0x392bdb8cdadbe111),

    );

    let up = r.hi + (r.lo + err);

    let lb = r.hi + (r.lo - err);

    if up != lb {

        return i0e_7p5_to_9p5_hard(x);

    }

    r.to_f64()

}

/**

Mid-interval [7.5;9.5] generated by Wolfram Mathematica:

I0(x)=R(1/x)/sqrt(x)*Exp(x)

Polynomial generated by Wolfram Mathematica:

```text

<<FunctionApproximations`

ClearAll["Global`*"]

f[x_]:=Sqrt[x] Exp[-x] BesselI[0,x]

g[z_]:=f[1/z]

{err,approx}=MiniMaxApproximation[g[z],{z,{1/9.5,1/7.5},13,13},WorkingPrecision->120]

poly=Numerator[approx][[1]];

coeffs=CoefficientList[poly,z];

TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]

poly=Denominator[approx][[1]];

coeffs=CoefficientList[poly,z];

TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]

```

**/

#[cold]

#[inline(never)]

fn i0e_7p5_to_9p5_hard(x: f64) -> f64 {

    static P: [DyadicFloat128; 14] = [

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -129,

            mantissa: 0xcc422a04_45cde144_75a3800b_45c38460_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -124,

            mantissa: 0xada66144_fcccc1a3_036f76b2_cabd6281_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -120,

            mantissa: 0xeabdda02_fa201d98_10e58d1f_7eb62bd7_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -116,

            mantissa: 0xbbfd3297_6f88a7df_5924587b_d5bdcdb8_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -113,

            mantissa: 0xfca29453_efe393bf_1651627b_7d543875_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -110,

            mantissa: 0xee7c7220_bbbd248e_bb6adac6_f9a5ce95_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -107,

            mantissa: 0xc07455dd_830ba705_414408c6_06732a5a_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -105,

            mantissa: 0xe2247793_b50cd0f0_80e8981d_933f75da_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -103,

            mantissa: 0xe14a9831_82582a0b_dd27e8b6_4ed9aac2_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -101,

            mantissa: 0xa3b2ae2f_5b64f37e_c1538435_34f02faf_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -100,

            mantissa: 0xbab73503_5b7e38d9_bbe4a84b_9007c6e8_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -99,

            mantissa: 0xa68911fc_5d87bbe7_0d4fe854_5c681ac5_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -99,

            mantissa: 0x9e997222_55ef4045_fa9f311d_57d082a2_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -99,

            mantissa: 0xbe93656a_b0a4f32d_3ebbfdeb_b1cbb839_u128,

        },

    ];

    static Q: [DyadicFloat128; 14] = [

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -127,

            mantissa: 0x80000000_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -123,

            mantissa: 0xdaa34a7e_861dddff_a0642080_cd83dd65_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -118,

            mantissa: 0x93f05740_f4758772_bb9992f9_91e72795_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -115,

            mantissa: 0xeddcb810_054c9aab_fa7e4214_d59d18b0_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -111,

            mantissa: 0xa0180fcd_831ff6c0_ac2b8f02_37f3cfd1_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -108,

            mantissa: 0x97d25106_3b66907e_90b4f786_26daa0bb_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -106,

            mantissa: 0xf595ce38_aac16c11_001b874a_99603b45_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -103,

            mantissa: 0x912b3715_4aca68f6_5821c2ed_43d77111_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -101,

            mantissa: 0x90f97141_b896e2b6_38b87354_8945a43c_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -100,

            mantissa: 0xd3e5f967_89097d6b_3a3060fe_852ff580_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -99,

            mantissa: 0xf0d6de35_939da009_9ced21fd_48af7281_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -98,

            mantissa: 0xd2a0b183_6ac613b2_6745ce1d_8ed1c323_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -98,

            mantissa: 0xbb9c089a_b7e939a2_732b3fb5_2e66cd77_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -98,

            mantissa: 0xcb2107c3_736bef81_609718c0_ba82cd8e_u128,

        },

    ];

    let recip = DyadicFloat128::accurate_reciprocal(x);

    let mut p_num = P[13];

    for i in (0..13).rev() {

        p_num = recip * p_num + P[i];

    }

    let mut p_den = Q[13];

    for i in (0..13).rev() {

        p_den = recip * p_den + Q[i];

    }

    let z = p_num * p_den.reciprocal();

    let r_sqrt = bessel_rsqrt_hard(x, recip);

    (z * r_sqrt).fast_as_f64()

}

/**

I0(x)=R(1/x)*Exp(x)/sqrt(x)

Generated in Wolfram:

```text

<<FunctionApproximations`

ClearAll["Global`*"]

f[x_]:=Sqrt[x] Exp[-x] BesselI[0,x]

g[z_]:=f[1/z]

{err,approx, err1}=MiniMaxApproximation[g[z],{z,{1/709.3,1/9.5},11,11},WorkingPrecision->120]

poly=Numerator[approx];

coeffs=CoefficientList[poly,z];

TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]

poly=Denominator[approx];

coeffs=CoefficientList[poly,z];

TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]]

```

**/

#[inline]

fn i0e_asympt(x: f64) -> f64 {

    let dx = x;

    let recip = DoubleDouble::from_quick_recip(x);

    const P: [(u64, u64); 12] = [

        (0xbc7ca19c5d824c54, 0x3fd9884533d43651),

        (0x3cca32eb734e010e, 0xc03b7ca35caf42eb),

        (0x3d03af8238d6f25e, 0x408b92cfcaa7070f),

        (0xbd7a8ff7fdebed70, 0xc0d0a3be432cce93),

        (0xbdababdb579bb076, 0x410a77dc51f1804d),

        (0x3dc5e4e3c972832a, 0xc13cb0be2f74839c),

        (0x3e01076f9b102e38, 0x41653b064cc61661),

        (0xbe2157e700d445f4, 0xc184e1b076927460),

        (0xbdfa4577156dde56, 0x41999e9653f9dc5f),

        (0xbe47aa238a739ffe, 0xc1a130f6ded40c00),

        (0xbe331b560b6fbf4a, 0x419317f11e674cae),

        (0xbe0765596077d1e3, 0xc16024404db59d3f),

    ];

    let x2 = DoubleDouble::quick_mult(recip, recip);

    let x4 = DoubleDouble::quick_mult(x2, x2);

    let x8 = DoubleDouble::quick_mult(x4, x4);

    let e0 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(P[1]),

        DoubleDouble::from_bit_pair(P[0]),

    );

    let e1 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(P[3]),

        DoubleDouble::from_bit_pair(P[2]),

    );

    let e2 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(P[5]),

        DoubleDouble::from_bit_pair(P[4]),

    );

    let e3 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(P[7]),

        DoubleDouble::from_bit_pair(P[6]),

    );

    let e4 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(P[9]),

        DoubleDouble::from_bit_pair(P[8]),

    );

    let e5 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(P[11]),

        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),

        (0xbcf687703e843d07, 0xc051418f1c4dd0b9),

        (0x3d468ab92cb87a0b, 0x40a15891d823e48d),

        (0x3d8bfc17e5183376, 0xc0e4fce40dd82796),

        (0xbdccbbcc2ecf8d4c, 0x4120beef869c61ec),

        (0xbdf42170b4d5d150, 0xc1523ad18834a7ed),

        (0xbe0eaa32f405afd4, 0x417b24ec57a8f480),

        (0x3e3ec900705e7252, 0xc19af2a91d23d62e),

        (0x3e3e220e274fa46b, 0x41b0cb905cc99ff5),

        (0xbe46c6c61dee11f6, 0xc1b7452e50518520),

        (0x3e3ed0fc063187bf, 0x41ac1772d1749896),

        (0xbe11c578f04f4be1, 0xc180feb5b2ca47cb),

    ];

    let e0 = DoubleDouble::mul_add_f64(

        recip,

        DoubleDouble::from_bit_pair(Q[1]),

        f64::from_bits(0x3ff0000000000000),

    );

    let e1 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(Q[3]),

        DoubleDouble::from_bit_pair(Q[2]),

    );

    let e2 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(Q[5]),

        DoubleDouble::from_bit_pair(Q[4]),

    );

    let e3 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(Q[7]),

        DoubleDouble::from_bit_pair(Q[6]),

    );

    let e4 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(Q[9]),

        DoubleDouble::from_bit_pair(Q[8]),

    );

    let e5 = DoubleDouble::mul_add(

        recip,

        DoubleDouble::from_bit_pair(Q[11]),

        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);

    let z = DoubleDouble::div(p_num, p_den);

    let r_sqrt = DoubleDouble::from_rsqrt_fast(dx);

    let r = z * r_sqrt;

    let err = f_fmla(

        r.hi,

        f64::from_bits(0x3c40000000000000), // 2^-59

        f64::from_bits(0x3be0000000000000), // 2^-65

    );

    let up = r.hi + (r.lo + err);

    let lb = r.hi + (r.lo - err);

    if up != lb {

        return i0e_asympt_hard(x);

    }

    lb

}

/**

I0(x)=R(1/x)*Exp(x)/sqrt(x)

Generated in Wolfram:

```text

<<FunctionApproximations`

ClearAll["Global`*"]

f[x_]:=Sqrt[x] Exp[-x] BesselI[0,x]

g[z_]:=f[1/z]

{err,approx, err1}=MiniMaxApproximation[g[z],{z,{1/709.3,1/9.5},15,15},WorkingPrecision->120]

poly=Numerator[approx];

coeffs=CoefficientList[poly,z];

TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]

poly=Denominator[approx];

coeffs=CoefficientList[poly,z];

TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]

```

**/

#[cold]

#[inline(never)]

fn i0e_asympt_hard(x: f64) -> f64 {

    static P: [DyadicFloat128; 16] = [

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -129,

            mantissa: 0xcc42299e_a1b28468_7e5aad4a_70b749c4_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -122,

            mantissa: 0xabb4209d_ca11bdaa_186bef7f_390e6b77_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -116,

            mantissa: 0x8a2725e2_4749db25_625ad1f2_12a2a16c_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -111,

            mantissa: 0x8b4c2cd4_9e5d0c8b_c9be4d3e_781bb676_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -107,

            mantissa: 0xc33fad7c_40599f7d_713e3081_6b5ad791_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -103,

            mantissa: 0xc81da271_b623ba88_0be032b5_827d92fa_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -99,

            mantissa: 0x99ec4975_b6aa7cae_7692a287_ed8ae47c_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -96,

            mantissa: 0xb3aa4745_fc2dd441_2dbd3e3c_d4539687_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -93,

            mantissa: 0x9f14edc2_6882afca_29d2a067_dc459729_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -91,

            mantissa: 0xd35c4d01_78d8cec6_fc8ae0ee_834da837_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -89,

            mantissa: 0xcdc529c7_6e082342_faad3073_07a9b61f_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -87,

            mantissa: 0x8ccac88f_2598c8a6_423b1f42_63591cb9_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -87,

            mantissa: 0xfc044cfb_a20f0885_93d58660_17819ed5_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -86,

            mantissa: 0x813a700c_74d23f52_f81b179d_7ff0da9f_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -88,

            mantissa: 0xe6c43da4_297216bf_bdd987cb_636906cf_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -91,

            mantissa: 0xa4998323_649c3cf2_64477869_3d2c6afd_u128,

        },

    ];

    static Q: [DyadicFloat128; 16] = [

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -127,

            mantissa: 0x80000000_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -121,

            mantissa: 0xd772d5fd_a7077638_6e007274_d83b013c_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -115,

            mantissa: 0xad914ef0_451ced2e_515657ef_fc7eeb53_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -110,

            mantissa: 0xaf41180c_dffe96e5_f192fa40_0b1bff1e_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -106,

            mantissa: 0xf60dc728_241f71fd_5b93e653_ccbedace_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -102,

            mantissa: 0xfcaefef9_39cf96e7_3cb75a98_da5e9761_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -98,

            mantissa: 0xc2d2c837_5789587a_13ef38c6_a24c3413_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -95,

            mantissa: 0xe41725c3_51d14486_a650044e_e8588f7b_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -92,

            mantissa: 0xcabeed9b_5e2e888d_81d32b11_d289a624_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -89,

            mantissa: 0x8764876d_11ad6607_8a8d5382_3ffe82d9_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -87,

            mantissa: 0x84c9f9e5_6a5f5034_ad2c8512_16cb1ba1_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -86,

            mantissa: 0xb7c1d143_a15d8aab_03a7fa3e_b7d07a36_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -85,

            mantissa: 0xa78f8257_4658040f_7a1ad39c_91ea9483_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -85,

            mantissa: 0xb231e0ca_b729a404_44c38f52_be208507_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -86,

            mantissa: 0xae317981_42349081_8bc68b28_f69b8e49_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -89,

            mantissa: 0xb451abd3_5cd79c6d_7e578164_32f16da1_u128,

        },

    ];

    let recip = DyadicFloat128::accurate_reciprocal(x);

    let mut p_num = P[15];

    for i in (0..15).rev() {

        p_num = recip * p_num + P[i];

    }

    let mut p_den = Q[15];

    for i in (0..15).rev() {

        p_den = recip * p_den + Q[i];

    }

    let z = p_num * p_den.reciprocal();

    let r_sqrt = bessel_rsqrt_hard(x, recip);

    (z * r_sqrt).fast_as_f64()

}

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn test_i0e() {

        assert_eq!(f_i0e(0.00000000000000000000000000052342), 1.0);

        assert_eq!(f_i0e(f64::EPSILON), 0.9999999999999998);

        assert_eq!(f_i0e(9.500000000005492,), 0.13125126081422883);

        assert!(f_i0e(f64::NAN).is_nan());

        assert_eq!(f_i0e(f64::INFINITY), 0.);

        assert_eq!(f_i0e(f64::NEG_INFINITY), 0.);

        assert_eq!(f_i0e(7.500000000788034), 0.14831583006929877);

        assert_eq!(f_i0e(715.), 0.014922205745802662);

        assert_eq!(f_i0e(12.), 0.11642622121344044);

        assert_eq!(f_i0e(16.), 0.10054412736125203);

        assert_eq!(f_i0e(18.432), 0.09357372647647);

        assert_eq!(f_i0e(26.432), 0.07797212360059241);

        assert_eq!(f_i0e(0.2), 0.8269385516343293);

        assert_eq!(f_i0e(7.5), 0.14831583007739552);

        assert_eq!(f_i0e(-1.5), 0.36743360905415834);

        assert_eq!(i0e_asympt_hard(18.432), 0.09357372647647);

    }

}