/*

 * // Copyright (c) Radzivon Bartoshyk 7/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.

 */

#![allow(clippy::excessive_precision)]

use crate::bessel::alpha1::bessel_1_asympt_alpha_fast;

use crate::bessel::beta1::bessel_1_asympt_beta_fast;

use crate::bessel::j1_coeffs::{J1_COEFFS, J1_ZEROS, J1_ZEROS_VALUE};

use crate::bessel::j1_coeffs_taylor::J1_COEFFS_TAYLOR;

use crate::common::f_fmla;

use crate::double_double::DoubleDouble;

use crate::polyeval::{f_polyeval9, f_polyeval19};

use crate::rounding::CpuCeil;

use crate::rounding::CpuRound;

use crate::sin_helper::sin_dd_small_fast;

/// Normalized jinc 2*J1(PI\*x)/(pi\*x)

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

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

    if ux >= 0x7ffu64 << 53 || ux <= 0x7960000000000000u64 {

        // |x| <= f64::EPSILON, |x| == inf, x == NaN

        if ux <= 0x7960000000000000u64 {

            // |x| <= f64::EPSILON

            return 1.0;

        }

        if x.is_infinite() {

            return 0.;

        }

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

    }

    let ax: u64 = x.to_bits() & 0x7fff_ffff_ffff_ffff;

    if ax < 0x4052a6784230fcf8u64 {

        // |x| < 74.60109

        if ax < 0x3fd3333333333333 {

            // |x| < 0.3

            return jincpi_near_zero(f64::from_bits(ax));

        }

        let scaled_pix = f64::from_bits(ax) * std::f64::consts::PI; // just test boundaries

        if scaled_pix < 74.60109 {

            return jinc_small_argument_fast(f64::from_bits(ax));

        }

    }

    jinc_asympt_fast(f64::from_bits(ax))

}

/*

   Evaluates:

   J1 = sqrt(2/(PI*x)) * beta(x) * cos(x - 3*PI/4 - alpha(x))

   discarding 1*PI/2 using identities gives:

   J1 = sqrt(2/(PI*x)) * beta(x) * sin(x - PI/4 - alpha(x))

   to avoid squashing small (-PI/4 - alpha(x)) into a large x actual expansion is:

   J1 = sqrt(2/(PI*x)) * beta(x) * sin((x mod 2*PI) - PI/4 - alpha(x))

*/

#[inline]

fn jinc_asympt_fast(ox: f64) -> f64 {

    const PI: DoubleDouble = DoubleDouble::new(

        f64::from_bits(0x3ca1a62633145c07),

        f64::from_bits(0x400921fb54442d18),

    );

    let px = DoubleDouble::quick_mult_f64(PI, ox);

    // 2^(3/2)/(Pi^2)

    // reduce argument 2*sqrt(2/(PI*(x*PI))) = 2*sqrt(2)/PI

    // adding additional pi from division then 2*sqrt(2)/PI^2

    const Z2_3_2_OVER_PI_SQR: DoubleDouble =

        DoubleDouble::from_bit_pair((0xbc76213a285b8094, 0x3fd25751e5614413));

    const MPI_OVER_4: DoubleDouble = DoubleDouble::new(

        f64::from_bits(0xbc81a62633145c07),

        f64::from_bits(0xbfe921fb54442d18),

    );

    // argument reduction assuming x here value is already multiple of PI.

    // k = round((x*Pi) / (pi*2))

    let kd = (ox * 0.5).cpu_round();

    //  y = (x * Pi) - k * 2

    let rem = f_fmla(kd, -2., ox);

    let angle = DoubleDouble::quick_mult_f64(PI, rem);

    let recip = px.recip();

    let alpha = bessel_1_asympt_alpha_fast(recip);

    let beta = bessel_1_asympt_beta_fast(recip);

    // Without full subtraction cancellation happens sometimes

    let x0pi34 = DoubleDouble::full_dd_sub(MPI_OVER_4, alpha);

    let r0 = DoubleDouble::full_dd_add(angle, x0pi34);

    let m_sin = sin_dd_small_fast(r0);

    let z0 = DoubleDouble::quick_mult(beta, m_sin);

    let ox_recip = DoubleDouble::from_quick_recip(ox);

    let dx_sqrt = ox_recip.fast_sqrt();

    let scale = DoubleDouble::quick_mult(Z2_3_2_OVER_PI_SQR, dx_sqrt);

    let p = DoubleDouble::quick_mult(scale, z0);

    DoubleDouble::quick_mult(p, ox_recip).to_f64()

}

#[inline]

