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

pxfm/bessel/
y0f.rs
