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

pxfm/bessel/
j0f.rs
