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

pxfm/bessel/
jincpi.rs
