pxfm/bessel/
jincpi.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 */
29
30#![allow(clippy::excessive_precision)]
31
32use crate::bessel::alpha1::bessel_1_asympt_alpha_fast;
33use crate::bessel::beta1::bessel_1_asympt_beta_fast;
34use crate::bessel::j1_coeffs::{J1_COEFFS, J1_ZEROS, J1_ZEROS_VALUE};
35use crate::bessel::j1_coeffs_taylor::J1_COEFFS_TAYLOR;
36use crate::common::f_fmla;
37use crate::double_double::DoubleDouble;
38use crate::polyeval::{f_polyeval9, f_polyeval19};
39use crate::round::RoundFinite;
40use crate::sin_helper::sin_dd_small_fast;
41
42/// Normalized jinc 2*J1(PI\*x)/(pi\*x)
43pub fn f_jincpi(x: f64) -> f64 {
44    let ux = x.to_bits().wrapping_shl(1);
45
46    if ux >= 0x7ffu64 << 53 || ux <= 0x7960000000000000u64 {
47        // |x| <= f64::EPSILON, |x| == inf, x == NaN
48        if ux <= 0x7960000000000000u64 {
49            // |x| <= f64::EPSILON
50            return 1.0;
51        }
52        if x.is_infinite() {
53            return 0.;
54        }
55        return x + f64::NAN; // x = NaN
56    }
57
58    let ax: u64 = x.to_bits() & 0x7fff_ffff_ffff_ffff;
59
60    if ax < 0x4052a6784230fcf8u64 {
61        // |x| < 74.60109
62        if ax < 0x3fd3333333333333 {
63            // |x| < 0.3
64            return jincpi_near_zero(f64::from_bits(ax));
65        }
66        let scaled_pix = f64::from_bits(ax) * std::f64::consts::PI; // just test boundaries
67        if scaled_pix < 74.60109 {
68            return jinc_small_argument_fast(f64::from_bits(ax));
69        }
70    }
71
72    jinc_asympt_fast(f64::from_bits(ax))
73}
74
75/*
76   Evaluates:
77   J1 = sqrt(2/(PI*x)) * beta(x) * cos(x - 3*PI/4 - alpha(x))
78   discarding 1*PI/2 using identities gives:
79   J1 = sqrt(2/(PI*x)) * beta(x) * sin(x - PI/4 - alpha(x))
80
81   to avoid squashing small (-PI/4 - alpha(x)) into a large x actual expansion is:
82
83   J1 = sqrt(2/(PI*x)) * beta(x) * sin((x mod 2*PI) - PI/4 - alpha(x))
84*/
85#[inline]
86fn jinc_asympt_fast(ox: f64) -> f64 {
87    const PI: DoubleDouble = DoubleDouble::new(
88        f64::from_bits(0x3ca1a62633145c07),
89        f64::from_bits(0x400921fb54442d18),
90    );
91
92    let x = DoubleDouble::quick_mult_f64(PI, ox);
93
94    const SQRT_2_OVER_PI: DoubleDouble = DoubleDouble::new(
95        f64::from_bits(0xbc8cbc0d30ebfd15),
96        f64::from_bits(0x3fe9884533d43651),
97    );
98    const MPI_OVER_4: DoubleDouble = DoubleDouble::new(
99        f64::from_bits(0xbc81a62633145c07),
100        f64::from_bits(0xbfe921fb54442d18),
101    );
102
103    // argument reduction assuming x here value is already multiple of PI.
104    // k = round((x*Pi) / (pi*2))
105    let kd = (ox * 0.5).round_finite();
106    //  y = (x * Pi) - k * 2
107    let rem = f_fmla(kd, -2., ox);
108    let angle = DoubleDouble::quick_mult_f64(PI, rem);
109
110    let recip = x.recip();
111
112    let alpha = bessel_1_asympt_alpha_fast(recip);
113    let beta = bessel_1_asympt_beta_fast(recip);
114
115    // Without full subtraction cancellation happens sometimes
116    let x0pi34 = DoubleDouble::full_dd_sub(MPI_OVER_4, alpha);
117    let r0 = DoubleDouble::full_dd_add(angle, x0pi34);
118
119    let m_sin = sin_dd_small_fast(r0);
120    let z0 = DoubleDouble::quick_mult(beta, m_sin);
121    let dx_sqrt = x.fast_sqrt();
122    let scale = DoubleDouble::div(SQRT_2_OVER_PI, dx_sqrt);
123    let p = DoubleDouble::quick_mult(scale, z0);
124
125    DoubleDouble::quick_mult(p, recip).to_f64() * 2.
126}
127
128#[inline]
129pub(crate) fn jincpi_near_zero(x: f64) -> f64 {
130    // Polynomial Generated by Wolfram Mathematica:
131    // <<FunctionApproximations`
132    // ClearAll["Global`*"]
133    // f[x_]:=BesselJ[1,x*Pi]/(x*Pi)
134    // {err,approx}=MiniMaxApproximation[f[z],{z,{2^-23,0.3},7,7},WorkingPrecision->60]
135    // poly=Numerator[approx][[1]];
136    // coeffs=CoefficientList[poly,z];
137    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
138    // poly=Denominator[approx][[1]];
139    // coeffs=CoefficientList[poly,z];
140    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
141    const P: [(u64, u64); 8] = [
142        (0xbb2bddffe9450ca6, 0x3fe0000000000000),
143        (0x3c3b0b0a7393eccb, 0xbfce4cd3c3c87615),
144        (0xbc7f9f784e0594a6, 0xbfe043283b1e383f),
145        (0xbc6af77bca466875, 0x3fcee46673cf919f),
146        (0xbc0b62837b038ea8, 0x3fc0b7cc55c9a4af),
147        (0x3c5c08841871f124, 0xbfb002b65231dcdd),
148        (0xbc26cf2d89ea63bc, 0xbf849022a7a0712b),
149        (0xbbe535d492c0ac1c, 0x3f740b48910d5105),
150    ];
151
152    const Q: [(u64, u64); 8] = [
153        (0x0000000000000000, 0x3ff0000000000000),
154        (0x3c4aba6577f3253e, 0xbfde4cd3c3c87615),
155        (0x3c52f58f82e3438c, 0x3fcbd0a475006cf9),
156        (0x3c36e496237d6b49, 0xbfb9f4cea13b06e9),
157        (0xbbbbf3e4ef3a28fe, 0x3f967ed0cee85392),
158        (0x3c267ac442bb3bcf, 0xbf846e192e22f862),
159        (0x3bd84e9888993cb0, 0x3f51e0fff3cfddee),
160        (0x3bd7c0285797bd8e, 0xbf3ea7a621fa1c8c),
161    ];
162
163    let x2 = DoubleDouble::from_exact_mult(x, x);
164    let x4 = x2 * x2;
165
166    let p0 = DoubleDouble::mul_f64_add(
167        DoubleDouble::from_bit_pair(P[1]),
168        x,
169        DoubleDouble::from_bit_pair(P[0]),
170    );
171    let p1 = DoubleDouble::mul_f64_add(
172        DoubleDouble::from_bit_pair(P[3]),
173        x,
174        DoubleDouble::from_bit_pair(P[2]),
175    );
176    let p2 = DoubleDouble::mul_f64_add(
177        DoubleDouble::from_bit_pair(P[5]),
178        x,
179        DoubleDouble::from_bit_pair(P[4]),
180    );
181    let p3 = DoubleDouble::mul_f64_add(
182        DoubleDouble::from_bit_pair(P[7]),
183        x,
184        DoubleDouble::from_bit_pair(P[6]),
185    );
186
187    let q0 = DoubleDouble::mul_add(x2, p1, p0);
188    let q1 = DoubleDouble::mul_add(x2, p3, p2);
189
190    let p_num = DoubleDouble::mul_add(x4, q1, q0);
191
192    let p0 = DoubleDouble::mul_f64_add(
193        DoubleDouble::from_bit_pair(Q[1]),
194        x,
195        DoubleDouble::from_bit_pair(Q[0]),
196    );
197    let p1 = DoubleDouble::mul_f64_add(
198        DoubleDouble::from_bit_pair(Q[3]),
199        x,
200        DoubleDouble::from_bit_pair(Q[2]),
201    );
202    let p2 = DoubleDouble::mul_f64_add(
203        DoubleDouble::from_bit_pair(Q[5]),
204        x,
205        DoubleDouble::from_bit_pair(Q[4]),
206    );
207    let p3 = DoubleDouble::mul_f64_add(
208        DoubleDouble::from_bit_pair(Q[7]),
209        x,
210        DoubleDouble::from_bit_pair(Q[6]),
211    );
212
213    let q0 = DoubleDouble::mul_add(x2, p1, p0);
214    let q1 = DoubleDouble::mul_add(x2, p3, p2);
215
216    let p_den = DoubleDouble::mul_add(x4, q1, q0);
217
218    DoubleDouble::quick_mult_f64(DoubleDouble::div(p_num, p_den), 2.).to_f64()
219}
220
221/// This method on small range searches for nearest zero or extremum.
222/// Then picks stored series expansion at the point end evaluates the poly at the point.
223#[inline]
224pub(crate) fn jinc_small_argument_fast(x: f64) -> f64 {
225    const PI: DoubleDouble = DoubleDouble::new(
226        f64::from_bits(0x3ca1a62633145c07),
227        f64::from_bits(0x400921fb54442d18),
228    );
229
230    // let avg_step = 74.60109 / 47.0;
231    // let inv_step = 1.0 / avg_step;
232
233    let dx = DoubleDouble::quick_mult_f64(PI, x);
234
235    const INV_STEP: f64 = 0.6300176043004198;
236
237    let fx = dx.hi * INV_STEP;
238    const J1_ZEROS_COUNT: f64 = (J1_ZEROS.len() - 1) as f64;
239    let idx0 = unsafe { fx.min(J1_ZEROS_COUNT).to_int_unchecked::<usize>() };
240    let idx1 = unsafe { fx.ceil().min(J1_ZEROS_COUNT).to_int_unchecked::<usize>() };
241
242    let found_zero0 = DoubleDouble::from_bit_pair(J1_ZEROS[idx0]);
243    let found_zero1 = DoubleDouble::from_bit_pair(J1_ZEROS[idx1]);
244
245    let dist0 = (found_zero0.hi - dx.hi).abs();
246    let dist1 = (found_zero1.hi - dx.hi).abs();
247
248    let (found_zero, idx, dist) = if dist0 < dist1 {
249        (found_zero0, idx0, dist0)
250    } else {
251        (found_zero1, idx1, dist1)
252    };
253
254    if idx == 0 {
255        return jincpi_near_zero(x);
256    }
257
258    let r = DoubleDouble::quick_dd_sub(dx, found_zero);
259
260    // We hit exact zero, value, better to return it directly
261    if dist == 0. {
262        return DoubleDouble::quick_mult_f64(
263            DoubleDouble::from_f64_div_dd(f64::from_bits(J1_ZEROS_VALUE[idx]), dx),
264            2.,
265        )
266        .to_f64();
267    }
268
269    let is_zero_too_close = dist.abs() < 1e-3;
270
271    let c = if is_zero_too_close {
272        &J1_COEFFS_TAYLOR[idx - 1]
273    } else {
274        &J1_COEFFS[idx - 1]
275    };
276
277    let p = f_polyeval19(
278        r.hi,
279        f64::from_bits(c[5].1),
280        f64::from_bits(c[6].1),
281        f64::from_bits(c[7].1),
282        f64::from_bits(c[8].1),
283        f64::from_bits(c[9].1),
284        f64::from_bits(c[10].1),
285        f64::from_bits(c[11].1),
286        f64::from_bits(c[12].1),
287        f64::from_bits(c[13].1),
288        f64::from_bits(c[14].1),
289        f64::from_bits(c[15].1),
290        f64::from_bits(c[16].1),
291        f64::from_bits(c[17].1),
292        f64::from_bits(c[18].1),
293        f64::from_bits(c[19].1),
294        f64::from_bits(c[20].1),
295        f64::from_bits(c[21].1),
296        f64::from_bits(c[22].1),
297        f64::from_bits(c[23].1),
298    );
299
300    let mut z = DoubleDouble::mul_f64_add(r, p, DoubleDouble::from_bit_pair(c[4]));
301    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[3]));
302    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[2]));
303    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[1]));
304    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[0]));
305
306    z = DoubleDouble::quick_mult_f64(DoubleDouble::div(z, dx), 2.);
307
308    let err = f_fmla(
309        z.hi,
310        f64::from_bits(0x3c70000000000000), // 2^-56
311        f64::from_bits(0x3bf0000000000000), // 2^-64
312    );
313    let ub = z.hi + (z.lo + err);
314    let lb = z.hi + (z.lo - err);
315
316    if ub == lb {
317        return z.to_f64();
318    }
319
320    j1_small_argument_dd(r, c, dx)
321}
322
323fn j1_small_argument_dd(r: DoubleDouble, c0: &[(u64, u64); 24], inv_scale: DoubleDouble) -> f64 {
324    let c = &c0[15..];
325
326    let p0 = f_polyeval9(
327        r.to_f64(),
328        f64::from_bits(c[0].1),
329        f64::from_bits(c[1].1),
330        f64::from_bits(c[2].1),
331        f64::from_bits(c[3].1),
332        f64::from_bits(c[4].1),
333        f64::from_bits(c[5].1),
334        f64::from_bits(c[6].1),
335        f64::from_bits(c[7].1),
336        f64::from_bits(c[8].1),
337    );
338
339    let c = c0;
340
341    let mut p_e = DoubleDouble::mul_f64_add(r, p0, DoubleDouble::from_bit_pair(c[14]));
342    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[13]));
343    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[12]));
344    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[11]));
345    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[10]));
346    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[9]));
347    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[8]));
348    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[7]));
349    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[6]));
350    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[5]));
351    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[4]));
352    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[3]));
353    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[2]));
354    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[1]));
355    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[0]));
356
357    let p = DoubleDouble::from_exact_add(p_e.hi, p_e.lo);
358    let z = DoubleDouble::div(p, inv_scale);
359    DoubleDouble::quick_mult_f64(z, 2.).to_f64()
360}
361
362#[cfg(test)]
363mod tests {
364    use super::*;
365
366    #[test]
367    fn test_jincpi() {
368        assert_eq!(f_jincpi(f64::EPSILON), 1.0);
369        assert_eq!(f_jincpi(0.5000000000020244), 0.7217028449014163);
370        assert_eq!(f_jincpi(73.81695991658546), -0.0004417546638317049);
371        assert_eq!(f_jincpi(0.01), 0.9998766350182722);
372        assert_eq!(f_jincpi(0.9), 0.28331697846510623);
373        assert_eq!(f_jincpi(3.831705970207517), -0.036684415010255086);
374        assert_eq!(f_jincpi(-3.831705970207517), -0.036684415010255086);
375        assert_eq!(
376            f_jincpi(0.000000000000000000000000000000000000008827127),
377            1.0
378        );
379        assert_eq!(
380            f_jincpi(-0.000000000000000000000000000000000000008827127),
381            1.0
382        );
383        assert_eq!(f_jincpi(5.4), -0.010821736808448256);
384        assert_eq!(
385            f_jincpi(77.743162408196766932633181568235159),
386            -0.00041799098646950523
387        );
388        assert_eq!(
389            f_jincpi(84.027189586293545175976760219782591),
390            -0.00023927934929850555
391        );
392        assert_eq!(f_jincpi(f64::NEG_INFINITY), 0.0);
393        assert_eq!(f_jincpi(f64::INFINITY), 0.0);
394        assert!(f_jincpi(f64::NAN).is_nan());
395    }
396}
pxfm/bessel/jincpi.rs

pxfm/bessel/
jincpi.rs
