pxfm 0.1.28

/*
 * // 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::i0::bessel_rsqrt_hard;
use crate::bessel::i0_exp;
use crate::common::f_fmla;
use crate::double_double::DoubleDouble;
use crate::dyadic_float::{DyadicFloat128, DyadicSign};
use crate::exponents::rational128_exp;
use crate::polyeval::{f_estrin_polyeval5, f_polyeval6};

/// Modified Bessel of the first kind of order 1
///
/// Max found ULP 0.5
pub fn f_i1(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 {
            // Power series of I1(x) ~ x/2 + O(x^3)
            // |x| <= 2.2204460492503131e-24
            return x * 0.5;
        }
        if ux <= 0x7960000000000000u64 {
            // |x| <= f64::EPSILON
            // Power series of I1(x) ~ x/2 + x^3/16 + O(x^4)
            const A0: f64 = 1. / 2.;
            const A1: f64 = 1. / 16.;
            let r0 = f_fmla(x, x * A1, A0);
            return r0 * x;
        }
        if x.is_infinite() {
            return if x.is_sign_positive() {
                f64::INFINITY
            } else {
                f64::NEG_INFINITY
            };
        }
        return x + f64::NAN; // x == NaN
    }

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

    if xb >= 0x40864fe69ff9fec8u64 {
        // |x| >= 713.9876098185423
        return if x.is_sign_negative() {
            f64::NEG_INFINITY
        } else {
            f64::INFINITY
        };
    }

    static SIGN: [f64; 2] = [1., -1.];

    let sign_scale = SIGN[x.is_sign_negative() as usize];

    if xb < 0x401f000000000000u64 {
        // |x| <= 7.75
        return f64::copysign(i1_0_to_7p75(f64::from_bits(xb)).to_f64(), sign_scale);
    }

    i1_asympt(f64::from_bits(xb), sign_scale)
}

/**
Computes
I1(x) = x/2 * (1 + 1 * (x/2)^2 + (x/2)^4 * P((x/2)^2))

Polynomial generated by Wolfram Mathematica:
```text
<<FunctionApproximations`
ClearAll["Global`*"]
f[x_]:=(BesselI[1,x]*2/x-1-1/2(x/2)^2)/(x/2)^4
g[z_]:=f[2 Sqrt[z]]
{err,approx}=MiniMaxApproximation[g[z],{z,{0.000000001,7.5},9,9},WorkingPrecision->60]
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}]]
```
**/
#[inline]
pub(crate) fn i1_0_to_7p75(x: f64) -> DoubleDouble {
    let half_x = x * 0.5; // this is exact
    let eval_x = DoubleDouble::from_exact_mult(half_x, half_x);

    const P: [(u64, u64); 5] = [
        (0x3c55555555555555, 0x3fb5555555555555),
        (0x3c1253e1df138479, 0x3f7304597c4fbd4c),
        (0x3bcec398b7059ee9, 0x3f287b5b01f6b9c0),
        (0xbb7354e7c92c4f77, 0x3ed21de117470d10),
        (0xbb1d35ac0d7923cc, 0x3e717f3714dddc04),
    ];

    let ps_num = f_estrin_polyeval5(
        eval_x.hi,
        f64::from_bits(0x3e063684ca1944a4),
        f64::from_bits(0x3d92ac4a0e48a9bb),
        f64::from_bits(0x3d1425988b0b0aea),
        f64::from_bits(0x3c899839e74ddefc),
        f64::from_bits(0x3bed8747bcdd1e61),
    );

    let mut p_num = DoubleDouble::mul_f64_add(eval_x, ps_num, DoubleDouble::from_bit_pair(P[4]));
    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[3]));
    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[2]));
    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[1]));
    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[0]));

    const Q: [(u64, u64); 4] = [
        (0x0000000000000000, 0x3ff0000000000000),
        (0xbc3e59afb81ac7ea, 0xbf9c4848e0661d70),
        (0x3bd62fa3dbc1a19c, 0x3f38a9eafcd7e674),
        (0x3b6f4688b9eab7d0, 0xbecbfdec51454533),
    ];

    let ps_den = f_polyeval6(
        eval_x.hi,
        f64::from_bits(0x3e56e7cde9266a32),
        f64::from_bits(0xbddc316dff4a672f),
        f64::from_bits(0x3d5a43aaee30ebb5),
        f64::from_bits(0xbcd1fb023f4f1fa0),
        f64::from_bits(0x3c4089ede324209f),
        f64::from_bits(0xbb9f64f47ba69604),
    );

    let mut p_den = DoubleDouble::mul_f64_add(eval_x, ps_den, DoubleDouble::from_bit_pair(Q[3]));
    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[2]));
    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[1]));
    p_den = DoubleDouble::mul_add_f64(eval_x, p_den, f64::from_bits(0x3ff0000000000000));

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

    let eval_sqr = DoubleDouble::quick_mult(eval_x, eval_x);

    let mut z = DoubleDouble::mul_f64_add_f64(eval_x, 0.5, 1.);
    z = DoubleDouble::mul_add(p, eval_sqr, z);

    let x_over_05 = DoubleDouble::from_exact_mult(x, 0.5);

    let r = DoubleDouble::quick_mult(z, x_over_05);

    let err = f_fmla(
        r.hi,
        f64::from_bits(0x3c40000000000000), // 2^-59
        f64::from_bits(0x3be0000000000000), // 2^-65
    );
    let ub = r.hi + (r.lo + err);
    let lb = r.hi + (r.lo - err);
    if ub == lb {
        return r;
    }
    i1_0_to_7p5_hard(x)
}

