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

pxfm/bessel/
y1f.rs