pub(crate) fn jincpi_near_zero(x: f64) -> f64 {

    // Polynomial Generated by Wolfram Mathematica:

    // <<FunctionApproximations`

    // ClearAll["Global`*"]

    // f[x_]:=2*BesselJ[1,x*Pi]/(x*Pi)

    // {err,approx}=MiniMaxApproximation[f[z],{z,{2^-23,0.3},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}]]

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

        (0xbb3bddffe9450ca6, 0x3ff0000000000000),

        (0x3c4b0b0a7393eccb, 0xbfde4cd3c3c87615),

        (0xbc8f9f784e0594a6, 0xbff043283b1e383f),

        (0xbc7af77bca466875, 0x3fdee46673cf919f),

        (0xbc1b62837b038ea8, 0x3fd0b7cc55c9a4af),

        (0x3c6c08841871f124, 0xbfc002b65231dcdd),

        (0xbc36cf2d89ea63bc, 0xbf949022a7a0712b),

        (0xbbf535d492c0ac1c, 0x3f840b48910d5105),

    ];

    const Q: [(u64, u64); 8] = [

        (0x0000000000000000, 0x3ff0000000000000),

        (0x3c4aba6577f3253e, 0xbfde4cd3c3c87615),

        (0x3c52f58f82e3438c, 0x3fcbd0a475006cf9),

        (0x3c36e496237d6b49, 0xbfb9f4cea13b06e9),

        (0xbbbbf3e4ef3a28fe, 0x3f967ed0cee85392),

        (0x3c267ac442bb3bcf, 0xbf846e192e22f862),

        (0x3bd84e9888993cb0, 0x3f51e0fff3cfddee),

        (0x3bd7c0285797bd8e, 0xbf3ea7a621fa1c8c),

    ];

    let x2 = DoubleDouble::from_exact_mult(x, x);

    let x4 = x2 * x2;

    let p0 = DoubleDouble::mul_f64_add(

        DoubleDouble::from_bit_pair(P[1]),

        x,

        DoubleDouble::from_bit_pair(P[0]),

    );

    let p1 = DoubleDouble::mul_f64_add(

        DoubleDouble::from_bit_pair(P[3]),

        x,

        DoubleDouble::from_bit_pair(P[2]),

    );

    let p2 = DoubleDouble::mul_f64_add(

        DoubleDouble::from_bit_pair(P[5]),

        x,

        DoubleDouble::from_bit_pair(P[4]),

    );

    let p3 = DoubleDouble::mul_f64_add(

        DoubleDouble::from_bit_pair(P[7]),

        x,

        DoubleDouble::from_bit_pair(P[6]),

    );

    let q0 = DoubleDouble::mul_add(x2, p1, p0);

    let q1 = DoubleDouble::mul_add(x2, p3, p2);

    let p_num = DoubleDouble::mul_add(x4, q1, q0);

    let p0 = DoubleDouble::mul_f64_add(

        DoubleDouble::from_bit_pair(Q[1]),

        x,

        DoubleDouble::from_bit_pair(Q[0]),

    );

    let p1 = DoubleDouble::mul_f64_add(

        DoubleDouble::from_bit_pair(Q[3]),

        x,

        DoubleDouble::from_bit_pair(Q[2]),

    );

    let p2 = DoubleDouble::mul_f64_add(

        DoubleDouble::from_bit_pair(Q[5]),

        x,

        DoubleDouble::from_bit_pair(Q[4]),

    );

    let p3 = DoubleDouble::mul_f64_add(

        DoubleDouble::from_bit_pair(Q[7]),

        x,

        DoubleDouble::from_bit_pair(Q[6]),

    );

    let q0 = DoubleDouble::mul_add(x2, p1, p0);

    let q1 = DoubleDouble::mul_add(x2, p3, p2);

    let p_den = DoubleDouble::mul_add(x4, q1, q0);

    DoubleDouble::div(p_num, p_den).to_f64()

}

/// This method on small range searches for nearest zero or extremum.

/// Then picks stored series expansion at the point end evaluates the poly at the point.

#[inline]

