/*

 * // 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::j0f::j1f_rsqrt;

use crate::exponents::core_expf;

use crate::polyeval::{f_estrin_polyeval8, f_estrin_polyeval9};

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

///

/// ULP 0.5

pub fn f_i2f(x: f32) -> f32 {

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

    if ux >= 0xffu32 << 24 || ux == 0 {

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

        if ux == 0 {

            // |x| == 0

            return 0.;

        }

        if x.is_infinite() {

            return f32::INFINITY;

        }

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

    }

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

    if xb >= 0x42b7d875u32 {

        // |x| >= 91.92277 it's infinity

        return f32::INFINITY;

    }

    if xb <= 0x40f80000u32 {

        // |x| <= 7.75

        if xb <= 0x34000000u32 {

            // |x| <= f32::EPSILON

            let dx = x as f64;

            const R: f64 = 1. / 8.;

            return (dx * dx * R) as f32;

        }

        return i2f_small(f32::from_bits(xb));

    }

    i2f_asympt(f32::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},8,7},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 i2f_small(x: f32) -> f32 {

    let dx = x as f64;

    let x_sqr = dx * dx;

    let p_num = f_estrin_polyeval9(

        x_sqr,

        f64::from_bits(0x3fc0000000000000),

        f64::from_bits(0x3f831469a38d72c7),

        f64::from_bits(0x3f2f453dd3dd98f4),

        f64::from_bits(0x3ec8af52ee8fce9b),

        f64::from_bits(0x3e5589f2cb4e0ec9),

        f64::from_bits(0x3dd60fa268a4206d),

        f64::from_bits(0x3d4ab3091ee18d6b),

        f64::from_bits(0x3cb1efec43b15186),

        f64::from_bits(0x3c050992c6e9e63f),

    );

    let p_den = f_estrin_polyeval8(

        x_sqr,

        f64::from_bits(0x3ff0000000000000),

        f64::from_bits(0xbf82075d8e3f1476),

        f64::from_bits(0x3f03ef86564a284b),

        f64::from_bits(0xbe7c498fab4a57d8),

        f64::from_bits(0x3dec162ca0f68486),

        f64::from_bits(0xbd53bb6398461540),

        f64::from_bits(0x3cb265215261e64a),

        f64::from_bits(0xbc01cf52cc350e81),

    );

    let p = p_num / p_den;

    (p * x_sqr) as f32

}

/**

Asymptotic expansion for I2.

Computes:

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

hence:

I2(x) = Pn(1/x)/Qm(1/x)*exp(x)/sqrt(x)

Generated by Mathematica:

```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/92.3,1/7.5},8,8},WorkingPrecision->70]

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 i2f_asympt(x: f32) -> f32 {

    let dx = x as f64;

    let recip = 1. / dx;

    let p_num = f_estrin_polyeval9(

        recip,

        f64::from_bits(0x3fd9884533d45f46),

        f64::from_bits(0xc02b979526807e1e),

        f64::from_bits(0x406b1dd3e795bbed),

        f64::from_bits(0xc09e43629031ec91),

        f64::from_bits(0x40c48c03a39aec1d),

        f64::from_bits(0xc0e0f022ccb8807a),

        f64::from_bits(0x40f0302eeb22a776),

        f64::from_bits(0xc0f02b01549d38b8),

        f64::from_bits(0x40dad4e70f2bc264),

    );

    let p_den = f_estrin_polyeval9(

        recip,

        f64::from_bits(0x3ff0000000000000),

        f64::from_bits(0xc0405a71a88b191c),

        f64::from_bits(0x407e19f7d247d098),

        f64::from_bits(0xc0aeaac6e0ca17fe),

        f64::from_bits(0x40d2301702f40a98),

        f64::from_bits(0xc0e7e6c6c01841b3),

        f64::from_bits(0x40ed67317e9e46cc),

        f64::from_bits(0xc0d13786aa1ef416),

        f64::from_bits(0xc0a6c9cfe579ae22),

    );

    let z = p_num / p_den;

    let e = core_expf(x);

    let r_sqrt = j1f_rsqrt(dx);

    (z * r_sqrt * e) as f32

}

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn test_i2f() {

        assert_eq!(f_i2f(0.), 0.);

        assert_eq!(f_i2f(f32::INFINITY), f32::INFINITY);

        assert_eq!(f_i2f(f32::NEG_INFINITY), f32::INFINITY);

        assert_eq!(f_i2f(1.), 0.13574767);

        assert_eq!(f_i2f(-1.), 0.13574767);

        assert_eq!(f_i2f(9.432), 1314.6553);

        assert_eq!(f_i2f(-9.432), 1314.6553);

    }

}