pxfm/bessel/
j0f.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 7/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::bessel::j0f_coeffs::{J0_ZEROS, J0_ZEROS_VALUE, J0F_COEFFS};
30use crate::bessel::trigo_bessel::cos_small;
31use crate::double_double::DoubleDouble;
32use crate::polyeval::{f_polyeval9, f_polyeval10, f_polyeval12, f_polyeval14};
33use crate::sincos_reduce::rem2pif_any;
34
35/// Bessel of the first kind of order 0
36///
37/// Max ulp 0.5
38pub fn f_j0f(x: f32) -> f32 {
39    let ux = x.to_bits().wrapping_shl(1);
40    if ux >= 0xffu32 << 24 || ux <= 0x6800_0000u32 {
41        // |x| == 0, |x| == inf, |x| == NaN, |x| <= f32::EPSILON
42        if ux == 0 {
43            // |x| == 0
44            return f64::from_bits(0x3ff0000000000000) as f32;
45        }
46        if x.is_infinite() {
47            return 0.;
48        }
49
50        if ux <= 0x6800_0000u32 {
51            // |x| < f32::EPSILON
52            // taylor series for J0(x) ~ 1 - x^2/4 + O(x^4)
53            #[cfg(any(
54                all(
55                    any(target_arch = "x86", target_arch = "x86_64"),
56                    target_feature = "fma"
57                ),
58                all(target_arch = "aarch64", target_feature = "neon")
59            ))]
60            {
61                use crate::common::f_fmlaf;
62                return f_fmlaf(x, -x * 0.25, 1.);
63            }
64            #[cfg(not(any(
65                all(
66                    any(target_arch = "x86", target_arch = "x86_64"),
67                    target_feature = "fma"
68                ),
69                all(target_arch = "aarch64", target_feature = "neon")
70            )))]
71            {
72                use crate::common::f_fmla;
73                let dx = x as f64;
74                return f_fmla(dx, -dx * 0.25, 1.) as f32;
75            }
76        }
77
78        return x + f32::NAN; // x == NaN
79    }
80
81    let x_abs = x.to_bits() & 0x7fff_ffff;
82
83    if x_abs <= 0x4295999au32 {
84        // |x| <= 74.8
85        if x_abs <= 0x3e800000u32 {
86            // |x| <= 0.25
87            return j0f_maclaurin_series(x);
88        }
89
90        if x_abs == 0x401a42e8u32 {
91            return f32::from_bits(0xbb3b2f69u32);
92        }
93
94        return small_argument_path(x);
95    }
96
97    // Exceptions
98    if x_abs == 0x65ce46e4 {
99        return f32::from_bits(0x1eed85c4);
100    } else if x_abs == 0x7e3dcda0 {
101        return f32::from_bits(0x92b81111);
102    } else if x_abs == 0x76d84625 {
103        return f32::from_bits(0x95d7a68b);
104    } else if x_abs == 0x6bf68a7b {
105        return f32::from_bits(0x1dc70a09);
106    } else if x_abs == 0x7842c820 {
107        return f32::from_bits(0x17ebf13e);
108    } else if x_abs == 0x4ba332e9 {
109        return f32::from_bits(0x27250206);
110    }
111
112    j0f_asympt(f32::from_bits(x_abs))
113}
114
115/**
116Generated by SageMath:
117```python
118# Maclaurin series for j0
119def print_expansion_at_0_f():
120    print(f"pub(crate) const J0_MACLAURIN_SERIES: [u64; 9] = [")
121    from mpmath import mp, j0, taylor
122    mp.prec = 60
123    poly = taylor(lambda val: j0(val), 0, 18)
124    z = 0
125    for i in range(0, 18, 2):
126        print(f"{double_to_hex(poly[i])},")
127    print("];")
128
129    print(f"poly {poly}")
130
131print_expansion_at_0_f()
132```
133**/
134#[inline]
135fn j0f_maclaurin_series(x: f32) -> f32 {
136    pub(crate) const C: [u64; 9] = [
137        0x3ff0000000000000,
138        0xbfd0000000000000,
139        0x3f90000000000000,
140        0xbf3c71c71c71c71c,
141        0x3edc71c71c71c71c,
142        0xbe723456789abcdf,
143        0x3e002e85c0898b71,
144        0xbd8522a43f65486a,
145        0x3d0522a43f65486a,
146    ];
147    let dx = x as f64;
148    f_polyeval9(
149        dx * dx,
150        f64::from_bits(C[0]),
151        f64::from_bits(C[1]),
152        f64::from_bits(C[2]),
153        f64::from_bits(C[3]),
154        f64::from_bits(C[4]),
155        f64::from_bits(C[5]),
156        f64::from_bits(C[6]),
157        f64::from_bits(C[7]),
158        f64::from_bits(C[8]),
159    ) as f32
160}
161
162/// This method on small range searches for nearest zero or extremum.
163/// Then picks stored series expansion at the point end evaluates the poly at the point.
164#[inline]
165fn small_argument_path(x: f32) -> f32 {
166    let x_abs = f32::from_bits(x.to_bits() & 0x7fff_ffff) as f64;
167
168    // let avg_step = 74.6145 / 47.0;
169    // let inv_step = 1.0 / avg_step;
170
171    const INV_STEP: f64 = 0.6299043751549631;
172
173    let fx = x_abs * INV_STEP;
174    const J0_ZEROS_COUNT: f64 = (J0_ZEROS.len() - 1) as f64;
175    let idx0 = unsafe { fx.min(J0_ZEROS_COUNT).to_int_unchecked::<usize>() };
176    let idx1 = unsafe { fx.ceil().min(J0_ZEROS_COUNT).to_int_unchecked::<usize>() };
177
178    let found_zero0 = DoubleDouble::from_bit_pair(J0_ZEROS[idx0]);
179    let found_zero1 = DoubleDouble::from_bit_pair(J0_ZEROS[idx1]);
180
181    let dist0 = (found_zero0.hi - x_abs).abs();
182    let dist1 = (found_zero1.hi - x_abs).abs();
183
184    let (found_zero, idx, dist) = if dist0 < dist1 {
185        (found_zero0, idx0, dist0)
186    } else {
187        (found_zero1, idx1, dist1)
188    };
189
190    if idx == 0 {
191        return j0f_maclaurin_series(x);
192    }
193
194    // We hit exact zero, value, better to return it directly
195    if dist == 0. {
196        return f64::from_bits(J0_ZEROS_VALUE[idx]) as f32;
197    }
198
199    let c = &J0F_COEFFS[idx - 1];
200
201    let r = (x_abs - found_zero.hi) - found_zero.lo;
202
203    let p = f_polyeval14(
204        r,
205        f64::from_bits(c[0]),
206        f64::from_bits(c[1]),
207        f64::from_bits(c[2]),
208        f64::from_bits(c[3]),
209        f64::from_bits(c[4]),
210        f64::from_bits(c[5]),
211        f64::from_bits(c[6]),
212        f64::from_bits(c[7]),
213        f64::from_bits(c[8]),
214        f64::from_bits(c[9]),
215        f64::from_bits(c[10]),
216        f64::from_bits(c[11]),
217        f64::from_bits(c[12]),
218        f64::from_bits(c[13]),
219    );
220
221    p as f32
222}
223
224#[inline]
225pub(crate) fn j1f_rsqrt(x: f64) -> f64 {
226    (1. / x) * x.sqrt()
227}
228
229/*
230   Evaluates:
231   J1 = sqrt(2/(PI*x)) * beta(x) * cos(x - PI/4 - alpha(x))
232*/
233#[inline]
234fn j0f_asympt(x: f32) -> f32 {
235    let dx = x as f64;
236
237    let alpha = j0f_asympt_alpha(dx);
238    let beta = j0f_asympt_beta(dx);
239
240    let angle = rem2pif_any(x);
241
242    const SQRT_2_OVER_PI: f64 = f64::from_bits(0x3fe9884533d43651);
243    const MPI_OVER_4: f64 = f64::from_bits(0xbfe921fb54442d18);
244
245    let x0pi34 = MPI_OVER_4 - alpha;
246    let r0 = angle + x0pi34;
247
248    let m_cos = cos_small(r0);
249
250    let z0 = beta * m_cos;
251    let scale = SQRT_2_OVER_PI * j1f_rsqrt(dx);
252
253    (scale * z0) as f32
254}
255
256/**
257Note expansion generation below: this is negative series expressed in Sage as positive,
258so before any real evaluation `x=1/x` should be applied.
259
260Generated by SageMath:
261```python
262def binomial_like(n, m):
263    prod = QQ(1)
264    z = QQ(4)*(n**2)
265    for k in range(1,m + 1):
266        prod *= (z - (2*k - 1)**2)
267    return prod / (QQ(2)**(2*m) * (ZZ(m).factorial()))
268
269R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
270x = R.gen()
271
272def Pn_asymptotic(n, y, terms=10):
273    # now y = 1/x
274    return sum( (-1)**m * binomial_like(n, 2*m) / (QQ(2)**(2*m)) * y**(QQ(2)*m) for m in range(terms) )
275
276def Qn_asymptotic(n, y, terms=10):
277    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) )
278
279P = Pn_asymptotic(0, x, 50)
280Q = Qn_asymptotic(0, x, 50)
281
282R_series = (-Q/P)
283
284# alpha is atan(R_series) so we're doing Taylor series atan expansion on R_series
285
286arctan_series_Z = sum([QQ(-1)**k * x**(QQ(2)*k+1) / RealField(700)(RealField(700)(2)*k+1) for k in range(25)])
287alpha_series = arctan_series_Z(R_series)
288
289# see the series
290print(alpha_series)
291```
292**/
293#[inline]
294pub(crate) fn j0f_asympt_alpha(x: f64) -> f64 {
295    const C: [u64; 12] = [
296        0x3fc0000000000000,
297        0xbfb0aaaaaaaaaaab,
298        0x3fcad33333333333,
299        0xbffa358492492492,
300        0x403779a1f8e38e39,
301        0xc080bd1fc8b1745d,
302        0x40d16b51e66c789e,
303        0xc128ecc3af33ab37,
304        0x418779dae2b8512f,
305        0xc1ec296336955c7f,
306        0x4254f5ee683b6432,
307        0xc2c2f51eced6693f,
308    ];
309    let recip = 1. / x;
310    let x2 = recip * recip;
311    let p = f_polyeval12(
312        x2,
313        f64::from_bits(C[0]),
314        f64::from_bits(C[1]),
315        f64::from_bits(C[2]),
316        f64::from_bits(C[3]),
317        f64::from_bits(C[4]),
318        f64::from_bits(C[5]),
319        f64::from_bits(C[6]),
320        f64::from_bits(C[7]),
321        f64::from_bits(C[8]),
322        f64::from_bits(C[9]),
323        f64::from_bits(C[10]),
324        f64::from_bits(C[11]),
325    );
326    p * recip
327}
328
329/**
330Beta series
331
332Generated by SageMath:
333```python
334#generate b series
335def binomial_like(n, m):
336    prod = QQ(1)
337    z = QQ(4)*(n**2)
338    for k in range(1,m + 1):
339        prod *= (z - (2*k - 1)**2)
340    return prod / (QQ(2)**(2*m) * (ZZ(m).factorial()))
341
342R = LaurentSeriesRing(RealField(300), 'x', default_prec=300)
343x = R.gen()
344
345def Pn_asymptotic(n, y, terms=10):
346    # now y = 1/x
347    return sum( (-1)**m * binomial_like(n, 2*m) / (QQ(2)**(2*m)) * y**(QQ(2)*m) for m in range(terms) )
348
349def Qn_asymptotic(n, y, terms=10):
350    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) )
351
352P = Pn_asymptotic(0, x, 50)
353Q = Qn_asymptotic(0, x, 50)
354
355def sqrt_series(s):
356    val = S.valuation()
357    lc = S[val]  # Leading coefficient
358    b = lc.sqrt() * x**(val // 2)
359
360    for _ in range(5):
361        b = (b + S / b) / 2
362        b = b
363    return b
364
365S = (P**2 + Q**2).truncate(50)
366
367b_series = sqrt_series(S).truncate(30)
368#see the series
369print(b_series)
370```
371**/
372#[inline]
373pub(crate) fn j0f_asympt_beta(x: f64) -> f64 {
374    const C: [u64; 10] = [
375        0x3ff0000000000000,
376        0xbfb0000000000000,
377        0x3fba800000000000,
378        0xbfe15f0000000000,
379        0x4017651180000000,
380        0xc05ab8c13b800000,
381        0x40a730492f262000,
382        0xc0fc73a7acd696f0,
383        0x41577458dd9fce68,
384        0xc1b903ab9b27e18f,
385    ];
386    let recip = 1. / x;
387    let x2 = recip * recip;
388    f_polyeval10(
389        x2,
390        f64::from_bits(C[0]),
391        f64::from_bits(C[1]),
392        f64::from_bits(C[2]),
393        f64::from_bits(C[3]),
394        f64::from_bits(C[4]),
395        f64::from_bits(C[5]),
396        f64::from_bits(C[6]),
397        f64::from_bits(C[7]),
398        f64::from_bits(C[8]),
399        f64::from_bits(C[9]),
400    )
401}
402
403#[cfg(test)]
404mod tests {
405    use crate::f_j0f;
406
407    #[test]
408    fn test_j0f() {
409        println!("0x{:8x}", f32::EPSILON.to_bits().wrapping_shl(1));
410        assert_eq!(f_j0f(-3123.), 0.012329336);
411        assert_eq!(f_j0f(-0.1), 0.99750155);
412        assert_eq!(f_j0f(-15.1), -0.03456193);
413        assert_eq!(f_j0f(3123.), 0.012329336);
414        assert_eq!(f_j0f(0.1), 0.99750155);
415        assert_eq!(f_j0f(15.1), -0.03456193);
416        assert_eq!(f_j0f(f32::INFINITY), 0.);
417        assert_eq!(f_j0f(f32::NEG_INFINITY), 0.);
418        assert!(f_j0f(f32::NAN).is_nan());
419    }
420}
pxfm/bessel/j0f.rs

pxfm/bessel/
j0f.rs
