/*

 * // 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::double_double::DoubleDouble;

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

use crate::polyeval::f_polyeval9;

/**

Note expansion generation below: this is negative series expressed in Sage as positive,

so before any real evaluation `x=1/x` should be applied

Generated by SageMath:

```python

def binomial_like(n, m):

    prod = QQ(1)

    z = QQ(4)*(n**2)

    for k in range(1,m + 1):

        prod *= (z - (2*k - 1)**2)

    return prod / (QQ(2)**(2*m) * (ZZ(m).factorial()))

R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)

x = R.gen()

def Pn_asymptotic(n, y, terms=10):

    # now y = 1/x

    return sum( (-1)**m * binomial_like(n, 2*m) / (QQ(2)**(2*m)) * y**(QQ(2)*m) for m in range(terms) )

def Qn_asymptotic(n, y, terms=10):

    return sum( (-1)**m * binomial_like(n, 2*m + 1) / (QQ(2)**(2*m + 1)) * y**(QQ(2)*m + 1) for m in range(terms) )

P = Pn_asymptotic(1, x, 50)

Q = Qn_asymptotic(1, x, 50)

def sqrt_series(s):

    val = S.valuation()

    lc = S[val]  # Leading coefficient

    b = lc.sqrt() * x**(val // 2)

    for _ in range(5):

        b = (b + S / b) / 2

        b = b

    return b

S = (P**2 + Q**2).truncate(50)

b_series = sqrt_series(S).truncate(30)

# see the beta series

print(b_series)

```

See notes/bessel_asympt.ipynb for generation

**/

#[inline]

pub(crate) fn bessel_1_asympt_beta_fast(recip: DoubleDouble) -> DoubleDouble {

    const C: [u64; 10] = [

        0x3ff0000000000000,

        0x3fc8000000000000,

        0xbfc8c00000000000,

        0x3fe9c50000000000,

        0xc01ef5b680000000,

        0x40609860dd400000,

        0xc0abae9b7a06e000,

        0x41008711d41c1428,

        0xc15ab70164c8be6e,

        0x41bc1055e24f297f,

    ];

    // Doing (1/x)*(1/x) instead (1/(x*x)) to avoid spurious overflow/underflow

    let x2 = DoubleDouble::quick_mult(recip, recip);

    let p = f_polyeval9(

        x2.hi,

        f64::from_bits(C[1]),

        f64::from_bits(C[2]),

        f64::from_bits(C[3]),

        f64::from_bits(C[4]),

        f64::from_bits(C[5]),

        f64::from_bits(C[6]),

        f64::from_bits(C[7]),

        f64::from_bits(C[8]),

        f64::from_bits(C[9]),

    );

    DoubleDouble::mul_f64_add_f64(x2, p, f64::from_bits(C[0]))

}

/**

Note expansion generation below: this is negative series expressed in Sage as positive,

so before any real evaluation `x=1/x` should be applied

Generated by SageMath:

```python

def binomial_like(n, m):

    prod = QQ(1)

    z = QQ(4)*(n**2)

    for k in range(1,m + 1):

        prod *= (z - (2*k - 1)**2)

    return prod / (QQ(2)**(2*m) * (ZZ(m).factorial()))

R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)

x = R.gen()

def Pn_asymptotic(n, y, terms=10):

    # now y = 1/x

    return sum( (-1)**m * binomial_like(n, 2*m) / (QQ(2)**(2*m)) * y**(QQ(2)*m) for m in range(terms) )

def Qn_asymptotic(n, y, terms=10):

    return sum( (-1)**m * binomial_like(n, 2*m + 1) / (QQ(2)**(2*m + 1)) * y**(QQ(2)*m + 1) for m in range(terms) )

P = Pn_asymptotic(1, x, 50)

Q = Qn_asymptotic(1, x, 50)

def sqrt_series(s):

    val = S.valuation()

    lc = S[val]  # Leading coefficient

    b = lc.sqrt() * x**(val // 2)

    for _ in range(5):

        b = (b + S / b) / 2

        b = b

    return b

S = (P**2 + Q**2).truncate(50)

b_series = sqrt_series(S).truncate(30)

# see the beta series

print(b_series)

```

See notes/bessel_asympt.ipynb for generation

**/

#[inline]

pub(crate) fn bessel_1_asympt_beta(recip: DoubleDouble) -> DoubleDouble {

    const C: [(u64, u64); 10] = [

        (0x0000000000000000, 0x3ff0000000000000), // 1

        (0x0000000000000000, 0x3fc8000000000000), // 2

        (0x0000000000000000, 0xbfc8c00000000000), // 3

        (0x0000000000000000, 0x3fe9c50000000000), // 4

        (0x0000000000000000, 0xc01ef5b680000000), // 5

        (0x0000000000000000, 0x40609860dd400000), // 6

        (0x0000000000000000, 0xc0abae9b7a06e000), // 7

        (0x0000000000000000, 0x41008711d41c1428), // 8

        (0xbdf7a00000000000, 0xc15ab70164c8be6e),

        (0xbe40e1f000000000, 0x41bc1055e24f297f),

    ];

    // Doing (1/x)*(1/x) instead (1/(x*x)) to avoid spurious overflow/underflow

    let x2 = DoubleDouble::quick_mult(recip, recip);

    let mut p = DoubleDouble::mul_add(

        x2,

        DoubleDouble::from_bit_pair(C[9]),

        DoubleDouble::from_bit_pair(C[8]),

    );

    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[7].1)); // 8

    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[6].1)); // 7

    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[5].1)); // 6

    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[4].1)); // 5

    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[3].1)); // 4

    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[2].1)); // 3

    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[1].1)); // 2

    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[0].1)); // 1

    p

}

/// see [bessel_1_asympt_beta] for more info

pub(crate) fn bessel_1_asympt_beta_hard(recip: DyadicFloat128) -> DyadicFloat128 {

    static C: [DyadicFloat128; 12] = [

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -127,

            mantissa: 0x80000000_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -130,

            mantissa: 0xc0000000_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -130,

            mantissa: 0xc6000000_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -128,

            mantissa: 0xce280000_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -125,

            mantissa: 0xf7adb400_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -120,

            mantissa: 0x84c306ea_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -116,

            mantissa: 0xdd74dbd0_37000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -110,

            mantissa: 0x84388ea0_e0a14000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -105,

            mantissa: 0xd5b80b26_45f372f4_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -99,

            mantissa: 0xe082af12_794bf6f1_e1000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -92,

            mantissa: 0x94a06149_f30146bc_fe8ed000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -86,

            mantissa: 0xf212edfc_42a62526_4fac2b0c_00000000_u128,

        },

    ];

    let x2 = recip * recip;

    let mut p = C[11];

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

        p = x2 * p + C[i];

    }

    p

}