pxfm 0.1.29

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