// Polynomial generated by Wolfram Mathematica:
// I1(x) = x/2 * (1 + 1 * (x/2)^2 + (x/2)^4 * P((x/2)^2))
// <<FunctionApproximations`
// ClearAll["Global`*"]
// f[x_]:=(BesselI[1,x]*2/x-1-1/2(x/2)^2)/(x/2)^4
// g[z_]:=f[2 Sqrt[z]]
// {err,approx}=MiniMaxApproximation[g[z],{z,{0.000000001,7.5},9,9},WorkingPrecision->60]
// 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)]
pub(crate) fn i1_0_to_7p5_hard(x: f64) -> DoubleDouble {
    const ONE_OVER_4: f64 = 1. / 4.;
    let eval_x = DoubleDouble::quick_mult_f64(DoubleDouble::from_exact_mult(x, x), ONE_OVER_4);

    const P: [(u64, u64); 10] = [
        (0x3c55555555555555, 0x3fb5555555555555),
        (0x3c1253e1df138479, 0x3f7304597c4fbd4c),
        (0x3bcec398b7059ee9, 0x3f287b5b01f6b9c0),
        (0xbb7354e7c92c4f77, 0x3ed21de117470d10),
        (0xbb1d35ac0d7923cc, 0x3e717f3714dddc04),
        (0xba928dee24678e32, 0x3e063684ca1944a4),
        (0xba36aa59912fcbed, 0x3d92ac4a0e48a9bb),
        (0x39bad76f18b5ad37, 0x3d1425988b0b0aea),
        (0xb923a6bab6928df4, 0x3c899839e74ddefc),
        (0x3864356cdfa7b321, 0x3bed8747bcdd1e61),
    ];

    let mut p_num = DoubleDouble::mul_add(
        eval_x,
        DoubleDouble::from_bit_pair(P[9]),
        DoubleDouble::from_bit_pair(P[8]),
    );
    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[7]));
    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[6]));
    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[5]));
    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[4]));
    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[3]));
    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[2]));
    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[1]));
    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[0]));

    const Q: [(u64, u64); 10] = [
        (0x0000000000000000, 0x3ff0000000000000),
        (0xbc3e59afb81ac7ea, 0xbf9c4848e0661d70),
        (0x3bd62fa3dbc1a19c, 0x3f38a9eafcd7e674),
        (0x3b6f4688b9eab7d0, 0xbecbfdec51454533),
        (0x3af0fb4a17103ef8, 0x3e56e7cde9266a32),
        (0xba71755779c6d4bd, 0xbddc316dff4a672f),
        (0x39cf8ed8d449e2c6, 0x3d5a43aaee30ebb5),
        (0x39704e900a373874, 0xbcd1fb023f4f1fa0),
        (0xb8e33e87e4bffbb1, 0x3c4089ede324209f),
        (0x380fb09b3fd49d5c, 0xbb9f64f47ba69604),
    ];

    let mut p_den = DoubleDouble::mul_add(
        eval_x,
        DoubleDouble::from_bit_pair(Q[9]),
        DoubleDouble::from_bit_pair(Q[8]),
    );
    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[7]));
    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[6]));
    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[5]));
    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[4]));
    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[3]));
    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[2]));
    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[1]));
    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[0]));

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

    let eval_sqr = DoubleDouble::quick_mult(eval_x, eval_x);

    let mut z = DoubleDouble::mul_f64_add_f64(eval_x, 0.5, 1.);
    z = DoubleDouble::mul_add(p, eval_sqr, z);

    let x_over_05 = DoubleDouble::from_exact_mult(x, 0.5);

    DoubleDouble::quick_mult(z, x_over_05)
}