pub(crate) fn jinc_small_argument_fast(x: f64) -> f64 {

    const PI: DoubleDouble = DoubleDouble::new(

        f64::from_bits(0x3ca1a62633145c07),

        f64::from_bits(0x400921fb54442d18),

    );

    // let avg_step = 74.60109 / 47.0;

    // let inv_step = 1.0 / avg_step;

    let dx = DoubleDouble::quick_mult_f64(PI, x);

    const INV_STEP: f64 = 0.6300176043004198;

    let fx = dx.hi * INV_STEP;

    const J1_ZEROS_COUNT: f64 = (J1_ZEROS.len() - 1) as f64;

    let idx0 = unsafe { fx.min(J1_ZEROS_COUNT).to_int_unchecked::<usize>() };

    let idx1 = unsafe {

        fx.cpu_ceil()

            .min(J1_ZEROS_COUNT)

            .to_int_unchecked::<usize>()

    };

    let found_zero0 = DoubleDouble::from_bit_pair(J1_ZEROS[idx0]);

    let found_zero1 = DoubleDouble::from_bit_pair(J1_ZEROS[idx1]);

    let dist0 = (found_zero0.hi - dx.hi).abs();

    let dist1 = (found_zero1.hi - dx.hi).abs();

    let (found_zero, idx, dist) = if dist0 < dist1 {

        (found_zero0, idx0, dist0)

    } else {

        (found_zero1, idx1, dist1)

    };

    if idx == 0 {

        return jincpi_near_zero(x);

    }

    let r = DoubleDouble::quick_dd_sub(dx, found_zero);

    // We hit exact zero, value, better to return it directly

    if dist == 0. {

        return DoubleDouble::quick_mult_f64(

            DoubleDouble::from_f64_div_dd(f64::from_bits(J1_ZEROS_VALUE[idx]), dx),

            2.,

        )

        .to_f64();

    }

    let is_zero_too_close = dist.abs() < 1e-3;

    let c = if is_zero_too_close {

        &J1_COEFFS_TAYLOR[idx - 1]

    } else {

        &J1_COEFFS[idx - 1]

    };

    let p = f_polyeval19(

        r.hi,

        f64::from_bits(c[5].1),

        f64::from_bits(c[6].1),

        f64::from_bits(c[7].1),

        f64::from_bits(c[8].1),

        f64::from_bits(c[9].1),

        f64::from_bits(c[10].1),

        f64::from_bits(c[11].1),

        f64::from_bits(c[12].1),

        f64::from_bits(c[13].1),

        f64::from_bits(c[14].1),

        f64::from_bits(c[15].1),

        f64::from_bits(c[16].1),

        f64::from_bits(c[17].1),

        f64::from_bits(c[18].1),

        f64::from_bits(c[19].1),

        f64::from_bits(c[20].1),

        f64::from_bits(c[21].1),

        f64::from_bits(c[22].1),

        f64::from_bits(c[23].1),

    );

    let mut z = DoubleDouble::mul_f64_add(r, p, DoubleDouble::from_bit_pair(c[4]));

    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[3]));

    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[2]));

    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[1]));

    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[0]));

    z = DoubleDouble::div(z, dx);

    z.hi *= 2.;

    z.lo *= 2.;

    let err = f_fmla(

        z.hi,

        f64::from_bits(0x3c70000000000000), // 2^-56

        f64::from_bits(0x3bf0000000000000), // 2^-64

    );

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

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

    if ub == lb {

        return z.to_f64();

    }

    j1_small_argument_dd(r, c, dx)

}

fn j1_small_argument_dd(r: DoubleDouble, c0: &[(u64, u64); 24], inv_scale: DoubleDouble) -> f64 {

    let c = &c0[15..];

    let p0 = f_polyeval9(

        r.to_f64(),

        f64::from_bits(c[0].1),

        f64::from_bits(c[1].1),

        f64::from_bits(c[2].1),

        f64::from_bits(c[3].1),

        f64::from_bits(c[4].1),

        f64::from_bits(c[5].1),

        f64::from_bits(c[6].1),

        f64::from_bits(c[7].1),

        f64::from_bits(c[8].1),

    );

    let c = c0;

    let mut p_e = DoubleDouble::mul_f64_add(r, p0, DoubleDouble::from_bit_pair(c[14]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[13]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[12]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[11]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[10]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[9]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[8]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[7]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[6]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[5]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[4]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[3]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[2]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[1]));

    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[0]));

    let p = DoubleDouble::from_exact_add(p_e.hi, p_e.lo);

    let mut z = DoubleDouble::div(p, inv_scale);

    z.hi *= 2.;

    z.lo *= 2.;

    z.to_f64()

}

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn test_jincpi() {

        assert_eq!(f_jincpi(f64::EPSILON), 1.0);

        assert_eq!(f_jincpi(0.000043242121), 0.9999999976931268);

        assert_eq!(f_jincpi(0.5000000000020244), 0.7217028449014163);

        assert_eq!(f_jincpi(73.81695991658546), -0.0004417546638317049);

        assert_eq!(f_jincpi(0.01), 0.9998766350182722);

        assert_eq!(f_jincpi(0.9), 0.28331697846510623);

        assert_eq!(f_jincpi(3.831705970207517), -0.036684415010255086);

        assert_eq!(f_jincpi(-3.831705970207517), -0.036684415010255086);

        assert_eq!(

            f_jincpi(0.000000000000000000000000000000000000008827127),

            1.0

        );

        assert_eq!(

            f_jincpi(-0.000000000000000000000000000000000000008827127),

            1.0

        );

        assert_eq!(f_jincpi(5.4), -0.010821736808448256);

        assert_eq!(

            f_jincpi(77.743162408196766932633181568235159),

            -0.00041799098646950523

        );

        assert_eq!(

            f_jincpi(84.027189586293545175976760219782591),

            -0.00023927934929850555

        );

        assert_eq!(f_jincpi(f64::NEG_INFINITY), 0.0);

        assert_eq!(f_jincpi(f64::INFINITY), 0.0);

        assert!(f_jincpi(f64::NAN).is_nan());

    }

}