/*

 * // Copyright (c) Radzivon Bartoshyk 8/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;

/// Modified bessel of the first kind of order 2

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

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

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

        // |x| == 0, |x| == inf, x == NaN

        if ux == 0 {

            // |x| == 0

            return 0.;

        }

        if x.is_infinite() {

            return f64::INFINITY;

        }

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

    }

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

    if xb < 0x401f000000000000u64 {

        // |x| < 7.75

        if xb <= 0x3cb0000000000000u64 {

            // |x| <= f64::EPSILON

            // Power series of I2(x) ~ x^2/8 + O(x^4)

            const R: f64 = 1. / 8.;

            let x2 = x * x * R;

            return x2;

        }

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

    }

    if xb >= 0x40864feaeefb23b8 {

        // x >= 713.9897136326099

        return f64::INFINITY;

    }

    i2_asympt(f64::from_bits(xb))

}

/**

Computes

I2(x) = x^2 * R(x^2)

Generated by Wolfram Mathematica:

```text

<<FunctionApproximations`

ClearAll["Global`*"]

f[x_]:=BesselI[2,x]/x^2

g[z_]:=f[Sqrt[z]]

{err,approx}=MiniMaxApproximation[g[z],{z,{0.000000000001,7.75},11,11},WorkingPrecision->75]

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 i2_small(x: f64) -> f64 {

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

        (0x0000000000000000, 0x3fc0000000000000),

        (0x3c247833fda9de9a, 0x3f8387c6e72a1b5f),

        (0xbbccaf0be91261a6, 0x3f30ba88efff56fa),

        (0x3b57c911bfebe1d7, 0x3ecc62e53d061300),

        (0x3af3b963f26a3d05, 0x3e5bb090327a14e1),

        (0x3a898bff9d40e030, 0x3de0d29c3d37e5b5),

        (0xb9f2f63c80d377db, 0x3d5a9e365f1bf6e0),

        (0xb965e6d78e1c2b65, 0x3ccbf7ef0929b813),

        (0xb8da83d7d40e7310, 0x3c33737520046f4d),

        (0xb83f811d5aa3f36e, 0x3b91506558dab318),

        (0xb78e601bf5c998c3, 0x3ae2013b3e858bd1),

        (0xb6c8185c51734ed8, 0x3a20cc277a5051ba),

    ];

    let x_sqr = DoubleDouble::from_exact_mult(x, x);

    let x2 = x_sqr * x_sqr;

    let x4 = x2 * x2;

    let x8 = x4 * x4;

    let e0 = DoubleDouble::mul_add_f64(

        x_sqr,

        DoubleDouble::from_bit_pair(P[1]),

        f64::from_bits(0x3fc0000000000000),

    );

    let e1 = DoubleDouble::mul_add(

        x_sqr,

        DoubleDouble::from_bit_pair(P[3]),

        DoubleDouble::from_bit_pair(P[2]),

    );

    let e2 = DoubleDouble::mul_add(

        x_sqr,

        DoubleDouble::from_bit_pair(P[5]),

        DoubleDouble::from_bit_pair(P[4]),

    );

    let e3 = DoubleDouble::mul_add(

        x_sqr,

        DoubleDouble::from_bit_pair(P[7]),

        DoubleDouble::from_bit_pair(P[6]),

    );

    let e4 = DoubleDouble::mul_add(

        x_sqr,

        DoubleDouble::from_bit_pair(P[9]),

        DoubleDouble::from_bit_pair(P[8]),

    );

    let e5 = DoubleDouble::mul_add(

        x_sqr,

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

        (0xbc0ba42af56ed76b, 0xbf7cd8e6e2b39f60),

        (0x3b90697aa005e598, 0x3efa0260394e1a3d),

        (0xbb0c7ccde1f63c82, 0xbe6f1766ec64e492),

        (0x3a63f18409bc336f, 0x3ddb80b6b5abad98),

        (0xb9e0cd49f22132fe, 0xbd42ff9b55d553da),

        (0xb934bfe64905d309, 0x3ca50814fa258ebc),

        (0x38a1e35c2d6860f4, 0xbc02c4f2faca2195),

        (0xb7ff39e268277e4e, 0x3b5aa545a2c1f16d),

        (0xb71053f58545760c, 0xbaacde4c133d42d1),

        (0xb68d0c2ccab0ae5b, 0x39f5a965b92b06bc),

        (0xb5dc35bda16bee7b, 0xb931375b1c9cfbc7),

    ];

    let e0 = DoubleDouble::mul_add_f64(

        x_sqr,

        DoubleDouble::from_bit_pair(Q[1]),

        f64::from_bits(0x3ff0000000000000),

    );

    let e1 = DoubleDouble::mul_add(

        x_sqr,

        DoubleDouble::from_bit_pair(Q[3]),

        DoubleDouble::from_bit_pair(Q[2]),

    );

    let e2 = DoubleDouble::mul_add(

        x_sqr,

        DoubleDouble::from_bit_pair(Q[5]),

        DoubleDouble::from_bit_pair(Q[4]),

    );

    let e3 = DoubleDouble::mul_add(

        x_sqr,

        DoubleDouble::from_bit_pair(Q[7]),

        DoubleDouble::from_bit_pair(Q[6]),

    );

    let e4 = DoubleDouble::mul_add(

        x_sqr,

        DoubleDouble::from_bit_pair(Q[9]),

        DoubleDouble::from_bit_pair(Q[8]),

    );

    let e5 = DoubleDouble::mul_add(

        x_sqr,

        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 p = DoubleDouble::div(p_num, p_den);

    DoubleDouble::quick_mult(p, x_sqr).to_f64()

}

/**

Asymptotic expansion for I2.

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

Generated in Wolfram:

```text

<<FunctionApproximations`

ClearAll["Global`*"]

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

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

{err,approx}=MiniMaxApproximation[g[z],{z,{1/714.0,1/7.5},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 i2_asympt(x: f64) -> f64 {

    let dx = x;

    let recip = DoubleDouble::from_quick_recip(x);

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

        (0x3c718bb28ebc5f4e, 0x3fd9884533d43650),

        (0x3c96e15a87b6e1d1, 0xc0350acc9e5cb0f9),

        (0xbd20b212a79e08f5, 0x40809251af67598a),

        (0xbd563b7397df3a54, 0xc0bfc09ede682c8b),

        (0xbd5eb872cb057d91, 0x40f44253a9e1e1ab),

        (0x3d7614735e566fc5, 0xc121cbcd96dc8765),

        (0xbddc4f8df2010026, 0x4145a592e8ec74ad),

        (0x3dea227617b678a7, 0xc161df96fb6a9df9),

        (0x3e17c9690d906194, 0x41732c71397757f8),

        (0x3e0638226ce0b938, 0xc178893fde0e6ed7),

        (0xbe09d8dc4e7930ce, 0x417066abe24b31df),

        (0xbde152007ee29e54, 0xc150531da3f31b16),

    ];

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

        (0xbcd0d33e9e73b503, 0xc0496f5a09751d50),

        (0x3d2f9c44a069dc4b, 0x40934427187ac370),

        (0xbd69e2e5a3618381, 0xc0d19983f74fdf52),

        (0x3d88c69a62ae8b44, 0x410524fcaa71e85a),

        (0xbdc0345b806dd0bf, 0xc13120daf531b66b),

        (0xbdd35875712fff6f, 0x4152943a4f9f1c7f),

        (0xbdf8dd50e92553fd, 0xc169b83aeede08ea),

        (0x3e0800ecaa77f79e, 0x41746c61554a08ce),

        (0x3dd74fbc32c5f696, 0xc16ba2febd1932a3),

        (0x3dc23eb2c943b539, 0x413574ae68b6b378),

        (0xbd95d86c5c94cd65, 0xc104adac99eaa90c),

    ];

    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 mut e = i0_exp(dx * 0.5);

    e = DoubleDouble::from_exact_add(e.hi, e.lo);

    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 r.to_f64();

    }

    i2_asympt_hard(x)

}

/**

Asymptotic expansion for I2.

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

Generated in Wolfram:

```text

<<FunctionApproximations`

ClearAll["Global`*"]

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

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

{err,approx}=MiniMaxApproximation[g[z],{z,{1/714.0,1/7.5},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 i2_asympt_hard(x: f64) -> f64 {

    static P: [DyadicFloat128; 16] = [

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -129,

            mantissa: 0xcc42299e_a1b28468_3bb16645_ba1dc793_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -123,

            mantissa: 0xe202abf7_de10e93f_2a2e6a0f_af69c788_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -118,

            mantissa: 0xf70296c3_ad33bde6_866cfd01_0e846cfc_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -113,

            mantissa: 0xa83df971_736c4e6c_1a35479b_ad6d9172_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -109,

            mantissa: 0x9baa2015_9c5ca461_0aff0b62_54a70fdb_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -106,

            mantissa: 0xc70af95d_f95d14ad_1094ea1b_e46b2d2f_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -103,

            mantissa: 0xa838fb48_e79fb706_642da604_6a73b4f8_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -101,

            mantissa: 0x8fe29f37_02b1e876_39e88664_1c8b3b5d_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -100,

            mantissa: 0xc8e9a474_0a03f93a_16d2e7a9_627eba4e_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -95,

            mantissa: 0x8807d1f6_6d646a08_8c7e8900_12d6a5ed_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -93,

            mantissa: 0xe5c25026_97518024_36878256_fd81c08f_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -91,

            mantissa: 0xeaa075f0_f5151bed_95ec612f_ab9834a7_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -89,

            mantissa: 0x9b267222_82d5c666_348d7d1d_0fedfba4_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -88,

            mantissa: 0x81b45c4c_3e828396_1d5bdede_869c3b84_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -89,

            mantissa: 0xf4495d43_4bc8dba6_42bdb5d6_c8ba2c9c_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -90,

            mantissa: 0xc9b29546_0c226270_bb428035_587b6d6a_u128,

        },

    ];

    static Q: [DyadicFloat128; 16] = [

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -127,

            mantissa: 0x80000000_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -121,

            mantissa: 0x89e18bae_ca9629a1_26927ba2_fbdd66ab_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -116,

            mantissa: 0x92a90fc2_e905f634_4946e8a0_dd8e3874_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -112,

            mantissa: 0xc1742696_d29e3846_3e183737_29db8b68_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -108,

            mantissa: 0xabf61cc0_236a0e90_2572113d_fa339591_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -105,

            mantissa: 0xcff0fe90_dac1b08e_9a5740ae_b2984fc1_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -102,

            mantissa: 0x9ff36729_e407c538_cfcea3a7_63f39043_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -101,

            mantissa: 0xc86ff6a3_9b803a31_d385e9ea_83f9d751_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -98,

            mantissa: 0xb4a125b1_6cab70f3_0f314558_708843df_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -94,

            mantissa: 0x9670fd33_f83bcaa7_85cf2d82_c0bf8cd5_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -92,

            mantissa: 0xd70b4ea5_32fedb9d_78a3c047_05e650f4_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -90,

            mantissa: 0xb9c7904c_3f97b633_c2c0ad9b_ad573ede_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -89,

            mantissa: 0xc2023c21_5155e9fe_6fb17bb2_c865becd_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -89,

            mantissa: 0xd9400a5e_27c58803_22948cf3_6154ac49_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -90,

            mantissa: 0x87aa272d_6a9700b4_449a9db8_1a93b0ee_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -93,

            mantissa: 0xd1a86655_5b259611_dfc7affc_6ffb0e20_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()

}

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn test_i2() {

        assert_eq!(f_i2(7.750000000757874), 257.0034362785801);

        assert_eq!(f_i2(7.482812501363189), 198.26765887136534);

        assert_eq!(f_i2(-7.750000000757874), 257.0034362785801);

        assert_eq!(f_i2(-7.482812501363189), 198.26765887136534);

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

        assert_eq!(f_i2(f64::INFINITY), f64::INFINITY);

        assert_eq!(f_i2(f64::NEG_INFINITY), f64::INFINITY);

        assert_eq!(f_i2(0.01), 1.2500104166992188e-5);

        assert_eq!(f_i2(-0.01), 1.2500104166992188e-5);

        assert_eq!(f_i2(-9.01), 872.9250699638584);

        assert_eq!(f_i2(9.01), 872.9250699638584);

    }

}