pxfm/bessel/
y0f.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::{j0f_asympt_alpha, j0f_asympt_beta, j1f_rsqrt};
30use crate::bessel::trigo_bessel::sin_small;
31use crate::bessel::y0f_coeffs::{Y0_ZEROS, Y0_ZEROS_VALUES, Y0F_COEFFS};
32use crate::common::f_fmla;
33use crate::double_double::DoubleDouble;
34use crate::logs::fast_logf;
35use crate::polyeval::{f_polyeval10, f_polyeval18};
36use crate::rounding::CpuCeil;
37use crate::sincos_reduce::rem2pif_any;
38
39/// Bessel of the second kind of order 0 (Y0)
40///
41/// Max ULP 0.5
42pub fn f_y0f(x: f32) -> f32 {
43    let ux = x.to_bits();
44    if ux >= 0xffu32 << 23 || ux == 0 {
45        // |x| == 0, |x| == inf, |x| == NaN, x < 0
46        if ux.wrapping_shl(1) == 0 {
47            return f32::NEG_INFINITY;
48        }
49
50        if x.is_infinite() {
51            if x.is_sign_negative() {
52                return f32::NAN;
53            }
54            return 0.;
55        }
56
57        return x + f32::NAN; // x == NaN
58    }
59
60    let xb = x.to_bits();
61
62    if xb <= 0x4296999au32 {
63        // x <= 75.3
64        if xb <= 0x40000000u32 {
65            // x <= 2
66            if xb <= 0x3faccccdu32 {
67                // x <= 1.35
68                return y0f_near_zero(f32::from_bits(xb));
69            }
70            // transient zone from 1.35 to 2 have bad behavior for log poly already,
71            // and not yet good to be easily covered, thus it use its own poly
72            return y0_transient_area(x);
73        }
74
75        return y0f_small_argument_path(f32::from_bits(xb));
76    }
77
78    // Exceptions:
79    let xb = x.to_bits();
80    if xb == 0x5023e87f {
81        return f32::from_bits(0x28085b2d);
82    } else if xb == 0x48171521 {
83        return f32::from_bits(0x2bd244ba);
84    } else if xb == 0x4398c299 {
85        return f32::from_bits(0x32c730db);
86    } else if xb == 0x7f0e5a38 {
87        return f32::from_bits(0x131f680b);
88    } else if xb == 0x6ef9be45 {
89        return f32::from_bits(0x987d8a8f);
90    }
91
92    y0f_asympt(x)
93}
94
95/**
96Generated by SageMath:
97Evaluates:
98Y0(x) = 2/pi*(euler_gamma + log(x/2))*J0(x) - sum((-1)^m*(x/2)^(2*m)/(m!)^2*sum(1+1/2 + ... 1/m))
99expressed as:
100Y0(x)=log(x)*W0(x) - Z0(x)
101```python
102from sage.all import *
103
104R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
105x = R.gen()
106N = 10  # Number of terms (adjust as needed)
107gamma = RealField(300)(euler_gamma)
108d2 = RealField(300)(2)
109pi = RealField(300).pi()
110
111# Define J0(x) Taylor expansion at x = 0
112def j_series(n, x):
113    return sum([(-1)**m * (x/2)**(ZZ(n) + ZZ(2)*ZZ(m)) / (ZZ(m).factorial() * (ZZ(m) + ZZ(n)).factorial()) for m in range(N)])
114
115J0_series = j_series(0, x)
116
117def z_series(x):
118    return sum([(-1)**m * (x/2)**(ZZ(2)*ZZ(m)) / ZZ(m).factorial()**ZZ(2) * sum(RealField(300)(1)/RealField(300)(k) for k in range(1, m+1)) for m in range(1, N)])
119
120W0 = (d2/pi) * J0_series
121Z0 = -gamma * (d2/pi) * J0_series + RealField(300)(2).log() * (d2/pi) * J0_series + (d2/pi) * z_series(x)
122
123# see the series
124print(W0)
125print(Z0)
126```
127**/
128#[inline]
129fn y0f_near_zero(x: f32) -> f32 {
130    const W: [u64; 10] = [
131        0x3fe45f306dc9c883,
132        0xbfc45f306dc9c883,
133        0x3f845f306dc9c883,
134        0xbf321bb945252402,
135        0x3ed21bb945252402,
136        0xbe672db9f21b0f5f,
137        0x3df49a6c656d62ff,
138        0xbd7ae90af76a4d0f,
139        0x3cfae90af76a4d0f,
140        0xbc754331c053fdad,
141    ];
142    let dx = x as f64;
143    let x2 = dx * dx;
144    let w0 = f_polyeval10(
145        x2,
146        f64::from_bits(W[0]),
147        f64::from_bits(W[1]),
148        f64::from_bits(W[2]),
149        f64::from_bits(W[3]),
150        f64::from_bits(W[4]),
151        f64::from_bits(W[5]),
152        f64::from_bits(W[6]),
153        f64::from_bits(W[7]),
154        f64::from_bits(W[8]),
155        f64::from_bits(W[9]),
156    );
157    const Z: [u64; 10] = [
158        0x3fb2e4d699cbd01f,
159        0xbfc6bbcb41034286,
160        0x3f9075b1bbf41364,
161        0xbf41a6206b7b973d,
162        0x3ee3e99794203bbd,
163        0xbe7bce4a600d3ea4,
164        0x3e0a6ee796b871b6,
165        0xbd92393d82c6b2e4,
166        0x3d131085da82054c,
167        0xbc8f4ed4b492ebcc,
168    ];
169    let z0 = f_polyeval10(
170        x2,
171        f64::from_bits(Z[0]),
172        f64::from_bits(Z[1]),
173        f64::from_bits(Z[2]),
174        f64::from_bits(Z[3]),
175        f64::from_bits(Z[4]),
176        f64::from_bits(Z[5]),
177        f64::from_bits(Z[6]),
178        f64::from_bits(Z[7]),
179        f64::from_bits(Z[8]),
180        f64::from_bits(Z[9]),
181    );
182    let w_log = fast_logf(x);
183    f_fmla(w0, w_log, -z0) as f32
184}
185
186#[inline]
187fn y0_transient_area(x: f32) -> f32 {
188    let dx = x as f64;
189    // first Y0 bessel zero
190    const ZERO: DoubleDouble =
191        DoubleDouble::from_bit_pair((0xbc8bd1e50d219bfd, 0x400193bed4dff243));
192    let r = (dx - ZERO.hi) - ZERO.lo;
193    /*
194        Poly generated by Wolfram Matematica:
195        <<FunctionApproximations`
196        ClearAll["Global`*"]
197        f[x_]:= BesselY[0,x + 2.1971413260310170351490335626990]
198        {approx,error} = MiniMaxApproximation[f[x],{x,{ 1.35 - 2.1971413260310170351490335626990, 2 - 2.1971413260310170351490335626990 },17,0},WorkingPrecision->120]
199        poly=error[[1]];
200        coeffs=CoefficientList[poly,x];
201        TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
202    */
203    let p = f_polyeval18(
204        r,
205        f64::from_bits(0x3fe0aa48442f8375),
206        f64::from_bits(0x3de601d3b959b8d8),
207        f64::from_bits(0xbfd0aa4840bb8529),
208        f64::from_bits(0x3fa439fc16d4835e),
209        f64::from_bits(0x3f80d2dcd97d2b4f),
210        f64::from_bits(0x3f4f833368f9f047),
211        f64::from_bits(0xbf541a702ee92277),
212        f64::from_bits(0x3f3abc113cf0f4da),
213        f64::from_bits(0xbefac1ded6f17ba8),
214        f64::from_bits(0x3f33ef372e24df82),
215        f64::from_bits(0x3f3bf8b42322df40),
216        f64::from_bits(0x3f4582f9daec9ca7),
217        f64::from_bits(0x3f479fc07175494e),
218        f64::from_bits(0x3f4477a5e32b723a),
219        f64::from_bits(0x3f39fbfd6a6d6f0c),
220        f64::from_bits(0x3f2760a66816527b),
221        f64::from_bits(0x3f0a68fdeeba224f),
222        f64::from_bits(0x3edd78c6c87089e1),
223    );
224    p as f32
225}
226
227/// This method on small range searches for nearest zero or extremum.
228/// Then picks stored series expansion at the point end evaluates the poly at the point.
229#[inline]
230fn y0f_small_argument_path(x: f32) -> f32 {
231    let x_abs = x as f64;
232
233    // let avg_step = 74.607799 / 47.0;
234    // let inv_step = 1.0 / avg_step;
235
236    const INV_STEP: f64 = 0.6299609508652038;
237
238    let fx = x_abs * INV_STEP;
239    const Y0_ZEROS_COUNT: f64 = (Y0_ZEROS.len() - 1) as f64;
240    let idx0 = unsafe { fx.min(Y0_ZEROS_COUNT).to_int_unchecked::<usize>() };
241    let idx1 = unsafe {
242        fx.cpu_ceil()
243            .min(Y0_ZEROS_COUNT)
244            .to_int_unchecked::<usize>()
245    };
246
247    let found_zero0 = DoubleDouble::from_bit_pair(Y0_ZEROS[idx0]);
248    let found_zero1 = DoubleDouble::from_bit_pair(Y0_ZEROS[idx1]);
249
250    let dist0 = (found_zero0.hi - x_abs).abs();
251    let dist1 = (found_zero1.hi - x_abs).abs();
252
253    let (found_zero, idx, dist) = if dist0 < dist1 {
254        (found_zero0, idx0, dist0)
255    } else {
256        (found_zero1, idx1, dist1)
257    };
258
259    if idx == 0 {
260        // Really should not happen here, but if it is then to log expansion
261        return y0f_near_zero(x);
262    }
263
264    // We hit exact zero, value, better to return it directly
265    if dist == 0. {
266        return f64::from_bits(Y0_ZEROS_VALUES[idx]) as f32;
267    }
268
269    let c = &Y0F_COEFFS[idx - 1];
270
271    let r = (x_abs - found_zero.hi) - found_zero.lo;
272
273    let p = f_polyeval18(
274        r,
275        f64::from_bits(c[0]),
276        f64::from_bits(c[1]),
277        f64::from_bits(c[2]),
278        f64::from_bits(c[3]),
279        f64::from_bits(c[4]),
280        f64::from_bits(c[5]),
281        f64::from_bits(c[6]),
282        f64::from_bits(c[7]),
283        f64::from_bits(c[8]),
284        f64::from_bits(c[9]),
285        f64::from_bits(c[10]),
286        f64::from_bits(c[11]),
287        f64::from_bits(c[12]),
288        f64::from_bits(c[13]),
289        f64::from_bits(c[14]),
290        f64::from_bits(c[15]),
291        f64::from_bits(c[16]),
292        f64::from_bits(c[17]),
293    );
294
295    p as f32
296}
297
298/*
299   Evaluates:
300   Y0 = sqrt(2/(PI*x)) * beta(x) * sin(x - PI/4 - alpha(x))
301*/
302#[inline]
303fn y0f_asympt(x: f32) -> f32 {
304    let dx = x as f64;
305
306    let alpha = j0f_asympt_alpha(dx);
307    let beta = j0f_asympt_beta(dx);
308
309    let angle = rem2pif_any(x);
310
311    const SQRT_2_OVER_PI: f64 = f64::from_bits(0x3fe9884533d43651);
312    const MPI_OVER_4: f64 = f64::from_bits(0xbfe921fb54442d18);
313
314    let x0pi34 = MPI_OVER_4 - alpha;
315    let r0 = angle + x0pi34;
316
317    let m_cos = sin_small(r0);
318
319    let z0 = beta * m_cos;
320    let scale = SQRT_2_OVER_PI * j1f_rsqrt(dx);
321
322    (scale * z0) as f32
323}
324
325#[cfg(test)]
326mod tests {
327    use crate::f_y0f;
328
329    #[test]
330    fn test_y0f() {
331        assert_eq!(f_y0f(90.5), 0.08254846);
332        assert_eq!(f_y0f(77.5), 0.087678276);
333        assert_eq!(f_y0f(1.5), 0.3824489);
334        assert_eq!(f_y0f(0.5), -0.44451874);
335        assert!(f_y0f(-1.).is_nan());
336        assert_eq!(f_y0f(0.), f32::NEG_INFINITY);
337        assert_eq!(f_y0f(-0.), f32::NEG_INFINITY);
338        assert_eq!(f_y0f(f32::INFINITY), 0.);
339        assert!(f_y0f(f32::NEG_INFINITY).is_nan());
340    }
341}
pxfm/bessel/y0f.rs

pxfm/bessel/
y0f.rs
