pxfm/bessel/
y1f.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::j1f::{j1f_asympt_alpha, j1f_asympt_beta};
31use crate::bessel::trigo_bessel::cos_small;
32use crate::bessel::y1f_coeffs::{Y1_ZEROS, Y1_ZEROS_VALUES, Y1F_COEFFS};
33use crate::common::f_fmla;
34use crate::double_double::DoubleDouble;
35use crate::logs::fast_logf;
36use crate::polyeval::{f_polyeval10, f_polyeval18, f_polyeval19};
37use crate::sincos_reduce::rem2pif_any;
38
39/// Bessel of the second kind of order 1 (Y1)
40///
41/// Max ULP 0.5
42pub fn f_y1f(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            // |x| == 0
48            return f32::NEG_INFINITY;
49        }
50
51        if x.is_infinite() {
52            if x.is_sign_negative() {
53                return f32::NAN;
54            }
55            return 0.;
56        }
57        return x + f32::NAN; // x == NaN
58    }
59
60    let xb = x.to_bits();
61
62    if xb <= 0x424e0000u32 {
63        // x <= 51.5
64        if xb <= 0x40000000u32 {
65            // x <= 2
66            if xb <= 0x3fb5c28fu32 {
67                // x <= 1.42
68                return y1f_near_zero(x);
69            }
70            // transient zone from 1.42 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 y1_transient_area(x);
73        }
74        return y1f_small_argument_path(x);
75    }
76
77    // Exceptions
78    let bx = x.to_bits();
79    if bx == 0x47037a3d {
80        return f32::from_bits(0x2deededb);
81    } else if bx == 0x65ce46e4 {
82        return f32::from_bits(0x9eed85c4);
83    } else if bx == 0x6bf68a7b {
84        return f32::from_bits(0x9dc70a09);
85    } else if bx == 0x76d84625 {
86        return f32::from_bits(0x15d7a68b);
87    } else if bx == 0x7e3dcda0 {
88        return f32::from_bits(0x12b81111);
89    }
90
91    y1f_asympt(x)
92}
93
94/**
95Generated by SageMath:
96Evaluates:
97y2 = -J1(x)*log(x) + 1/x * (1 - sum((-1)^m*(H(m)+H(m-1))/(2^m*m!*(m-1)!)*x^(2*m))
98Y1(x) = 2/pi*(-y2(x)+(euler_gamma - log(2))*J1(x))
99expressed as:
100Y1(x)=log(x)*W1(x) - Z1(x) - 2/(pi*x)
101```python
102from sage.all import *
103
104R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
105x = R.gen()
106N = 16  # Number of terms (adjust as needed)
107gamma = RealField(300)(euler_gamma)
108d2 = RealField(300)(2)
109pi = RealField(300).pi()
110log2 = RealField(300)(2).log()
111
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
115J1_series = j_series(1, x)
116
117def harmony(m):
118    return sum(RealField(300)(1)/RealField(300)(k) for k in range(1, m+1))
119
120def z_series(x):
121    return sum([(-1)**m * (x)**(ZZ(2)*ZZ(m)) / (ZZ(2)**(2*m) * ZZ(m).factorial() * (ZZ(m) - ZZ(1)).factorial()) * (harmony(m) + harmony(m - 1)) for m in range(1, N)])
122
123W1 = d2/pi * J1_series
124Z1 = -(d2/(x*pi) * z_series(x) + d2/pi * gamma * J1_series(x) - d2/pi * log2 * J1_series(x))
125
126# see the series
127print(W0)
128print(Z0)
129```
130See ./notes/bessel_y1_taylor.ipynb for generation
131**/
132#[inline]
133fn y1f_near_zero(x: f32) -> f32 {
134    const W: [u64; 10] = [
135        0x3fd45f306dc9c883,
136        0xbfa45f306dc9c883,
137        0x3f5b2995e7b7b604,
138        0xbf021bb945252402,
139        0x3e9cf9286ea1d337,
140        0xbe2ee7a29824147f,
141        0x3db78be9987d036d,
142        0xbd3ae90af76a4d0f,
143        0x3cb7eb97f85e7d62,
144        0xbc31028e3376648a,
145    ];
146    let dx = x as f64;
147    let x2 = dx * dx;
148    let w0 = f_polyeval10(
149        x2,
150        f64::from_bits(W[0]),
151        f64::from_bits(W[1]),
152        f64::from_bits(W[2]),
153        f64::from_bits(W[3]),
154        f64::from_bits(W[4]),
155        f64::from_bits(W[5]),
156        f64::from_bits(W[6]),
157        f64::from_bits(W[7]),
158        f64::from_bits(W[8]),
159        f64::from_bits(W[9]),
160    ) * dx;
161    const Z: [u64; 10] = [
162        0x3fc91866143cbc8a,
163        0xbfabd3975c75b4a7,
164        0x3f6835b97894be5b,
165        0xbf12c7dbffcde97d,
166        0x3eb0a780ac776eac,
167        0xbe432e5a4ddeea30,
168        0x3dcf0ce34d2066a6,
169        0xbd52a4e1aea45c18,
170        0x3cd1474ade9154ac,
171        0xbc4978ba84f218c0,
172    ];
173    let z0 = f_polyeval10(
174        x2,
175        f64::from_bits(Z[0]),
176        f64::from_bits(Z[1]),
177        f64::from_bits(Z[2]),
178        f64::from_bits(Z[3]),
179        f64::from_bits(Z[4]),
180        f64::from_bits(Z[5]),
181        f64::from_bits(Z[6]),
182        f64::from_bits(Z[7]),
183        f64::from_bits(Z[8]),
184        f64::from_bits(Z[9]),
185    ) * dx;
186    let w_log = fast_logf(x);
187    const TWO_OVER_PI: f64 = f64::from_bits(0x3fe45f306dc9c883);
188    let recip = 1. / dx;
189    let z = f_fmla(w0, w_log, -z0);
190    f_fmla(recip, -TWO_OVER_PI, z) as f32
191}
192
193#[inline]
194fn y1_transient_area(x: f32) -> f32 {
195    let dx = x as f64;
196    // first Y0 bessel zero
197    const ZERO: DoubleDouble =
198        DoubleDouble::from_bit_pair((0xbc8bd1e50d219bfd, 0x400193bed4dff243));
199    let r = (dx - ZERO.hi) - ZERO.lo;
200    /*
201        Poly generated by Wolfram Matematica:
202        <<FunctionApproximations`
203        ClearAll["Global`*"]
204        f[x_]:= BesselY[1,x + 2.1971413260310170351490335626990]
205        {approx,error} = MiniMaxApproximation[f[x],{x,{1.42 - 2.1971413260310170351490335626990, 2 - 2.1971413260310170351490335626990 },17,0},WorkingPrecision->120]
206        poly=error[[1]];
207        coeffs=CoefficientList[poly,x];
208        TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
209    */
210    let p = f_polyeval18(
211        r,
212        f64::from_bits(0x3d9b15a8283b069b),
213        f64::from_bits(0x3fe0aa484455fd09),
214        f64::from_bits(0xbfbe56f80802fa38),
215        f64::from_bits(0xbfa0d2ac9d0409ad),
216        f64::from_bits(0xbf73a619b3551650),
217        f64::from_bits(0x3f7e6c480057ecbb),
218        f64::from_bits(0xbf650dc773a5df4d),
219        f64::from_bits(0x3f531e9ccab7d4da),
220        f64::from_bits(0xbf29b76999169b0e),
221        f64::from_bits(0x3f509c829abceaf7),
222        f64::from_bits(0x3f575aee5697c4d8),
223        f64::from_bits(0x3f63f7f9598be176),
224        f64::from_bits(0x3f67a6ae61541282),
225        f64::from_bits(0x3f665e6d3de19021),
226        f64::from_bits(0x3f5ee8837b9197f6),
227        f64::from_bits(0x3f4e6924f270fd7e),
228        f64::from_bits(0x3f32ca61e5b74925),
229        f64::from_bits(0x3f0725735bc3890b),
230    );
231    p as f32
232}
233
234/// This method on small range searches for nearest zero or extremum.
235/// Then picks stored series expansion at the point end evaluates the poly at the point.
236#[inline]
237fn y1f_small_argument_path(x: f32) -> f32 {
238    let x_abs = x as f64;
239
240    // let avg_step = 51.03 / 33.0;
241    // let inv_step = 1.0 / avg_step;
242    //
243    // println!("inv_step {}", inv_step);
244
245    const INV_STEP: f64 = 0.6466784244562023;
246
247    let fx = x_abs * INV_STEP;
248    const Y1_ZEROS_COUNT: f64 = (Y1_ZEROS.len() - 1) as f64;
249    let idx0 = unsafe { fx.min(Y1_ZEROS_COUNT).to_int_unchecked::<usize>() };
250    let idx1 = unsafe { fx.ceil().min(Y1_ZEROS_COUNT).to_int_unchecked::<usize>() };
251
252    let found_zero0 = DoubleDouble::from_bit_pair(Y1_ZEROS[idx0]);
253    let found_zero1 = DoubleDouble::from_bit_pair(Y1_ZEROS[idx1]);
254
255    let dist0 = (found_zero0.hi - x_abs).abs();
256    let dist1 = (found_zero1.hi - x_abs).abs();
257
258    let (found_zero, idx, dist) = if dist0 < dist1 {
259        (found_zero0, idx0, dist0)
260    } else {
261        (found_zero1, idx1, dist1)
262    };
263
264    if idx == 0 {
265        // Really should not happen here, but if it is then to log expansion
266        return y1f_near_zero(x);
267    }
268
269    // We hit exact zero, value, better to return it directly
270    if dist == 0. {
271        return f64::from_bits(Y1_ZEROS_VALUES[idx]) as f32;
272    }
273
274    let c = &Y1F_COEFFS[idx - 1];
275
276    let r = (x_abs - found_zero.hi) - found_zero.lo;
277
278    let p = f_polyeval19(
279        r,
280        f64::from_bits(c[0]),
281        f64::from_bits(c[1]),
282        f64::from_bits(c[2]),
283        f64::from_bits(c[3]),
284        f64::from_bits(c[4]),
285        f64::from_bits(c[5]),
286        f64::from_bits(c[6]),
287        f64::from_bits(c[7]),
288        f64::from_bits(c[8]),
289        f64::from_bits(c[9]),
290        f64::from_bits(c[10]),
291        f64::from_bits(c[11]),
292        f64::from_bits(c[12]),
293        f64::from_bits(c[13]),
294        f64::from_bits(c[14]),
295        f64::from_bits(c[15]),
296        f64::from_bits(c[16]),
297        f64::from_bits(c[17]),
298        f64::from_bits(c[18]),
299    );
300
301    p as f32
302}
303
304/*
305   Evaluates:
306   Y1 = sqrt(2/(PI*x)) * beta(x) * sin(x - 3*PI/4 - alpha(x))
307
308   Discarding 1/2*PI gives:
309   Y1 = sqrt(2/(PI*x)) * beta(x) * (-cos(x - PI/4 - alpha(x)))
310*/
311#[inline]
312fn y1f_asympt(x: f32) -> f32 {
313    let dx = x as f64;
314
315    let alpha = j1f_asympt_alpha(dx);
316    let beta = j1f_asympt_beta(dx);
317
318    let angle = rem2pif_any(x);
319
320    const SQRT_2_OVER_PI: f64 = f64::from_bits(0x3fe9884533d43651);
321    const MPI_OVER_4: f64 = f64::from_bits(0xbfe921fb54442d18);
322
323    let x0pi34 = MPI_OVER_4 - alpha;
324    let r0 = angle + x0pi34;
325
326    let m_cos = -cos_small(r0);
327
328    let z0 = beta * m_cos;
329    let scale = SQRT_2_OVER_PI * j1f_rsqrt(dx);
330
331    (scale * z0) as f32
332}
333
334#[cfg(test)]
335mod tests {
336    use super::*;
337
338    #[test]
339    fn test_bessel_zero() {
340        assert_eq!(f_y1f(700.76), 0.024876066);
341        assert_eq!(f_y1f(35.76), 0.121432826);
342        assert_eq!(f_y1f(1.76), -0.24787569);
343        assert_eq!(f_y1f(0.87), -0.9030042);
344        assert_eq!(f_y1f(f32::INFINITY), 0.0);
345        assert!(f_y1f(f32::NEG_INFINITY).is_nan());
346        assert!(f_y1f(f32::NAN).is_nan());
347    }
348}
pxfm/bessel/y1f.rs

pxfm/bessel/
y1f.rs
