pxfm 0.1.28

/*
 * // 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());
    }
}