/*

 * // 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::bessel::k1::k1_small;

use crate::double_double::DoubleDouble;

use crate::dyadic_float::{DyadicFloat128, DyadicSign};

/// Modified exponentially scaled Bessel of the second kind of order 1

///

/// Computes K1(x)exp(x)

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

    let ix = x.to_bits();

    if ix >= 0x7ffu64 << 52 || ix == 0 {

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

        if ix.wrapping_shl(1) == 0 {

            // |x| == 0

            return f64::INFINITY;

        }

        if x.is_infinite() {

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

        }

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

    }

    let xb = x.to_bits();

    if xb <= 0x3ff0000000000000 {

        // x <= 1

        let v_exp = i0_exp(x);

        let v_k = k1_small(x);

        return DoubleDouble::quick_mult(v_exp, v_k).to_f64();

    }

    k1e_asympt(x)

}

/**

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.0000000000001,1},11,11},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 k1e_asympt(x: f64) -> f64 {

    let recip = DoubleDouble::from_quick_recip(x);

    let r_sqrt = DoubleDouble::from_sqrt(x);

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

        (0xbc9a6a0690becb3b, 0x3ff40d931ff62706),

        (0xbce573e1bbf2f0b7, 0x40402cebfab5721d),

        (0x3d11a739b7c11e7b, 0x4074f58abc0cfbf1),

        (0xbd2682a09ded0116, 0x409c8315f8facef2),

        (0xbd3a19e91a120168, 0x40b65f7a4caed8b9),

        (0x3d449c3d2b834543, 0x40c4fe41fdb4e7b8),

        (0xbd6bdd415ac7f7e1, 0x40c7aa402d035d03),

        (0x3d528412ff0d6b24, 0x40bf68faddd7d850),

        (0xbd48f4bb3f61dac6, 0x40a75f5650249952),

        (0xbd1dc534b275e309, 0x4081bddd259c0582),

        (0xbcce5103350bd226, 0x4046c7a049014484),

        (0x3c8935f8acd6c1d0, 0x3fef7524082b1859),

    ];

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

        (0x3cc0d2508437b3f4, 0x40396ff483adec14),

        (0xbd130a9c9f8a5338, 0x4070225588d8c15d),

        (0xbceceba8fa0e65a2, 0x4095481f6684e3bb),

        (0x3d4099f3c178fd2a, 0x40afedc8a778bf42),

        (0xbd3a7e6a6276a3e7, 0x40bc0c060112692e),

        (0x3d11538c155b16d8, 0x40bcb12bd1101782),

        (0xbd5f7b04cdea2c61, 0x40b07fa363202e10),

        (0xbce444ed035b66c6, 0x4093d6fe8f44f838),

        (0xbcf6f88fb942b610, 0x4065c99fa44030c3),

        (0xbcbd1d2aedee5bc9, 0x40207ffabeb00eea),

        (0xbc39a0c8091102c9, 0x3facff3d892cd57a),

    ];

    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 = DoubleDouble::div(z, r_sqrt);

    let err = r.hi * f64::from_bits(0x3c10000000000000); // 2^-61

    let ub = r.hi + (r.lo + err);

    let lb = r.hi + (r.lo - err);

    if ub != lb {

        return k1e_asympt_hard(x);

    }

    r.to_f64()

}

/**

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.0000000000001,1},14,14},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}]]

```

**/

#[cold]

#[inline(never)]

