/*

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

R_series = (-Q/P)

# alpha is atan(R_series) so we're doing Taylor series atan expansion on R_series

arctan_series_Z = sum([QQ(-1)**k * x**(QQ(2)*k+1) / RealField(700)(RealField(700)(2)*k+1) for k in range(25)])

alpha_series = arctan_series_Z(R_series)

# see the series

print(alpha_series)

```

See notes/bessel_asympt.ipynb for generation

**/

#[inline]

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

    const C: [u64; 12] = [

        0xbfd8000000000000,

        0x3fc5000000000000,

        0xbfd7bccccccccccd,

        0x4002f486db6db6db,

        0xc03e9fbf40000000,

        0x4084997b55945d17,

        0xc0d4a914195269d9,

        0x412cd1b53816aec1,

        0xc18aa4095d419351,

        0x41ef809305f11b9d,

        0xc2572e6809ed618b,

        0x42c4c5b6057839f9,

    ];

    // 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[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]),

        f64::from_bits(C[10]),

        f64::from_bits(C[11]),

    );

    let mut z = DoubleDouble::mul_f64_add_f64(x2, p, f64::from_bits(C[2]));

    z = DoubleDouble::mul_add_f64(x2, z, f64::from_bits(C[1]));

    z = DoubleDouble::mul_add_f64(x2, z, f64::from_bits(C[0]));

    DoubleDouble::quick_mult(z, recip)

}

/**

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)

R_series = (-Q/P)

# alpha is atan(R_series) so we're doing Taylor series atan expansion on R_series

arctan_series_Z = sum([QQ(-1)**k * x**(QQ(2)*k+1) / RealField(700)(RealField(700)(2)*k+1) for k in range(25)])

alpha_series = arctan_series_Z(R_series)

# see the series

print(alpha_series)

```

See notes/bessel_asympt.ipynb for generation

**/

#[inline]

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

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

        (0x0000000000000000, 0xbfd8000000000000),

        (0x0000000000000000, 0x3fc5000000000000),

        (0x3c6999999999999a, 0xbfd7bccccccccccd),

        (0x3cab6db6db6db6db, 0x4002f486db6db6db),

        (0x0000000000000000, 0xc03e9fbf40000000),

        (0x3d21745d1745d174, 0x4084997b55945d17),

        (0x3d789d89d89d89d9, 0xc0d4a914195269d9),

        (0xbdb999999999999a, 0x412cd1b53816aec1),

        (0xbdfe5a5a5a5a5a5a, 0xc18aa4095d419351),

        (0x3e7e0ca50d79435e, 0x41ef809305f11b9d),

        (0xbedff8b720000000, 0xc2572e6809ed618b),

        (0xbf64e5d8ae68b7a7, 0x42c4c5b6057839f9),

    ];

    // 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[11]),

        DoubleDouble::from_bit_pair(C[10]),

    );

    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[9]));

    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[8]));

    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[7]));

    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[6]));

    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[5]));

    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[4]));

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

    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[2]));

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

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

    let z = DoubleDouble::quick_mult(p, recip);

    DoubleDouble::from_exact_add(z.hi, z.lo)

}

//

/// See [bessel_1_asympt_alpha] for the info

pub(crate) fn bessel_1_asympt_alpha_hard(reciprocal: DyadicFloat128) -> DyadicFloat128 {

    static C: [DyadicFloat128; 18] = [

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -129,

            mantissa: 0xc0000000_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -130,

            mantissa: 0xa8000000_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -129,

            mantissa: 0xbde66666_66666666_66666666_66666666_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -126,

            mantissa: 0x97a436db_6db6db6d_b6db6db6_db6db6db_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -123,

            mantissa: 0xf4fdfa00_00000000_00000000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -118,

            mantissa: 0xa4cbdaac_a2e8ba2e_8ba2e8ba_2e8ba2e9_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -113,

            mantissa: 0xa548a0ca_934ec4ec_4ec4ec4e_c4ec4ec5_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -108,

            mantissa: 0xe68da9c0_b5760666_66666666_66666666_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -102,

            mantissa: 0xd5204aea_0c9a8879_69696969_69696969_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -96,

            mantissa: 0xfc04982f_88dce9e0_ca50d794_35e50d79_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -89,

            mantissa: 0xb973404f_6b0c58ff_c5b90000_00000000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -82,

            mantissa: 0xa62db02b_c1cfc563_44ea32e9_0b21642d_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -75,

            mantissa: 0xb220e7ff_443c1584_7e85f4e0_55eb851f_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -68,

            mantissa: 0xe10a255c_ca5e68cc_00c2d6c0_acdc8000_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -60,

            mantissa: 0xa573790c_5186f23b_5db502ea_d9fa5432_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -52,

            mantissa: 0x8c0ffedc_407a7015_453df84e_9c3f1d39_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Neg,

            exponent: -44,

            mantissa: 0x874226ed_c298a17a_d8c49a4e_dc9281a5_u128,

        },

        DyadicFloat128 {

            sign: DyadicSign::Pos,

            exponent: -36,

            mantissa: 0x93cab36c_9ab9495c_310fa9cd_4b065359_u128,

        },

    ];

    let x2 = reciprocal * reciprocal;

    let mut p = C[17];

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

        p = x2 * p + C[i];

    }

    p * reciprocal

}