/**
Asymptotic expansion for I1.

Computes:
sqrt(x) * exp(-x) * I1(x) = Pn(1/x)/Qn(1/x)
hence:
I1(x) = Pn(1/x)/Qm(1/x)*exp(x)/sqrt(x)

Generated by Wolfram Mathematica:
```text
<<FunctionApproximations`
ClearAll["Global`*"]
f[x_]:=Sqrt[x] Exp[-x] BesselI[1,x]
g[z_]:=f[1/z]
{err,approx}=MiniMaxApproximation[g[z],{z,{1/713.98,1/7.75},11,11},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}]]
```
**/
#[inline]
fn i1_asympt(x: f64, sign_scale: f64) -> f64 {
    let dx = x;
    let recip = DoubleDouble::from_quick_recip(x);
    const P: [(u64, u64); 12] = [
        (0xbc73a823f28a2f5e, 0x3fd9884533d43651),
        (0x3cc0d5bb78e674b3, 0xc0354325c8029263),
        (0x3d20e1155aaaa283, 0x4080c09b027c46a4),
        (0xbd5b90dcf81b99c1, 0xc0bfc1311090c839),
        (0xbd98f2fda9e8fa1b, 0x40f3bb81bb190ae2),
        (0xbdcec960752b60da, 0xc1207c0bbbc31cd9),
        (0x3dd3c9a299c9c41f, 0x414253e25c4584af),
        (0xbde82e7b9a3e1acc, 0xc159a656aece377c),
        (0x3e0d3d30d701a8ab, 0x416398df24c74ef2),
        (0xbdf57b85ab7006e2, 0xc151fd119be1702b),
        (0x3dd760928f4515fd, 0xc1508327e42639bc),
        (0x3dc09e71bc648589, 0x4143e4933afa621c),
    ];

    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),
        (0xbcb334d5a476d9ad, 0xc04a75f94c1a0c1a),
        (0xbd324d58ed98bfae, 0x4094b00e60301c42),
        (0x3d7c8725666c4360, 0xc0d36b9d28d45928),
        (0x3d7f8457c2945822, 0x4107d6c398a174ed),
        (0x3dbc655ea216594b, 0xc1339393e6776e38),
        (0xbdebb5dffef78272, 0x415537198d23f6a1),
        (0xbdb577f8abad883e, 0xc16c6c399dcd6949),
        (0x3e14261c5362f109, 0x4173c02446576949),
        (0x3dc382ededad42c5, 0xc1547dff5770f4ec),
        (0xbe05c0f74d4c7956, 0xc165c88046952562),
        (0xbdbf9145927aa2c7, 0x414395e46bc45d5b),
    ];

    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 e = i0_exp(dx * 0.5);
    let r_sqrt = DoubleDouble::from_rsqrt_fast(dx);

    let r = DoubleDouble::quick_mult(z * r_sqrt * e, e);

    let err = f_fmla(
        r.hi,
        f64::from_bits(0x3c40000000000000), // 2^-59
        f64::from_bits(0x3ba0000000000000), // 2^-69
    );
    let up = r.hi + (r.lo + err);
    let lb = r.hi + (r.lo - err);
    if up == lb {
        return f64::copysign(r.to_f64(), sign_scale);
    }
    i1_asympt_hard(x, sign_scale)
}

