/*

 * // 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::bessel::i1::i1_0_to_7p75;

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 1

///

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

pub fn f_i1e(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-24

            return x * 0.5;

        }

        if ux <= 0x7960000000000000u64 {

            // |x| <= f64::EPSILON

            // Power series of I1(x)*exp(-|x|) ~ x/2 - x^2/2 + O(x^3)

            return f_fmla(x, -x * 0.5, x * 0.5);

        }

        if x.is_infinite() {

            return 0.;

        }

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

    }

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

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

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

    if xb < 0x401f000000000000u64 {

        // |x| <= 7.75

        let v_exp = i0_exp(-f64::from_bits(xb));

        let vi1 = i1_0_to_7p75(f64::from_bits(xb));

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

        return f64::copysign(r.to_f64(), sign_scale);

    }

    i1e_asympt(f64::from_bits(xb), sign_scale)

}

/**

Asymptotic expansion for I1.

Computes:

sqrt(x) * exp(-x) * I1(x) = Pn(1/x)/Qn(1/x)

hence:

I1(x)exp(-|x|) = Pn(1/x)/Qm(1/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 i1e_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 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(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);

    }

    i1e_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)exp(-|x|) = Pn(1/x)/Qm(1/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 i1e_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);

    (z * r_sqrt).fast_as_f64() * sign_scale

}

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn test_fi1e() {

        assert_eq!(f_i1e(f64::EPSILON), 1.1102230246251563e-16);

        assert_eq!(f_i1e(7.750000000757874), 0.13605110007443239);

        assert_eq!(f_i1e(7.482812501363189), 0.13818116726273896);

        assert_eq!(f_i1e(-7.750000000757874), -0.13605110007443239);

        assert_eq!(f_i1e(-7.482812501363189), -0.13818116726273896);

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

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

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

        assert_eq!(f_i1e(0.01), 0.004950311047118276);

        assert_eq!(f_i1e(-0.01), -0.004950311047118276);

        assert_eq!(f_i1e(-9.01), -0.12716101566063667);

        assert_eq!(f_i1e(9.01), 0.12716101566063667);

        assert_eq!(f_i1e(763.), 0.014435579051182581);

        assert_eq!(i1e_asympt_hard(9.01, 1.), 0.12716101566063667);

    }

}