/*

 * // 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::common::f_fmla;

use crate::exponents::core_expf;

use crate::logs::fast_logf;

use crate::polyeval::{f_estrin_polyeval5, f_estrin_polyeval8, f_polyeval4, f_polyeval11};

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

///

/// ulp 0.5

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

    let ux = x.to_bits();

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

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

        if ux.wrapping_shl(1) == 0 {

            // |x| == 0

            return f32::INFINITY;

        }

        if x.is_infinite() {

            return if x.is_sign_positive() { 0. } else { f32::NAN };

        }

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

    }

    let xb = x.to_bits();

    if xb >= 0x42cbceefu32 {

        // |x| >= 101.90417

        return 0.;

    }

    if xb <= 0x3f800000u32 {

        // x <= 1.0

        if xb <= 0x3e9eb852u32 {

            // x <= 0.31

            if xb <= 0x34000000u32 {

                // x <= f32::EPSILON

                let dx = x as f64;

                let r = 2. / (dx * dx);

                return r as f32;

            }

            return k2f_tiny(x);

        }

        return k2f_small(x);

    }

    k2f_asympt(x)

}

#[inline]

fn k2f_tiny(x: f32) -> f32 {

    // Power series at zero for K2

    // 2.0000000000000000/x^2-0.50000000000000000-0.12500000000000000 (-0.86593151565841245+1.0000000000000000 Log[x]) x^2-0.010416666666666667 (-1.5325981823250791+1.0000000000000000 Log[x]) x^4-0.00032552083333333333 (-1.9075981823250791+1.0000000000000000 Log[x]) x^6-0.0000054253472222222222 (-2.1742648489917458+1.0000000000000000 Log[x]) x^8+O[x]^9

    //-0.50000000000000000+2.0000000000000000/x^2 + a3 * x^8 + x^6 * a2 + x^4 * a1 + x^2 * a0

    let dx = x as f64;

    let log_x = fast_logf(x);

    let a0 = f_fmla(-4.0000000000000000, log_x, 3.4637260626336498) * 0.031250000000000000;

    let a1 = f_fmla(-12.000000000000000, log_x, 18.391178187900949) * 0.00086805555555555556;

    let a2 = f_fmla(-24.000000000000000, log_x, 45.782356375801899) * 0.000013563368055555556;

    let a3 = (log_x - 2.1742648489917458) * (-0.0000054253472222222222);

    let dx_sqr = dx * dx;

    let two_over_dx = 2. / dx_sqr;

    let p = f_polyeval4(dx_sqr, a0, a1, a2, a3);

    let r = f_fmla(p, dx_sqr, two_over_dx) - 0.5;

    r as f32

}

/**

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.3,1},4,4},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) -> f64 {

    let dx = x as f64;

    let x_sqr = dx * dx;

    let p_num = f_estrin_polyeval5(

        x_sqr,

        f64::from_bits(0x3fc0000000000000),

        f64::from_bits(0x3f81520c0669099e),

        f64::from_bits(0x3f27310bf5c5e9b0),

        f64::from_bits(0x3eb8e2947e0a6098),

        f64::from_bits(0x3e336dfad46e2f35),

    );

    let p_den = f_estrin_polyeval5(

        x_sqr,

        f64::from_bits(0x3ff0000000000000),

        f64::from_bits(0xbf900d253bb12edc),

        f64::from_bits(0x3f1ed3d9ab228297),

        f64::from_bits(0xbea14e6660c00303),

        f64::from_bits(0x3e13eb951a6cf38f),

    );

    let p = p_num / p_den;

    p * x_sqr

}

/**

Series for

R(x^2) := (BesselK(2, x) - Log(x)*BesselI(2, x) - 2/x^2)/(1+x^2)

Generated by Wolfram Mathematica:

```text

<<FunctionApproximations`

ClearAll["Global`*"]

f[x_]:=(BesselK[2,x]-Log[x]BesselI[2,x]-2/(x^2))/(1+x^2)

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

{err,approx}=MiniMaxApproximation[g[z],{z,{0.3,1.0},10,10},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]

fn k2f_small(x: f32) -> f32 {

    let dx = x as f64;

    let dx_sqr = dx * dx;

    let p_num = f_polyeval11(

        dx_sqr,

        f64::from_bits(0xbfdff794c9ee3b5c),

        f64::from_bits(0xc047d3276f18e5d2),

        f64::from_bits(0xc09200ed3702875a),

        f64::from_bits(0xc0c39f395c47be27),

        f64::from_bits(0xc0e0ec95bd1a3192),

        f64::from_bits(0xc0e5973cb871c8d0),

        f64::from_bits(0xc0cdaf528de00d53),

        f64::from_bits(0xc0afe6d3009de17c),

        f64::from_bits(0xc098417b22844112),

        f64::from_bits(0x4025c45260bb1b6a),

        f64::from_bits(0x402f2bf6b95ffe0c),

    );

    let p_den = f_polyeval11(

        dx_sqr,

        f64::from_bits(0x3ff0000000000000),

        f64::from_bits(0x405879a43b253224),

        f64::from_bits(0x40a3a501408a0198),

        f64::from_bits(0x40d8172abc4a8ccc),

        f64::from_bits(0x40f9fcb05e98bdbd),

        f64::from_bits(0x4109c45b54be586b),

        f64::from_bits(0x4106ad7023dd0b90),

        f64::from_bits(0x40ed7e988d2ba5a9),

        f64::from_bits(0x40966305e1c1123a),

        f64::from_bits(0xc090832b6a87317c),

        f64::from_bits(0x403b48eb703f4644),

    );

    let p = p_num / p_den;

    let two_over_dx_sqr = 2. / dx_sqr;

    let lg = fast_logf(x);

    let v_i = i2f_small(x);

    let z = f_fmla(lg, v_i, two_over_dx_sqr);

    let z0 = f_fmla(p, f_fmla(dx, dx, 1.), z);

    z0 as f32

}

/**

Generated by Wolfram Mathematica:

```text

<<FunctionApproximations`

ClearAll["Global`*"]

f[x_]:=Sqrt[x] Exp[x] BesselK[1,x]

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

{err, approx}=MiniMaxApproximation[g[z],{z,{0.000000001,1},7,7},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]

fn k2f_asympt(x: f32) -> f32 {

    let dx = x as f64;

    let recip = 1. / dx;

    let e = core_expf(x);

    let r_sqrt = dx.sqrt();

    let p_num = f_estrin_polyeval8(

        recip,

        f64::from_bits(0x3ff40d931ff626f2),

        f64::from_bits(0x402d954dceb445df),

        f64::from_bits(0x405084ea6680d028),

        f64::from_bits(0x406242344a8ea488),

        f64::from_bits(0x406594aa56f50fea),

        f64::from_bits(0x405aa04eb4f0af1c),

        f64::from_bits(0x403dd3e8e63849ef),

        f64::from_bits(0x4004e85453648d43),

    );

    let p_den = f_estrin_polyeval8(

        recip,

        f64::from_bits(0x3ff0000000000000),

        f64::from_bits(0x4023da9f4e05358e),

        f64::from_bits(0x4040a4e4ceb523c9),

        f64::from_bits(0x404725c423c9f990),

        f64::from_bits(0x403a60c00deededc),

        f64::from_bits(0x40149975b84c3946),

        f64::from_bits(0x3fc69439846db871),

        f64::from_bits(0xbf6400819bac6f45),

    );

    let v = p_num / p_den;

    let pp = v / (e * r_sqrt);

    pp as f32

}

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn test_k2f() {

        assert!(f_k2f(-1.).is_nan());

        assert!(f_k2f(f32::NAN).is_nan());

        assert_eq!(f_k2f(0.), f32::INFINITY);

        assert_eq!(f_k2f(0.65), 4.3059196);

        assert_eq!(f_k2f(1.65), 0.44830766);

    }

}