/*

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

    }

}