pxfm 0.1.29

/*
 * // 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.
 */
use crate::bessel::j0f_coeffs::{J0_ZEROS, J0_ZEROS_VALUE, J0F_COEFFS};
use crate::bessel::trigo_bessel::cos_small;
use crate::double_double::DoubleDouble;
use crate::polyeval::{f_polyeval9, f_polyeval10, f_polyeval12, f_polyeval14};
use crate::rounding::CpuCeil;
use crate::sincos_reduce::rem2pif_any;

/// Bessel of the first kind of order 0
///
/// Max ulp 0.5
pub fn f_j0f(x: f32) -> f32 {
    let ux = x.to_bits().wrapping_shl(1);
    if ux >= 0xffu32 << 24 || ux <= 0x6800_0000u32 {
        // |x| == 0, |x| == inf, |x| == NaN, |x| <= f32::EPSILON
        if ux == 0 {
            // |x| == 0
            return f64::from_bits(0x3ff0000000000000) as f32;
        }
        if x.is_infinite() {
            return 0.;
        }

        if ux <= 0x6800_0000u32 {
            // |x| < f32::EPSILON
            // taylor series for J0(x) ~ 1 - x^2/4 + O(x^4)
            #[cfg(any(
                all(
                    any(target_arch = "x86", target_arch = "x86_64"),
                    target_feature = "fma"
                ),
                target_arch = "aarch64"
            ))]
            {
                use crate::common::f_fmlaf;
                return f_fmlaf(x, -x * 0.25, 1.);
            }
            #[cfg(not(any(
                all(
                    any(target_arch = "x86", target_arch = "x86_64"),
                    target_feature = "fma"
                ),
                target_arch = "aarch64"
            )))]
            {
                use crate::common::f_fmla;
                let dx = x as f64;
                return f_fmla(dx, -dx * 0.25, 1.) as f32;
            }
        }

        return x + f32::NAN; // x == NaN
    }

    let x_abs = x.to_bits() & 0x7fff_ffff;

    if x_abs <= 0x4295999au32 {
        // |x| <= 74.8
        if x_abs <= 0x3e800000u32 {
            // |x| <= 0.25
            return j0f_maclaurin_series(x);
        }

        if x_abs == 0x401a42e8u32 {
            return f32::from_bits(0xbb3b2f69u32);
        }

        return small_argument_path(x);
    }

    // Exceptions
    if x_abs == 0x65ce46e4 {
        return f32::from_bits(0x1eed85c4);
    } else if x_abs == 0x7e3dcda0 {
        return f32::from_bits(0x92b81111);
    } else if x_abs == 0x76d84625 {
        return f32::from_bits(0x95d7a68b);
    } else if x_abs == 0x6bf68a7b {
        return f32::from_bits(0x1dc70a09);
    } else if x_abs == 0x7842c820 {
        return f32::from_bits(0x17ebf13e);
    } else if x_abs == 0x4ba332e9 {
        return f32::from_bits(0x27250206);
    }

    j0f_asympt(f32::from_bits(x_abs))
}

/**
Generated by SageMath:
```python
# Maclaurin series for j0
def print_expansion_at_0_f():
    print(f"pub(crate) const J0_MACLAURIN_SERIES: [u64; 9] = [")
    from mpmath import mp, j0, taylor
    mp.prec = 60
    poly = taylor(lambda val: j0(val), 0, 18)
    z = 0
    for i in range(0, 18, 2):
        print(f"{double_to_hex(poly[i])},")
    print("];")

    print(f"poly {poly}")

print_expansion_at_0_f()
```
**/
#[inline]
fn j0f_maclaurin_series(x: f32) -> f32 {
    pub(crate) const C: [u64; 9] = [
        0x3ff0000000000000,
        0xbfd0000000000000,
        0x3f90000000000000,
        0xbf3c71c71c71c71c,
        0x3edc71c71c71c71c,
        0xbe723456789abcdf,
        0x3e002e85c0898b71,
        0xbd8522a43f65486a,
        0x3d0522a43f65486a,
    ];
    let dx = x as f64;
    f_polyeval9(
        dx * dx,
        f64::from_bits(C[0]),
        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]),
    ) as f32
}

