pxfm/bessel/
j0.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::alpha0::{
30    bessel_0_asympt_alpha, bessel_0_asympt_alpha_fast, bessel_0_asympt_alpha_hard,
31};
32use crate::bessel::beta0::{
33    bessel_0_asympt_beta, bessel_0_asympt_beta_fast, bessel_0_asympt_beta_hard,
34};
35use crate::bessel::i0::bessel_rsqrt_hard;
36use crate::bessel::j0_coeffs_remez::J0_COEFFS_REMEZ;
37use crate::bessel::j0_coeffs_taylor::J0_COEFFS_TAYLOR;
38use crate::bessel::j0f_coeffs::{J0_ZEROS, J0_ZEROS_VALUE};
39use crate::common::f_fmla;
40use crate::double_double::DoubleDouble;
41use crate::dyadic_float::{DyadicFloat128, DyadicSign};
42use crate::polyeval::{f_polyeval9, f_polyeval10, f_polyeval12, f_polyeval19};
43use crate::sin_helper::{cos_dd_small, cos_dd_small_fast, cos_f128_small};
44use crate::sincos_reduce::{AngleReduced, rem2pi_any, rem2pi_f128};
45
46/// Bessel of the first kind of order 0
47pub fn f_j0(x: f64) -> f64 {
48    let ux = x.to_bits().wrapping_shl(1);
49
50    if ux >= 0x7ffu64 << 53 || ux <= 0x7960000000000000u64 {
51        // |x| <= f64::EPSILON, |x| == inf, |x| == NaN
52        if ux <= 0x77723ef88126da90u64 {
53            // |x| <= 0.00000000000000000000532
54            return 1.;
55        }
56        if ux <= 0x7960000000000000u64 {
57            // |x| <= f64::EPSILON
58            // J0(x) ~ 1-x^2/4+O[x]^4
59            let half_x = 0.5 * x; // exact.
60            return f_fmla(half_x, -half_x, 1.);
61        }
62        if x.is_infinite() {
63            return 0.;
64        }
65        return x + f64::NAN; // x == NaN
66    }
67
68    let x_abs = x.to_bits() & 0x7fff_ffff_ffff_ffff;
69
70    let ax = f64::from_bits(x_abs);
71
72    if x_abs <= 0x4052b33333333333u64 {
73        // |x| <= 74.8
74        if x_abs <= 0x3ff199999999999au64 {
75            // |x| <= 1.1
76            return j0_maclaurin_series_fast(ax);
77        }
78        return j0_small_argument_fast(ax);
79    }
80
81    j0_asympt_fast(ax)
82}
83
84/**
85Generated by SageMath:
86```python
87mp.prec = 180
88def print_expansion_at_0():
89    print(f"const J0_MACLAURIN_SERIES: [(u64, u64); 12] = [")
90    from mpmath import mp, j0, taylor
91    poly = taylor(lambda val: j0(val), 0, 24)
92    real_i = 0
93    for i in range(0, 24, 2):
94        print_double_double("", DD(poly[i]))
95        real_i = real_i + 1
96    print("];")
97    print(poly)
98
99print_expansion_at_0()
100```
101**/
102#[inline]
103pub(crate) fn j0_maclaurin_series_fast(x: f64) -> f64 {
104    const C: [u64; 12] = [
105        0x3ff0000000000000,
106        0xbfd0000000000000,
107        0x3f90000000000000,
108        0xbf3c71c71c71c71c,
109        0x3edc71c71c71c71c,
110        0xbe723456789abcdf,
111        0x3e002e85c0898b71,
112        0xbd8522a43f65486a,
113        0x3d0522a43f65486a,
114        0xbc80b313289be0b9,
115        0x3bf5601885e63e5d,
116        0xbb669ca9cf3b7f54,
117    ];
118
119    let dx2 = DoubleDouble::from_exact_mult(x, x);
120
121    let p = f_polyeval10(
122        dx2.hi,
123        f64::from_bits(C[2]),
124        f64::from_bits(C[3]),
125        f64::from_bits(C[4]),
126        f64::from_bits(C[5]),
127        f64::from_bits(C[6]),
128        f64::from_bits(C[7]),
129        f64::from_bits(C[8]),
130        f64::from_bits(C[9]),
131        f64::from_bits(C[10]),
132        f64::from_bits(C[11]),
133    );
134    let mut z = DoubleDouble::mul_f64_add_f64(dx2, p, f64::from_bits(C[1]));
135    z = DoubleDouble::mul_add_f64(dx2, z, f64::from_bits(C[0]));
136
137    // squaring error (2^-56) + poly error 2^-75
138    let err = f_fmla(
139        dx2.hi,
140        f64::from_bits(0x3c70000000000000), // 2^-56
141        f64::from_bits(0x3b40000000000000), // 2^-75
142    );
143    let ub = z.hi + (z.lo + err);
144    let lb = z.hi + (z.lo - err);
145
146    if ub == lb {
147        return z.to_f64();
148    }
149    j0_maclaurin_series(x)
150}
151
152/**
153Generated by SageMath:
154```python
155mp.prec = 180
156def print_expansion_at_0():
157    print(f"const J0_MACLAURIN_SERIES: [(u64, u64); 12] = [")
158    from mpmath import mp, j0, taylor
159    poly = taylor(lambda val: j0(val), 0, 24)
160    real_i = 0
161    for i in range(0, 24, 2):
162        print_double_double("", DD(poly[i]))
163        real_i = real_i + 1
164    print("];")
165    print(poly)
166
167print_expansion_at_0()
168```
169**/
170#[cold]
171pub(crate) fn j0_maclaurin_series(x: f64) -> f64 {
172    const C: [(u64, u64); 12] = [
173        (0x0000000000000000, 0x3ff0000000000000),
174        (0x0000000000000000, 0xbfd0000000000000),
175        (0x0000000000000000, 0x3f90000000000000),
176        (0xbbdc71c71c71c71c, 0xbf3c71c71c71c71c),
177        (0x3b7c71c71c71c71c, 0x3edc71c71c71c71c),
178        (0xbab23456789abcdf, 0xbe723456789abcdf),
179        (0xba8b6edec0692e65, 0x3e002e85c0898b71),
180        (0x3a2604db055bd075, 0xbd8522a43f65486a),
181        (0xb9a604db055bd075, 0x3d0522a43f65486a),
182        (0x3928824198c6f6e1, 0xbc80b313289be0b9),
183        (0xb869b0b430eb27b8, 0x3bf5601885e63e5d),
184        (0x380ee6b4638f3a25, 0xbb669ca9cf3b7f54),
185    ];
186
187    let dx2 = DoubleDouble::from_exact_mult(x, x);
188
189    let p = f_polyeval12(
190        dx2,
191        DoubleDouble::from_bit_pair(C[0]),
192        DoubleDouble::from_bit_pair(C[1]),
193        DoubleDouble::from_bit_pair(C[2]),
194        DoubleDouble::from_bit_pair(C[3]),
195        DoubleDouble::from_bit_pair(C[4]),
196        DoubleDouble::from_bit_pair(C[5]),
197        DoubleDouble::from_bit_pair(C[6]),
198        DoubleDouble::from_bit_pair(C[7]),
199        DoubleDouble::from_bit_pair(C[8]),
200        DoubleDouble::from_bit_pair(C[9]),
201        DoubleDouble::from_bit_pair(C[10]),
202        DoubleDouble::from_bit_pair(C[11]),
203    );
204    let r = DoubleDouble::from_exact_add(p.hi, p.lo);
205    const ERR: f64 = f64::from_bits(0x39d0000000000000); // 2^-98
206    let ub = r.hi + (r.lo + ERR);
207    let lb = r.hi + (r.lo - ERR);
208    if ub == lb {
209        return r.to_f64();
210    }
211    j0_maclaurin_series_hard(x)
212}
213
214/**
215Generated by SageMath:
216
217```python
218mp.prec = 180
219def print_expansion_at_0():
220    print(f"const P: [DyadicFloat128; 12] = [")
221    from mpmath import mp, j0, taylor
222    poly = taylor(lambda val: j0(val), 0, 24)
223    # print(poly)
224    real_i = 0
225    for i in range(0, 24, 2):
226        print_dyadic(DD(poly[i]))
227        real_i = real_i + 1
228    print("];")
229    print(poly)
230
231print_expansion_at_0()
232```
233**/
234#[cold]
235#[inline(never)]
236pub(crate) fn j0_maclaurin_series_hard(x: f64) -> f64 {
237    static P: [DyadicFloat128; 12] = [
238        DyadicFloat128 {
239            sign: DyadicSign::Pos,
240            exponent: -127,
241            mantissa: 0x80000000_00000000_00000000_00000000_u128,
242        },
243        DyadicFloat128 {
244            sign: DyadicSign::Neg,
245            exponent: -129,
246            mantissa: 0x80000000_00000000_00000000_00000000_u128,
247        },
248        DyadicFloat128 {
249            sign: DyadicSign::Pos,
250            exponent: -133,
251            mantissa: 0x80000000_00000000_00000000_00000000_u128,
252        },
253        DyadicFloat128 {
254            sign: DyadicSign::Neg,
255            exponent: -139,
256            mantissa: 0xe38e38e3_8e38e38e_38e38e38_e38e38e4_u128,
257        },
258        DyadicFloat128 {
259            sign: DyadicSign::Pos,
260            exponent: -145,
261            mantissa: 0xe38e38e3_8e38e38e_38e38e38_e38e38e4_u128,
262        },
263        DyadicFloat128 {
264            sign: DyadicSign::Neg,
265            exponent: -151,
266            mantissa: 0x91a2b3c4_d5e6f809_1a2b3c4d_5e6f8092_u128,
267        },
268        DyadicFloat128 {
269            sign: DyadicSign::Pos,
270            exponent: -158,
271            mantissa: 0x81742e04_4c5b8724_8909fcb6_8cd4e410_u128,
272        },
273        DyadicFloat128 {
274            sign: DyadicSign::Neg,
275            exponent: -166,
276            mantissa: 0xa91521fb_2a434d3f_649f5485_f169a743_u128,
277        },
278        DyadicFloat128 {
279            sign: DyadicSign::Pos,
280            exponent: -174,
281            mantissa: 0xa91521fb_2a434d3f_649f5485_f169a743_u128,
282        },
283        DyadicFloat128 {
284            sign: DyadicSign::Neg,
285            exponent: -182,
286            mantissa: 0x85989944_df05c4ef_b7cce721_23e1b391_u128,
287        },
288        DyadicFloat128 {
289            sign: DyadicSign::Pos,
290            exponent: -191,
291            mantissa: 0xab00c42f_31f2e799_3d2f3c53_6120e5d8_u128,
292        },
293        DyadicFloat128 {
294            sign: DyadicSign::Neg,
295            exponent: -200,
296            mantissa: 0xb4e54e79_dbfa9c23_29738e18_bb602809_u128,
297        },
298    ];
299    let dx = DyadicFloat128::new_from_f64(x);
300    let x2 = dx * dx;
301
302    let mut p = P[11];
303    for i in (0..11).rev() {
304        p = x2 * p + P[i];
305    }
306
307    p.fast_as_f64()
308}
309
310/// This method on small range searches for nearest zero or extremum.
311/// Then picks stored series expansion at the point end evaluates the poly at the point.
312#[inline]
313pub(crate) fn j0_small_argument_fast(x: f64) -> f64 {
314    // let avg_step = 74.6145 / 47.0;
315    // let inv_step = 1.0 / avg_step;
316
317    const INV_STEP: f64 = 0.6299043751549631;
318
319    let fx = x * INV_STEP;
320    const J0_ZEROS_COUNT: f64 = (J0_ZEROS.len() - 1) as f64;
321    let idx0 = unsafe { fx.min(J0_ZEROS_COUNT).to_int_unchecked::<usize>() };
322    let idx1 = unsafe { fx.ceil().min(J0_ZEROS_COUNT).to_int_unchecked::<usize>() };
323
324    let found_zero0 = DoubleDouble::from_bit_pair(J0_ZEROS[idx0]);
325    let found_zero1 = DoubleDouble::from_bit_pair(J0_ZEROS[idx1]);
326
327    let dist0 = (found_zero0.hi - x).abs();
328    let dist1 = (found_zero1.hi - x).abs();
329
330    let (found_zero, idx, dist) = if dist0 < dist1 {
331        (found_zero0, idx0, dist0)
332    } else {
333        (found_zero1, idx1, dist1)
334    };
335
336    if idx == 0 {
337        return j0_maclaurin_series_fast(x);
338    }
339
340    let is_too_close_too_zero = dist.abs() < 1e-3;
341
342    let c = if is_too_close_too_zero {
343        &J0_COEFFS_TAYLOR[idx - 1]
344    } else {
345        &J0_COEFFS_REMEZ[idx - 1]
346    };
347
348    let r = DoubleDouble::full_add_f64(-found_zero, x.abs());
349
350    // We hit exact zero, value, better to return it directly
351    if dist == 0. {
352        return f64::from_bits(J0_ZEROS_VALUE[idx]);
353    }
354    let p = f_polyeval19(
355        r.hi,
356        f64::from_bits(c[5].1),
357        f64::from_bits(c[6].1),
358        f64::from_bits(c[7].1),
359        f64::from_bits(c[8].1),
360        f64::from_bits(c[9].1),
361        f64::from_bits(c[10].1),
362        f64::from_bits(c[11].1),
363        f64::from_bits(c[12].1),
364        f64::from_bits(c[13].1),
365        f64::from_bits(c[14].1),
366        f64::from_bits(c[15].1),
367        f64::from_bits(c[16].1),
368        f64::from_bits(c[17].1),
369        f64::from_bits(c[18].1),
370        f64::from_bits(c[19].1),
371        f64::from_bits(c[20].1),
372        f64::from_bits(c[21].1),
373        f64::from_bits(c[22].1),
374        f64::from_bits(c[23].1),
375    );
376
377    let mut z = DoubleDouble::mul_f64_add(r, p, DoubleDouble::from_bit_pair(c[4]));
378    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[3]));
379    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[2]));
380    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[1]));
381    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[0]));
382    let err = f_fmla(
383        z.hi,
384        f64::from_bits(0x3c70000000000000), // 2^-56
385        f64::from_bits(0x3bf0000000000000), // 2^-64
386    );
387    let ub = z.hi + (z.lo + err);
388    let lb = z.hi + (z.lo - err);
389
390    if ub == lb {
391        return z.to_f64();
392    }
393
394    j0_small_argument_dd(r, c)
395}
396
397#[cold]
398fn j0_small_argument_dd(r: DoubleDouble, c0: &[(u64, u64); 24]) -> f64 {
399    let c = &c0[15..];
400
401    let p0 = f_polyeval9(
402        r.to_f64(),
403        f64::from_bits(c[0].1),
404        f64::from_bits(c[1].1),
405        f64::from_bits(c[2].1),
406        f64::from_bits(c[3].1),
407        f64::from_bits(c[4].1),
408        f64::from_bits(c[5].1),
409        f64::from_bits(c[6].1),
410        f64::from_bits(c[7].1),
411        f64::from_bits(c[8].1),
412    );
413
414    let c = c0;
415
416    let mut p_e = DoubleDouble::mul_f64_add(r, p0, DoubleDouble::from_bit_pair(c[14]));
417    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[13]));
418    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[12]));
419    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[11]));
420    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[10]));
421    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[9]));
422    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[8]));
423    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[7]));
424    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[6]));
425    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[5]));
426    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[4]));
427    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[3]));
428    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[2]));
429    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[1]));
430    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[0]));
431
432    let p = DoubleDouble::from_exact_add(p_e.hi, p_e.lo);
433    let err = f_fmla(
434        p.hi,
435        f64::from_bits(0x3c10000000000000), // 2^-62
436        f64::from_bits(0x3a90000000000000), // 2^-86
437    );
438    let ub = p.hi + (p.lo + err);
439    let lb = p.hi + (p.lo - err);
440    if ub != lb {
441        return j0_small_argument_hard(r, c);
442    }
443    p.to_f64()
444}
445
446#[cold]
447#[inline(never)]
448fn j0_small_argument_hard(r: DoubleDouble, c: &[(u64, u64); 24]) -> f64 {
449    let mut p = DoubleDouble::from_bit_pair(c[23]);
450    for i in (0..23).rev() {
451        p = DoubleDouble::mul_add(r, p, DoubleDouble::from_bit_pair(c[i]));
452        p = DoubleDouble::from_exact_add(p.hi, p.lo);
453    }
454    p.to_f64()
455}
456
457/*
458   Evaluates:
459   J0 = sqrt(2/(PI*x)) * beta(x) * cos(x - PI/4 - alpha(x))
460*/
461#[inline]
462pub(crate) fn j0_asympt_fast(x: f64) -> f64 {
463    let x = x.abs();
464
465    const SQRT_2_OVER_PI: DoubleDouble = DoubleDouble::new(
466        f64::from_bits(0xbc8cbc0d30ebfd15),
467        f64::from_bits(0x3fe9884533d43651),
468    );
469
470    const MPI_OVER_4: DoubleDouble = DoubleDouble::new(
471        f64::from_bits(0xbc81a62633145c07),
472        f64::from_bits(0xbfe921fb54442d18),
473    );
474
475    let recip = if x.to_bits() > 0x7fd000000000000u64 {
476        DoubleDouble::quick_mult_f64(DoubleDouble::from_exact_safe_div(4.0, x), 0.25)
477    } else {
478        DoubleDouble::from_recip(x)
479    };
480
481    let alpha = bessel_0_asympt_alpha_fast(recip);
482    let beta = bessel_0_asympt_beta_fast(recip);
483
484    let AngleReduced { angle } = rem2pi_any(x);
485
486    // Without full subtraction cancellation happens sometimes
487    let x0pi34 = DoubleDouble::full_dd_sub(MPI_OVER_4, alpha);
488    let r0 = DoubleDouble::full_dd_add(angle, x0pi34);
489
490    let m_cos = cos_dd_small_fast(r0);
491    let z0 = DoubleDouble::quick_mult(beta, m_cos);
492    let r_sqrt = DoubleDouble::from_rsqrt_fast(x);
493    let scale = DoubleDouble::quick_mult(SQRT_2_OVER_PI, r_sqrt);
494    let p = DoubleDouble::quick_mult(scale, z0);
495
496    let err = f_fmla(
497        p.hi,
498        f64::from_bits(0x3be0000000000000), // 2^-65
499        f64::from_bits(0x3a60000000000000), // 2^-89
500    );
501    let ub = p.hi + (p.lo + err);
502    let lb = p.hi + (p.lo - err);
503
504    if ub == lb {
505        return p.to_f64();
506    }
507    j0_asympt(x, recip, r_sqrt, angle)
508}
509
510/*
511   Evaluates:
512   J0 = sqrt(2/(PI*x)) * beta(x) * cos(x - PI/4 - alpha(x))
513*/
514pub(crate) fn j0_asympt(
515    x: f64,
516    recip: DoubleDouble,
517    r_sqrt: DoubleDouble,
518    angle: DoubleDouble,
519) -> f64 {
520    const SQRT_2_OVER_PI: DoubleDouble = DoubleDouble::new(
521        f64::from_bits(0xbc8cbc0d30ebfd15),
522        f64::from_bits(0x3fe9884533d43651),
523    );
524
525    const MPI_OVER_4: DoubleDouble = DoubleDouble::new(
526        f64::from_bits(0xbc81a62633145c07),
527        f64::from_bits(0xbfe921fb54442d18),
528    );
529
530    let alpha = bessel_0_asympt_alpha(recip);
531    let beta = bessel_0_asympt_beta(recip);
532
533    // Without full subtraction cancellation happens sometimes
534    let x0pi34 = DoubleDouble::full_dd_sub(MPI_OVER_4, alpha);
535    let r0 = DoubleDouble::full_dd_add(angle, x0pi34);
536
537    let m_cos = cos_dd_small(r0);
538    let z0 = DoubleDouble::quick_mult(beta, m_cos);
539    let scale = DoubleDouble::quick_mult(SQRT_2_OVER_PI, r_sqrt);
540    let r = DoubleDouble::quick_mult(scale, z0);
541
542    let p = DoubleDouble::from_exact_add(r.hi, r.lo);
543    let err = f_fmla(
544        p.hi,
545        f64::from_bits(0x3bd0000000000000), // 2^-66
546        f64::from_bits(0x39e0000000000000), // 2^-97
547    );
548    let ub = p.hi + (p.lo + err);
549    let lb = p.hi + (p.lo - err);
550
551    if ub == lb {
552        return p.to_f64();
553    }
554    j0_asympt_hard(x)
555}
556
557/*
558   Evaluates:
559   J0 = sqrt(2/(PI*x)) * beta(x) * cos(x - PI/4 - alpha(x))
560*/
561#[cold]
562#[inline(never)]
563pub(crate) fn j0_asympt_hard(x: f64) -> f64 {
564    const SQRT_2_OVER_PI: DyadicFloat128 = DyadicFloat128 {
565        sign: DyadicSign::Pos,
566        exponent: -128,
567        mantissa: 0xcc42299e_a1b28468_7e59e280_5d5c7180_u128,
568    };
569
570    const MPI_OVER_4: DyadicFloat128 = DyadicFloat128 {
571        sign: DyadicSign::Neg,
572        exponent: -128,
573        mantissa: 0xc90fdaa2_2168c234_c4c6628b_80dc1cd1_u128,
574    };
575
576    let x_dyadic = DyadicFloat128::new_from_f64(x);
577    let recip = DyadicFloat128::accurate_reciprocal(x);
578
579    let alpha = bessel_0_asympt_alpha_hard(recip);
580    let beta = bessel_0_asympt_beta_hard(recip);
581
582    let angle = rem2pi_f128(x_dyadic);
583
584    let x0pi34 = MPI_OVER_4 - alpha;
585    let r0 = angle + x0pi34;
586
587    let m_sin = cos_f128_small(r0);
588
589    let z0 = beta * m_sin;
590    let r_sqrt = bessel_rsqrt_hard(x, recip);
591    let scale = SQRT_2_OVER_PI * r_sqrt;
592    let p = scale * z0;
593    p.fast_as_f64()
594}
595
596#[cfg(test)]
597mod tests {
598    use super::*;
599
600    #[test]
601    fn test_j0() {
602        assert_eq!(f_j0(f64::EPSILON), 1.0);
603        assert_eq!(f_j0(-0.000000000000000000000532), 1.0);
604        assert_eq!(f_j0(0.0000000000000000000532), 1.0);
605        assert_eq!(f_j0(-2.000976555054876), 0.22332760641907712);
606        assert_eq!(f_j0(-2.3369499004222215E+304), -3.3630754230844632e-155);
607        assert_eq!(
608            f_j0(f64::from_bits(0xd71a31ffe2ff7e9f)),
609            f64::from_bits(0xb2e58532f95056ff)
610        );
611        assert_eq!(f_j0(6.1795701510782757E+307), 6.075192922402001e-155);
612        assert_eq!(f_j0(6.1795701510782757E+301), 4.118334155030934e-152);
613        assert_eq!(f_j0(6.1795701510782757E+157), 9.5371668900364e-80);
614        assert_eq!(f_j0(79.), -0.08501719554953485);
615        // Without FMA 2.703816901253004e-16
616        #[cfg(any(
617            all(target_arch = "x86_64", target_feature = "fma"),
618            target_arch = "aarch64"
619        ))]
620        assert_eq!(f_j0(93.463718781944774171190), 2.7038169012530046e-16);
621        assert_eq!(f_j0(99.746819858680596470279979), -8.419106281522749e-17);
622        assert_eq!(f_j0(f64::INFINITY), 0.);
623        assert_eq!(f_j0(f64::NEG_INFINITY), 0.);
624        assert!(f_j0(f64::NAN).is_nan());
625    }
626}
pxfm/bessel/j0.rs

pxfm/bessel/
j0.rs
