pxfm/bessel/
jincpif.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::j1f_rsqrt;
30use crate::bessel::j1_coeffs::{J1_ZEROS, J1_ZEROS_VALUE};
31use crate::bessel::j1f::{j1f_asympt_alpha, j1f_asympt_beta};
32use crate::bessel::j1f_coeffs::J1F_COEFFS;
33use crate::bessel::trigo_bessel::sin_small;
34use crate::common::f_fmla;
35use crate::double_double::DoubleDouble;
36use crate::polyeval::{f_polyeval6, f_polyeval14};
37use crate::rounding::CpuCeil;
38use crate::rounding::CpuRound;
39
40/// Normalized jinc 2*J1(PI\*x)/(pi\*x)
41///
42pub fn f_jincpif(x: f32) -> f32 {
43    let ux = x.to_bits().wrapping_shl(1);
44    if ux >= 0xffu32 << 24 || ux <= 0x6800_0000u32 {
45        // |x| <= f32::EPSILON, |x| == inf, |x| == NaN
46        if ux <= 0x6800_0000u32 {
47            // |x| == 0
48            return 1.;
49        }
50        if x.is_infinite() {
51            return 0.;
52        }
53        return x + f32::NAN; // x == NaN
54    }
55
56    let ax = x.to_bits() & 0x7fff_ffff;
57
58    if ax < 0x429533c2u32 {
59        // |x| < 74.60109
60        if ax <= 0x3e800000u32 {
61            // |x| < 0.25
62            return jincf_near_zero(f32::from_bits(ax));
63        }
64        let scaled_pix = f32::from_bits(ax) * std::f32::consts::PI; // just test boundaries
65        if scaled_pix < 74.60109 {
66            return jincpif_small_argument(f32::from_bits(ax));
67        }
68    }
69
70    jincpif_asympt(f32::from_bits(ax)) as f32
71}
72
73#[inline]
74fn jincf_near_zero(x: f32) -> f32 {
75    let dx = x as f64;
76    // Generated in Wolfram Mathematica:
77    // <<FunctionApproximations`
78    // ClearAll["Global`*"]
79    // f[x_]:=2*BesselJ[1,x*Pi]/(x*Pi)
80    // {err,approx}=MiniMaxApproximation[f[z],{z,{2^-23,0.3},5,5},WorkingPrecision->60]
81    // poly=Numerator[approx][[1]];
82    // coeffs=CoefficientList[poly,z];
83    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
84    // poly=Denominator[approx][[1]];
85    // coeffs=CoefficientList[poly,z];
86    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
87    let p_num = f_polyeval6(
88        dx,
89        f64::from_bits(0x3ff0000000000002),
90        f64::from_bits(0xbfe46cd1822a5aa0),
91        f64::from_bits(0xbfee583c923dc6f4),
92        f64::from_bits(0x3fe3834f47496519),
93        f64::from_bits(0x3fc8118468756e6f),
94        f64::from_bits(0xbfbfaff09f13df88),
95    );
96    let p_den = f_polyeval6(
97        dx,
98        f64::from_bits(0x3ff0000000000000),
99        f64::from_bits(0xbfe46cd1822a4cb0),
100        f64::from_bits(0x3fd2447a026f477a),
101        f64::from_bits(0xbfc6bdf2192404e5),
102        f64::from_bits(0x3fa0cf182218e448),
103        f64::from_bits(0xbf939ab46c3f7a7d),
104    );
105    (p_num / p_den) as f32
106}
107
108/// This method on small range searches for nearest zero or extremum.
109/// Then picks stored series expansion at the point end evaluates the poly at the point.
110#[inline]
111fn jincpif_small_argument(ox: f32) -> f32 {
112    const PI: f64 = f64::from_bits(0x400921fb54442d18);
113    let x = ox as f64 * PI;
114    let x_abs = f64::from_bits(x.to_bits() & 0x7fff_ffff_ffff_ffff);
115
116    // let avg_step = 74.60109 / 47.0;
117    // let inv_step = 1.0 / avg_step;
118
119    const INV_STEP: f64 = 0.6300176043004198;
120
121    let inv_scale = x;
122
123    let fx = x_abs * INV_STEP;
124    const J1_ZEROS_COUNT: f64 = (J1_ZEROS.len() - 1) as f64;
125    let idx0 = unsafe { fx.min(J1_ZEROS_COUNT).to_int_unchecked::<usize>() };
126    let idx1 = unsafe {
127        fx.cpu_ceil()
128            .min(J1_ZEROS_COUNT)
129            .to_int_unchecked::<usize>()
130    };
131
132    let found_zero0 = DoubleDouble::from_bit_pair(J1_ZEROS[idx0]);
133    let found_zero1 = DoubleDouble::from_bit_pair(J1_ZEROS[idx1]);
134
135    let dist0 = (found_zero0.hi - x_abs).abs();
136    let dist1 = (found_zero1.hi - x_abs).abs();
137
138    let (found_zero, idx, dist) = if dist0 < dist1 {
139        (found_zero0, idx0, dist0)
140    } else {
141        (found_zero1, idx1, dist1)
142    };
143
144    if idx == 0 {
145        return jincf_near_zero(ox);
146    }
147
148    // We hit exact zero, value, better to return it directly
149    if dist == 0. {
150        return (f64::from_bits(J1_ZEROS_VALUE[idx]) / inv_scale * 2.) as f32;
151    }
152
153    let c = &J1F_COEFFS[idx - 1];
154
155    let r = (x_abs - found_zero.hi) - found_zero.lo;
156
157    let p = f_polyeval14(
158        r,
159        f64::from_bits(c[0]),
160        f64::from_bits(c[1]),
161        f64::from_bits(c[2]),
162        f64::from_bits(c[3]),
163        f64::from_bits(c[4]),
164        f64::from_bits(c[5]),
165        f64::from_bits(c[6]),
166        f64::from_bits(c[7]),
167        f64::from_bits(c[8]),
168        f64::from_bits(c[9]),
169        f64::from_bits(c[10]),
170        f64::from_bits(c[11]),
171        f64::from_bits(c[12]),
172        f64::from_bits(c[13]),
173    );
174
175    (p / inv_scale * 2.) as f32
176}
177
178/*
179   Evaluates:
180   J1 = sqrt(2/(PI*x)) * beta(x) * cos(x - 3*PI/4 - alpha(x))
181   discarding 1*PI/2 using identities gives:
182   J1 = sqrt(2/(PI*x)) * beta(x) * sin(x - PI/4 - alpha(x))
183
184   to avoid squashing small (-PI/4 - alpha(x)) into a large x actual expansion is:
185
186   J1 = sqrt(2/(PI*x)) * beta(x) * sin((x mod 2*PI) - PI/4 - alpha(x))
187*/
188#[inline]
189pub(crate) fn jincpif_asympt(x: f32) -> f64 {
190    const PI: f64 = f64::from_bits(0x400921fb54442d18);
191
192    let dox = x as f64;
193    let dx = dox * PI;
194
195    let alpha = j1f_asympt_alpha(dx);
196    let beta = j1f_asympt_beta(dx);
197
198    // argument reduction assuming x here value is already multiple of PI.
199    // k = round((x*Pi) / (pi*2))
200    let kd = (dox * 0.5).cpu_round();
201    //  y = (x * Pi) - k * 2
202    let angle = f_fmla(kd, -2., dox) * PI;
203
204    // 2^(3/2)/(Pi^2)
205    // reduce argument 2*sqrt(2/(PI*(x*PI))) = 2*sqrt(2)/PI
206    // adding additional pi from division then 2*sqrt(2)/PI^2
207    const Z2_3_2_OVER_PI_SQR: f64 = f64::from_bits(0x3fd25751e5614413);
208    const MPI_OVER_4: f64 = f64::from_bits(0xbfe921fb54442d18);
209
210    let x0pi34 = MPI_OVER_4 - alpha;
211    let r0 = angle + x0pi34;
212
213    let m_sin = sin_small(r0);
214
215    let z0 = beta * m_sin;
216    let scale = Z2_3_2_OVER_PI_SQR * j1f_rsqrt(dox);
217
218    let j1pix = scale * z0;
219    j1pix / dox
220}
221
222#[cfg(test)]
223mod tests {
224    use super::*;
225
226    #[test]
227    fn test_jincpif() {
228        assert_eq!(f_jincpif(-102.59484), 0.00024380769);
229        assert_eq!(f_jincpif(102.59484), 0.00024380769);
230        assert_eq!(f_jincpif(100.08199), -0.00014386141);
231        assert_eq!(f_jincpif(0.27715185), 0.9081822);
232        assert_eq!(f_jincpif(0.007638072), 0.99992806);
233        assert_eq!(f_jincpif(-f32::EPSILON), 1.0);
234        assert_eq!(f_jincpif(f32::EPSILON), 1.0);
235        assert_eq!(
236            f_jincpif(0.000000000000000000000000000000000000008827127),
237            1.0
238        );
239        assert_eq!(f_jincpif(5.4), -0.010821743);
240        assert_eq!(
241            f_jincpif(77.743162408196766932633181568235159),
242            -0.00041799102
243        );
244        assert_eq!(
245            f_jincpif(-77.743162408196766932633181568235159),
246            -0.00041799102
247        );
248        assert_eq!(
249            f_jincpif(84.027189586293545175976760219782591),
250            -0.00023927793
251        );
252        assert_eq!(f_jincpif(f32::INFINITY), 0.);
253        assert_eq!(f_jincpif(f32::NEG_INFINITY), 0.);
254        assert!(f_jincpif(f32::NAN).is_nan());
255        assert_eq!(f_jincpif(-1.7014118e38), -0.0);
256    }
257}
pxfm/bessel/jincpif.rs

pxfm/bessel/
jincpif.rs