fn k1e_asympt_hard(x: f64) -> f64 {

    static P: [DyadicFloat128; 15] = [

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -127,

            mantissa: 0xa06c98ff_b1382cb2_be5210ac_f26f25d1_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -122,

            mantissa: 0xc5f546cb_659a39d0_fafbd188_36ca05b9_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -118,

            mantissa: 0xcd0b7cfa_de158d26_7084bbe9_f1bdb66d_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -115,

            mantissa: 0xeac7be2f_957d1260_8849508a_2a5a8972_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -112,

            mantissa: 0xa4d14fec_fecc6444_4c7b0287_dad71a86_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -110,

            mantissa: 0x94e3180c_01df9932_ad2acd8b_bab59c05_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -109,

            mantissa: 0xb0de10f8_74918442_94a96368_8eaa4d0d_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -108,

            mantissa: 0x8adfea76_d6dbe5d9_46bfaf83_9341f4b5_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -108,

            mantissa: 0x8f0a4337_b69b602c_cf187222_f3a3379f_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -109,

            mantissa: 0xbd4c3ebf_c2db0fad_1b425641_cc470043_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -110,

            mantissa: 0x9b14d29f_9b97e3c8_c1a7b9d0_787f0ddb_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -112,

            mantissa: 0x93e670d2_07a553ef_a90d4895_cf1b5011_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -115,

            mantissa: 0x93e0ee0a_cb4d8910_6b4d3e37_f4f9df49_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -120,

            mantissa: 0xff0ce10d_5585abd1_e8a53a12_65131ad4_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -126,

            mantissa: 0xf020536d_822cbe51_c8de095a_03367c83_u128,

        },

    ];

    static Q: [DyadicFloat128; 15] = [

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -127,

            mantissa: 0x80000000_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -122,

            mantissa: 0x9c729dd5_4828a918_42807f58_d485a511_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -118,

            mantissa: 0x9ff6f631_0794001d_433ab0c5_d4c682a9_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -115,

            mantissa: 0xb3f81e8b_1e0e85a6_3928342e_c83088a1_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -113,

            mantissa: 0xf6b1c203_a60d4294_239ad045_2c67c224_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -111,

            mantissa: 0xd7a98b14_7a499762_abde5c38_3a5b40e4_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -110,

            mantissa: 0xf4eb8b77_a2cdc686_afd1273f_d464c8b7_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -109,

            mantissa: 0xb4c1e12a_93ee86fc_930c6f94_cfa6ac3a_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -109,

            mantissa: 0xaaeaab88_32b776b7_fdd76b0f_24349f41_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -110,

            mantissa: 0xc8ec9d61_5bf2ee9b_878b4962_4a5cee85_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -111,

            mantissa: 0x8b97bab0_3351673f_22f10d40_fd1c9ff3_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -114,

            mantissa: 0xd31cb80a_bf8cbedc_b0dcf7e7_c599f79e_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -117,

            mantissa: 0x96b354c8_69197193_ea4f608f_81943988_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -122,

            mantissa: 0x989af1bb_e48b5c44_7cd09746_f15e935a_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -130,

            mantissa: 0xb7b51326_23c29ed5_8d3dcf5a_79bd9a4f_u128,

        },

    ];

    let recip = DyadicFloat128::accurate_reciprocal(x);

    let r_sqrt = bessel_rsqrt_hard(x, recip);

    let mut p0 = P[14];

    for i in (0..14).rev() {

        p0 = recip * p0 + P[i];

    }

    let mut q0 = Q[14];

    for i in (0..14).rev() {

        q0 = recip * q0 + Q[i];

    }

    let v = p0 * q0.reciprocal();

    let r = v * r_sqrt;

    r.fast_as_f64()

}

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn test_k1() {

        assert_eq!(f_k1e(0.643), 2.253195748291852);

        assert_eq!(f_k1e(0.964), 1.6787831013451477);

        assert_eq!(f_k1e(2.964), 0.8123854795542738);

        assert_eq!(f_k1e(8.43), 0.4502184086111872);

        assert_eq!(f_k1e(16.43), 0.3161307996938612);

        assert_eq!(f_k1e(423.43), 0.06096117017402597);

        assert_eq!(f_k1e(9044.431), 0.01317914752085687);

        assert_eq!(k1e_asympt_hard(16.43), 0.3161307996938612);

        assert_eq!(k1e_asympt_hard(423.43), 0.06096117017402597);

        assert_eq!(k1e_asympt_hard(9044.431), 0.01317914752085687);

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

        assert_eq!(f_k1e(-0.), f64::INFINITY);

        assert!(f_k1e(-0.5).is_nan());

        assert!(f_k1e(f64::NEG_INFINITY).is_nan());

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

    }

}