/// 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]
fn small_argument_path(x: f32) -> f32 {
    let x_abs = f32::from_bits(x.to_bits() & 0x7fff_ffff) as f64;

    // let avg_step = 74.6145 / 47.0;
    // let inv_step = 1.0 / avg_step;

    const INV_STEP: f64 = 0.6299043751549631;

    let fx = x_abs * INV_STEP;
    const J0_ZEROS_COUNT: f64 = (J0_ZEROS.len() - 1) as f64;
    let idx0 = unsafe { fx.min(J0_ZEROS_COUNT).to_int_unchecked::<usize>() };
    let idx1 = unsafe {
        fx.cpu_ceil()
            .min(J0_ZEROS_COUNT)
            .to_int_unchecked::<usize>()
    };

    let found_zero0 = DoubleDouble::from_bit_pair(J0_ZEROS[idx0]);
    let found_zero1 = DoubleDouble::from_bit_pair(J0_ZEROS[idx1]);

    let dist0 = (found_zero0.hi - x_abs).abs();
    let dist1 = (found_zero1.hi - x_abs).abs();

    let (found_zero, idx, dist) = if dist0 < dist1 {
        (found_zero0, idx0, dist0)
    } else {
        (found_zero1, idx1, dist1)
    };

    if idx == 0 {
        return j0f_maclaurin_series(x);
    }

    // We hit exact zero, value, better to return it directly
    if dist == 0. {
        return f64::from_bits(J0_ZEROS_VALUE[idx]) as f32;
    }

    let c = &J0F_COEFFS[idx - 1];

    let r = (x_abs - found_zero.hi) - found_zero.lo;

    let p = f_polyeval14(
        r,
        f64::from_bits(c[0]),
        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]),
        f64::from_bits(c[10]),
        f64::from_bits(c[11]),
        f64::from_bits(c[12]),
        f64::from_bits(c[13]),
    );

    p as f32
}

#[inline]
pub(crate) fn j1f_rsqrt(x: f64) -> f64 {
    (1. / x) * x.sqrt()
}

/*
   Evaluates:
   J1 = sqrt(2/(PI*x)) * beta(x) * cos(x - PI/4 - alpha(x))
*/
#[inline]
fn j0f_asympt(x: f32) -> f32 {
    let dx = x as f64;

    let alpha = j0f_asympt_alpha(dx);
    let beta = j0f_asympt_beta(dx);

    let angle = rem2pif_any(x);

    const SQRT_2_OVER_PI: f64 = f64::from_bits(0x3fe9884533d43651);
    const MPI_OVER_4: f64 = f64::from_bits(0xbfe921fb54442d18);

    let x0pi34 = MPI_OVER_4 - alpha;
    let r0 = angle + x0pi34;

    let m_cos = cos_small(r0);

    let z0 = beta * m_cos;
    let scale = SQRT_2_OVER_PI * j1f_rsqrt(dx);

    (scale * z0) as f32
}

/**
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(0, x, 50)
Q = Qn_asymptotic(0, 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)
```
**/
#[inline]
pub(crate) fn j0f_asympt_alpha(x: f64) -> f64 {
    const C: [u64; 12] = [
        0x3fc0000000000000,
        0xbfb0aaaaaaaaaaab,
        0x3fcad33333333333,
        0xbffa358492492492,
        0x403779a1f8e38e39,
        0xc080bd1fc8b1745d,
        0x40d16b51e66c789e,
        0xc128ecc3af33ab37,
        0x418779dae2b8512f,
        0xc1ec296336955c7f,
        0x4254f5ee683b6432,
        0xc2c2f51eced6693f,
    ];
    let recip = 1. / x;
    let x2 = recip * recip;
    let p = f_polyeval12(
        x2,
        f64::from_bits(C[0]),
        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]),
        f64::from_bits(C[10]),
        f64::from_bits(C[11]),
    );
    p * recip
}

/**
Beta series

Generated by SageMath:
```python
#generate b series
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(0, x, 50)
Q = Qn_asymptotic(0, 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 series
print(b_series)
```
**/
#[inline]
pub(crate) fn j0f_asympt_beta(x: f64) -> f64 {
    const C: [u64; 10] = [
        0x3ff0000000000000,
        0xbfb0000000000000,
        0x3fba800000000000,
        0xbfe15f0000000000,
        0x4017651180000000,
        0xc05ab8c13b800000,
        0x40a730492f262000,
        0xc0fc73a7acd696f0,
        0x41577458dd9fce68,
        0xc1b903ab9b27e18f,
    ];
    let recip = 1. / x;
    let x2 = recip * recip;
    f_polyeval10(
        x2,
        f64::from_bits(C[0]),
        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]),
    )
}

#[cfg(test)]
mod tests {
    use crate::f_j0f;

    #[test]
    fn test_j0f() {
        println!("0x{:8x}", f32::EPSILON.to_bits().wrapping_shl(1));
        assert_eq!(f_j0f(-3123.), 0.012329336);
        assert_eq!(f_j0f(-0.1), 0.99750155);
        assert_eq!(f_j0f(-15.1), -0.03456193);
        assert_eq!(f_j0f(3123.), 0.012329336);
        assert_eq!(f_j0f(0.1), 0.99750155);
        assert_eq!(f_j0f(15.1), -0.03456193);
        assert_eq!(f_j0f(f32::INFINITY), 0.);
        assert_eq!(f_j0f(f32::NEG_INFINITY), 0.);
        assert!(f_j0f(f32::NAN).is_nan());
    }
}