/**
Asymptotic expansion for I1.

Computes:
sqrt(x) * exp(-x) * I1(x) = Pn(1/x)/Qn(1/x)
hence:
I1(x) = Pn(1/x)/Qm(1/x)*exp(x)/sqrt(x)

Generated by Wolfram Mathematica:
```text
<<FunctionApproximations`
ClearAll["Global`*"]
f[x_]:=Sqrt[x] Exp[-x] BesselI[1,x]
g[z_]:=f[1/z]
{err,approx}=MiniMaxApproximation[g[z],{z,{1/713.98,1/7.75},15,15},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 i1_asympt_hard(x: f64, sign_scale: f64) -> f64 {
    static P: [DyadicFloat128; 16] = [
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -129,
            mantissa: 0xcc42299e_a1b28468_bea7da47_28f13acc_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -124,
            mantissa: 0xda979406_3df6e66f_cf31c3f5_f194b48c_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -120,
            mantissa: 0xd60b7b96_c958929b_cabe1d8c_3d874767_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -113,
            mantissa: 0xd27aad9a_8fb38d56_46ab4510_8479306e_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -108,
            mantissa: 0xe0167305_c451bd1f_d2f17b68_5c62e2ff_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -103,
            mantissa: 0x8f6d238f_c80d8e4a_08c130f6_24e1c925_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -100,
            mantissa: 0xfe32280f_2ea99024_d9924472_92d7ac8f_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -96,
            mantissa: 0xa48815ac_d265609f_da4ace94_811390b2_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -93,
            mantissa: 0x9ededfe5_833b4cc1_731efd5c_f8729c6c_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -91,
            mantissa: 0xe5b43203_2784ae6a_f7458556_0a8308ea_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -89,
            mantissa: 0xf5df279a_3fb4ef60_8d10adee_7dd2f47b_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -87,
            mantissa: 0xbdb59963_7a757ed1_87280e0e_7f93ca2b_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -86,
            mantissa: 0xc87fdea5_53177ca8_c91de5fb_3f8f78d3_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -85,
            mantissa: 0x846d16a7_4663ef5d_ad42d599_5bc726b8_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -86,
            mantissa: 0xb3ed2055_74262d95_389f33e4_2ac3774a_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -88,
            mantissa: 0xa511aa32_c18c34e4_3d029a90_a71b7a55_u128,
        },
    ];
    static Q: [DyadicFloat128; 16] = [
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -127,
            mantissa: 0x80000000_00000000_00000000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -122,
            mantissa: 0x877b771a_ad8f5fd3_5aacf5f9_f04ee9de_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -118,
            mantissa: 0x89475ecd_9c84361e_800c8a3a_c8af23bf_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -111,
            mantissa: 0x837d1196_cf2723f1_23b54da8_225efe05_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -106,
            mantissa: 0x8ae3aecb_15355751_a9ee12e5_a4dd9dde_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -102,
            mantissa: 0xb0886afa_bc13f996_ab45d252_75c8f587_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -98,
            mantissa: 0x9b37d7cd_b114b86b_7d14a389_26599aa1_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -95,
            mantissa: 0xc716bf54_09d5dd9f_bc16679b_93aaeca4_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -92,
            mantissa: 0xbe0cd82e_c8af8371_ab028ed9_c7902dd2_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -89,
            mantissa: 0x875f8d91_8ef5d434_a39d00f9_2aed3d2a_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -87,
            mantissa: 0x8e030781_5aa4ce7f_70156b82_8b216b7c_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -86,
            mantissa: 0xd4dd2687_92646fbd_5ea2d422_da64fc0b_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -85,
            mantissa: 0xd6d72ab3_64b4a827_0499af0f_13a51a80_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -84,
            mantissa: 0x828f4e8b_728747a9_2cebe54a_810e2681_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -85,
            mantissa: 0x91570096_36a3fcfb_6b936d44_68dda1be_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -89,
            mantissa: 0xf082ad00_86024ed4_dd31613b_ec41e3f8_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);
    let f_exp = rational128_exp(x);
    (z * r_sqrt * f_exp).fast_as_f64() * sign_scale
}

#[cfg(test)]
mod tests {
    use super::*;
    #[test]
    fn test_fi1() {
        assert_eq!(
            f_i1(0.0000000000000000000000000000000006423424234121),
            3.2117121170605e-34
        );
        assert_eq!(f_i1(7.750000000757874), 315.8524811496668);
        assert_eq!(f_i1(7.482812501363189), 245.58002285881892);
        assert_eq!(f_i1(-7.750000000757874), -315.8524811496668);
        assert_eq!(f_i1(-7.482812501363189), -245.58002285881892);
        assert!(f_i1(f64::NAN).is_nan());
        assert_eq!(f_i1(f64::INFINITY), f64::INFINITY);
        assert_eq!(f_i1(f64::NEG_INFINITY), f64::NEG_INFINITY);
        assert_eq!(f_i1(0.01), 0.005000062500260418);
        assert_eq!(f_i1(-0.01), -0.005000062500260418);
        assert_eq!(f_i1(-9.01), -1040.752038018038);
        assert_eq!(f_i1(9.01), 1040.752038018038);
    }
}