pxfm/bessel/
y0.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::alpha0::{
30    bessel_0_asympt_alpha, bessel_0_asympt_alpha_fast, bessel_0_asympt_alpha_hard,
31};
32use crate::bessel::beta0::{
33    bessel_0_asympt_beta, bessel_0_asympt_beta_fast, bessel_0_asympt_beta_hard,
34};
35use crate::bessel::i0::bessel_rsqrt_hard;
36use crate::bessel::j0::j0_maclaurin_series;
37use crate::bessel::y0_coeffs::Y0_COEFFS;
38use crate::bessel::y0_coeffs_taylor::Y0_COEFFS_TAYLOR;
39use crate::bessel::y0f_coeffs::{Y0_ZEROS, Y0_ZEROS_VALUES};
40use crate::common::f_fmla;
41use crate::double_double::DoubleDouble;
42use crate::dyadic_float::{DyadicFloat128, DyadicSign};
43use crate::logs::log_dd_fast;
44use crate::polyeval::{f_polyeval12, f_polyeval13, f_polyeval15, f_polyeval22, f_polyeval24};
45use crate::rounding::CpuCeil;
46use crate::sin_helper::{sin_dd_small, sin_dd_small_fast, sin_f128_small};
47use crate::sincos_reduce::{AngleReduced, rem2pi_any, rem2pi_f128};
48
49/// Bessel of the second kind of order 0 (Y0)
50pub fn f_y0(x: f64) -> f64 {
51    let ix = x.to_bits();
52
53    if ix >= 0x7ffu64 << 52 || ix == 0 {
54        // |x| == NaN, x == inf, |x| == 0, x < 0
55        if ix.wrapping_shl(1) == 0 {
56            // |x| == 0
57            return f64::NEG_INFINITY;
58        }
59
60        if x.is_infinite() {
61            if x.is_sign_negative() {
62                return f64::NAN;
63            }
64            return 0.;
65        }
66        return x + f64::NAN; // x == NaN
67    }
68
69    let xb = x.to_bits();
70
71    if xb <= 0x4052d9999999999au64 {
72        // x <= 75.4
73        if xb <= 0x4000000000000000u64 {
74            // x <= 2
75            if xb <= 0x3ff599999999999au64 {
76                // x < 1.35
77                return y0_near_zero_fast(x);
78            }
79            // transient zone from 1.46 to 2 have bad behavior for log poly already,
80            // and not yet good to be easily covered, thus it use its own poly
81            return y0_transient_area_fast(x);
82        }
83        return y0_small_argument_fast(x);
84    }
85
86    y0_asympt_fast(x)
87}
88
89/**
90Generated by SageMath:
91Evaluates:
92Y0(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))
93expressed as:
94Y0(x)=log(x)*W0(x) - Z0(x)
95```python
96from sage.all import *
97
98R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
99x = R.gen()
100N = 10  # Number of terms (adjust as needed)
101gamma = RealField(300)(euler_gamma)
102d2 = RealField(300)(2)
103pi = RealField(300).pi()
104
105# Define J0(x) Taylor expansion at x = 0
106def j_series(n, x):
107    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)])
108
109J0_series = j_series(0, x)
110
111def z_series(x):
112    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)])
113
114W0 = (d2/pi) * J0_series
115Z0 = -gamma * (d2/pi) * J0_series + RealField(300)(2).log() * (d2/pi) * J0_series + (d2/pi) * z_series(x)
116
117# see the series
118print(W0)
119print(Z0)
120```
121**/
122#[inline]
123fn y0_near_zero_fast(x: f64) -> f64 {
124    const W: [(u64, u64); 15] = [
125        (0xbc86b01ec5417056, 0x3fe45f306dc9c883),
126        (0x3c66b01ec5417056, 0xbfc45f306dc9c883),
127        (0xbc26b01ec5417056, 0x3f845f306dc9c883),
128        (0xbbd67fe4a5feb897, 0xbf321bb945252402),
129        (0x3b767fe4a5feb897, 0x3ed21bb945252402),
130        (0xbaf5c2495706f745, 0xbe672db9f21b0f5f),
131        (0x3a90c8209874dfad, 0x3df49a6c656d62ff),
132        (0x3a12921e91b07dd0, 0xbd7ae90af76a4d0f),
133        (0xb992921e91b07dd0, 0x3cfae90af76a4d0f),
134        (0x39089b0d8a9228ca, 0xbc754331c053fdad),
135        (0x3878d321ddfd3c6e, 0x3beb3749ebf0a0dd),
136        (0x37e77548130d809b, 0xbb5cca5ae46eae67),
137        (0xb73a848e7ca1c943, 0x3ac9976d3cd4293f),
138        (0xb6c884706195a054, 0xba336206ff1ce731),
139        (0xb6387a7d2389630d, 0x39995103e9f1818f),
140    ];
141    let x2 = DoubleDouble::from_exact_mult(x, x);
142
143    let w0 = f_polyeval12(
144        x2.hi,
145        f64::from_bits(W[3].1),
146        f64::from_bits(W[4].1),
147        f64::from_bits(W[5].1),
148        f64::from_bits(W[6].1),
149        f64::from_bits(W[7].1),
150        f64::from_bits(W[8].1),
151        f64::from_bits(W[9].1),
152        f64::from_bits(W[10].1),
153        f64::from_bits(W[11].1),
154        f64::from_bits(W[12].1),
155        f64::from_bits(W[13].1),
156        f64::from_bits(W[14].1),
157    );
158
159    let mut w = DoubleDouble::mul_f64_add(x2, w0, DoubleDouble::from_bit_pair(W[2]));
160    w = DoubleDouble::mul_add(x2, w, DoubleDouble::from_bit_pair(W[1]));
161    w = DoubleDouble::mul_add(x2, w, DoubleDouble::from_bit_pair(W[0]));
162
163    const Z: [(u64, u64); 15] = [
164        (0xbc5ddfd831a70821, 0x3fb2e4d699cbd01f),
165        (0xbc6d93e63489aea6, 0xbfc6bbcb41034286),
166        (0xbc1b88525c2e130b, 0x3f9075b1bbf41364),
167        (0x3be097334e26e578, 0xbf41a6206b7b973d),
168        (0x3b51c64a34c78cda, 0x3ee3e99794203bbd),
169        (0xbb1c407b0f5b2805, 0xbe7bce4a600d3ea4),
170        (0xbaa57d1e1e88c9ca, 0x3e0a6ee796b871b6),
171        (0x3a3b6e7030a77899, 0xbd92393d82c6b2e4),
172        (0x397fcfedacb03781, 0x3d131085da82054c),
173        (0xb8e45f51f6118e46, 0xbc8f4ed4b492ebcc),
174        (0xb89bd46046c3c8de, 0x3c04b7ac8a1b15d0),
175        (0x37d1a206fb205e32, 0xbb769201941d0d49),
176        (0x3782f38acbf23993, 0x3ae4987e587ab039),
177        (0x36b691bdabf5672b, 0xba4ff1953e0a7c5b),
178        (0x3636e1c8cd260e18, 0x39b55031dc5e1967),
179    ];
180    let z0 = f_polyeval12(
181        x2.hi,
182        f64::from_bits(Z[3].1),
183        f64::from_bits(Z[4].1),
184        f64::from_bits(Z[5].1),
185        f64::from_bits(Z[6].1),
186        f64::from_bits(Z[7].1),
187        f64::from_bits(Z[8].1),
188        f64::from_bits(Z[9].1),
189        f64::from_bits(Z[10].1),
190        f64::from_bits(Z[11].1),
191        f64::from_bits(Z[12].1),
192        f64::from_bits(Z[13].1),
193        f64::from_bits(Z[14].1),
194    );
195
196    let mut z = DoubleDouble::mul_f64_add(x2, z0, DoubleDouble::from_bit_pair(Z[2]));
197    z = DoubleDouble::mul_add(x2, z, DoubleDouble::from_bit_pair(Z[1]));
198    z = DoubleDouble::mul_add(x2, z, DoubleDouble::from_bit_pair(Z[0]));
199    let w_log = log_dd_fast(x);
200    let p = DoubleDouble::mul_add(w, w_log, -z);
201    let err = f_fmla(
202        p.hi,
203        f64::from_bits(0x3c50000000000000), // 2^-58
204        f64::from_bits(0x3c30000000000000), // 2^-60
205    );
206    let ub = p.hi + (p.lo + err);
207    let lb = p.hi + (p.lo - err);
208    if ub == lb {
209        return p.to_f64();
210    }
211    y0_near_zero(x, w_log)
212}
213
214/**
215Generated by SageMath:
216Evaluates:
217Y0(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))
218expressed as:
219Y0(x)=log(x)*W0(x) - Z0(x)
220```python
221from sage.all import *
222
223R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
224x = R.gen()
225N = 10  # Number of terms (adjust as needed)
226gamma = RealField(300)(euler_gamma)
227d2 = RealField(300)(2)
228pi = RealField(300).pi()
229
230# Define J0(x) Taylor expansion at x = 0
231def j_series(n, x):
232    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)])
233
234J0_series = j_series(0, x)
235
236def z_series(x):
237    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)])
238
239W0 = (d2/pi) * J0_series
240Z0 = -gamma * (d2/pi) * J0_series + RealField(300)(2).log() * (d2/pi) * J0_series + (d2/pi) * z_series(x)
241
242# see the series
243print(W0)
244print(Z0)
245```
246**/
247#[cold]
248#[inline(never)]
249fn y0_near_zero(x: f64, w_log: DoubleDouble) -> f64 {
250    const W: [(u64, u64); 15] = [
251        (0xbc86b01ec5417056, 0x3fe45f306dc9c883),
252        (0x3c66b01ec5417056, 0xbfc45f306dc9c883),
253        (0xbc26b01ec5417056, 0x3f845f306dc9c883),
254        (0xbbd67fe4a5feb897, 0xbf321bb945252402),
255        (0x3b767fe4a5feb897, 0x3ed21bb945252402),
256        (0xbaf5c2495706f745, 0xbe672db9f21b0f5f),
257        (0x3a90c8209874dfad, 0x3df49a6c656d62ff),
258        (0x3a12921e91b07dd0, 0xbd7ae90af76a4d0f),
259        (0xb992921e91b07dd0, 0x3cfae90af76a4d0f),
260        (0x39089b0d8a9228ca, 0xbc754331c053fdad),
261        (0x3878d321ddfd3c6e, 0x3beb3749ebf0a0dd),
262        (0x37e77548130d809b, 0xbb5cca5ae46eae67),
263        (0xb73a848e7ca1c943, 0x3ac9976d3cd4293f),
264        (0xb6c884706195a054, 0xba336206ff1ce731),
265        (0xb6387a7d2389630d, 0x39995103e9f1818f),
266    ];
267    let x2 = DoubleDouble::from_exact_mult(x, x);
268    let w = f_polyeval15(
269        x2,
270        DoubleDouble::from_bit_pair(W[0]),
271        DoubleDouble::from_bit_pair(W[1]),
272        DoubleDouble::from_bit_pair(W[2]),
273        DoubleDouble::from_bit_pair(W[3]),
274        DoubleDouble::from_bit_pair(W[4]),
275        DoubleDouble::from_bit_pair(W[5]),
276        DoubleDouble::from_bit_pair(W[6]),
277        DoubleDouble::from_bit_pair(W[7]),
278        DoubleDouble::from_bit_pair(W[8]),
279        DoubleDouble::from_bit_pair(W[9]),
280        DoubleDouble::from_bit_pair(W[10]),
281        DoubleDouble::from_bit_pair(W[11]),
282        DoubleDouble::from_bit_pair(W[12]),
283        DoubleDouble::from_bit_pair(W[13]),
284        DoubleDouble::from_bit_pair(W[14]),
285    );
286
287    const Z: [(u64, u64); 15] = [
288        (0xbc5ddfd831a70821, 0x3fb2e4d699cbd01f),
289        (0xbc6d93e63489aea6, 0xbfc6bbcb41034286),
290        (0xbc1b88525c2e130b, 0x3f9075b1bbf41364),
291        (0x3be097334e26e578, 0xbf41a6206b7b973d),
292        (0x3b51c64a34c78cda, 0x3ee3e99794203bbd),
293        (0xbb1c407b0f5b2805, 0xbe7bce4a600d3ea4),
294        (0xbaa57d1e1e88c9ca, 0x3e0a6ee796b871b6),
295        (0x3a3b6e7030a77899, 0xbd92393d82c6b2e4),
296        (0x397fcfedacb03781, 0x3d131085da82054c),
297        (0xb8e45f51f6118e46, 0xbc8f4ed4b492ebcc),
298        (0xb89bd46046c3c8de, 0x3c04b7ac8a1b15d0),
299        (0x37d1a206fb205e32, 0xbb769201941d0d49),
300        (0x3782f38acbf23993, 0x3ae4987e587ab039),
301        (0x36b691bdabf5672b, 0xba4ff1953e0a7c5b),
302        (0x3636e1c8cd260e18, 0x39b55031dc5e1967),
303    ];
304    let z = f_polyeval15(
305        x2,
306        DoubleDouble::from_bit_pair(Z[0]),
307        DoubleDouble::from_bit_pair(Z[1]),
308        DoubleDouble::from_bit_pair(Z[2]),
309        DoubleDouble::from_bit_pair(Z[3]),
310        DoubleDouble::from_bit_pair(Z[4]),
311        DoubleDouble::from_bit_pair(Z[5]),
312        DoubleDouble::from_bit_pair(Z[6]),
313        DoubleDouble::from_bit_pair(Z[7]),
314        DoubleDouble::from_bit_pair(Z[8]),
315        DoubleDouble::from_bit_pair(Z[9]),
316        DoubleDouble::from_bit_pair(Z[10]),
317        DoubleDouble::from_bit_pair(Z[11]),
318        DoubleDouble::from_bit_pair(Z[12]),
319        DoubleDouble::from_bit_pair(Z[13]),
320        DoubleDouble::from_bit_pair(Z[14]),
321    );
322    DoubleDouble::mul_add(w, w_log, -z).to_f64()
323}
324
325/**
326Path for transient area between 1.35 to 2.
327**/
328#[inline]
329pub(crate) fn y0_transient_area_fast(x: f64) -> f64 {
330    /**
331    Polynomial generated by Wolfram:
332    ```text
333    <<FunctionApproximations`
334    ClearAll["Global`*"]
335    f[x_]:= BesselY[0,x + 2.1971413260310170351490335626990]
336    {approx,error} = MiniMaxApproximation[f[x],{x,{1.35 - 2.1971413260310170351490335626990, 2 - 2.1971413260310170351490335626990 },27,0},WorkingPrecision->120]
337    poly=error[[1]];
338    coeffs=CoefficientList[poly,x];
339    TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
340    ```
341    **/
342    const C: [(u64, u64); 28] = [
343        (0xbc689e111675434b, 0x3fe0aa48442f014b),
344        (0x396ffb11562e8c70, 0x3cc1806f07aceb3c),
345        (0xbc6156edff56513d, 0xbfd0aa48442f0030),
346        (0xbc278dbff5ee7db4, 0x3fa439fac165db2d),
347        (0x3c1f9463c2023663, 0x3f80d2af4ebc8fc4),
348        (0xbbee09e5733a1236, 0x3f4f716488aebd9c),
349        (0xbbf261a1f255ddf4, 0xbf5444bd2a7e6254),
350        (0x3bd4f543544f1fe7, 0x3f384c300e8674d8),
351        (0xbb96ef8d95fe049b, 0xbf217a0fc83af41e),
352        (0x3ba2be82573ae98d, 0x3f0dbc664048e495),
353        (0xbb942b15646c85f2, 0xbef8522f83e4a3e3),
354        (0x3b7a127725ba4606, 0x3ee775c010ce4146),
355        (0x3ae4a02f0b2a18e2, 0x3e7d1d7cf40f9697),
356        (0x3b8fcf1a3d27236b, 0x3eea9c226c21712d),
357        (0xbb70b3aa0a1e9ffb, 0x3ef8f237831ec74b),
358        (0x3baa6c24261245f3, 0x3f08ebfea3ea469e),
359        (0x3bb5fa1b8c4587c4, 0x3f1474022b90cbda),
360        (0x3b69545db8d098d1, 0x3f1d153dc04c51c0),
361        (0x3bc68eab6520d21b, 0x3f2198a4578cb6ca),
362        (0xbbc255734bc49c8b, 0x3f2212febf7cecdd),
363        (0x3bb8dd02722339f5, 0x3f1f314deec17049),
364        (0x3bbb6ef8f04b26a2, 0x3f1657d051699088),
365        (0x3b878233fbf501dc, 0x3f0a1a422dafcef6),
366        (0xbb73730138d1dbc2, 0x3ef8423cd021f1dd),
367        (0x3b7e2a0d7009d709, 0x3ee145cae37afe1b),
368        (0x3b6af21aeaba4e57, 0x3ec1bc74380f6d7b),
369        (0xbb3607fb9242657f, 0x3e977341fc10cdc8),
370        (0xbac747923880f651, 0x3e5e30218bc1fee3),
371    ];
372
373    const ZERO: DoubleDouble =
374        DoubleDouble::from_bit_pair((0xbc8bd1e50d219bfd, 0x400193bed4dff243));
375
376    let r = DoubleDouble::full_add_f64(-ZERO, x);
377
378    let p0 = f_polyeval24(
379        r.to_f64(),
380        f64::from_bits(C[4].1),
381        f64::from_bits(C[5].1),
382        f64::from_bits(C[6].1),
383        f64::from_bits(C[7].1),
384        f64::from_bits(C[8].1),
385        f64::from_bits(C[9].1),
386        f64::from_bits(C[10].1),
387        f64::from_bits(C[11].1),
388        f64::from_bits(C[12].1),
389        f64::from_bits(C[13].1),
390        f64::from_bits(C[14].1),
391        f64::from_bits(C[15].1),
392        f64::from_bits(C[16].1),
393        f64::from_bits(C[17].1),
394        f64::from_bits(C[18].1),
395        f64::from_bits(C[19].1),
396        f64::from_bits(C[20].1),
397        f64::from_bits(C[21].1),
398        f64::from_bits(C[22].1),
399        f64::from_bits(C[23].1),
400        f64::from_bits(C[24].1),
401        f64::from_bits(C[25].1),
402        f64::from_bits(C[26].1),
403        f64::from_bits(C[27].1),
404    );
405
406    let mut p = DoubleDouble::mul_f64_add(r, p0, DoubleDouble::from_bit_pair(C[3]));
407    p = DoubleDouble::mul_add(p, r, DoubleDouble::from_bit_pair(C[2]));
408    p = DoubleDouble::mul_add(p, r, DoubleDouble::from_bit_pair(C[1]));
409    p = DoubleDouble::mul_add(p, r, DoubleDouble::from_bit_pair(C[0]));
410
411    let err = f_fmla(
412        p.hi,
413        f64::from_bits(0x3c50000000000000), // 2^-58
414        f64::from_bits(0x3b10000000000000), // 2^-78
415    );
416    let ub = p.hi + (p.lo + err);
417    let lb = p.hi + (p.lo - err);
418    if ub != lb {
419        return y0_transient_area_moderate(x);
420    }
421    p.to_f64()
422}
423
424/**
425Path for transient area between 1.35 to 2.
426**/
427fn y0_transient_area_moderate(x: f64) -> f64 {
428    /**
429    Polynomial generated by Wolfram:
430    ```text
431    <<FunctionApproximations`
432    ClearAll["Global`*"]
433    f[x_]:= BesselY[0,x + 2.1971413260310170351490335626990]
434    {approx,error} = MiniMaxApproximation[f[x],{x,{1.35 - 2.1971413260310170351490335626990, 2- 2.1971413260310170351490335626990 },27,0},WorkingPrecision->120]
435    poly=error[[1]];
436    coeffs=CoefficientList[poly,x];
437    TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
438    ```
439    **/
440    const C: [(u64, u64); 28] = [
441        (0xbc689e111675434b, 0x3fe0aa48442f014b),
442        (0x396ffb11562e8c70, 0x3cc1806f07aceb3c),
443        (0xbc6156edff56513d, 0xbfd0aa48442f0030),
444        (0xbc278dbff5ee7db4, 0x3fa439fac165db2d),
445        (0x3c1f9463c2023663, 0x3f80d2af4ebc8fc4),
446        (0xbbee09e5733a1236, 0x3f4f716488aebd9c),
447        (0xbbf261a1f255ddf4, 0xbf5444bd2a7e6254),
448        (0x3bd4f543544f1fe7, 0x3f384c300e8674d8),
449        (0xbb96ef8d95fe049b, 0xbf217a0fc83af41e),
450        (0x3ba2be82573ae98d, 0x3f0dbc664048e495),
451        (0xbb942b15646c85f2, 0xbef8522f83e4a3e3),
452        (0x3b7a127725ba4606, 0x3ee775c010ce4146),
453        (0x3ae4a02f0b2a18e2, 0x3e7d1d7cf40f9697),
454        (0x3b8fcf1a3d27236b, 0x3eea9c226c21712d),
455        (0xbb70b3aa0a1e9ffb, 0x3ef8f237831ec74b),
456        (0x3baa6c24261245f3, 0x3f08ebfea3ea469e),
457        (0x3bb5fa1b8c4587c4, 0x3f1474022b90cbda),
458        (0x3b69545db8d098d1, 0x3f1d153dc04c51c0),
459        (0x3bc68eab6520d21b, 0x3f2198a4578cb6ca),
460        (0xbbc255734bc49c8b, 0x3f2212febf7cecdd),
461        (0x3bb8dd02722339f5, 0x3f1f314deec17049),
462        (0x3bbb6ef8f04b26a2, 0x3f1657d051699088),
463        (0x3b878233fbf501dc, 0x3f0a1a422dafcef6),
464        (0xbb73730138d1dbc2, 0x3ef8423cd021f1dd),
465        (0x3b7e2a0d7009d709, 0x3ee145cae37afe1b),
466        (0x3b6af21aeaba4e57, 0x3ec1bc74380f6d7b),
467        (0xbb3607fb9242657f, 0x3e977341fc10cdc8),
468        (0xbac747923880f651, 0x3e5e30218bc1fee3),
469    ];
470
471    const ZERO: DoubleDouble =
472        DoubleDouble::from_bit_pair((0xbc8bd1e50d219bfd, 0x400193bed4dff243));
473
474    let mut r = DoubleDouble::full_add_f64(-ZERO, x);
475    r = DoubleDouble::from_exact_add(r.hi, r.lo);
476
477    let p0 = f_polyeval13(
478        r.to_f64(),
479        f64::from_bits(C[15].1),
480        f64::from_bits(C[16].1),
481        f64::from_bits(C[17].1),
482        f64::from_bits(C[18].1),
483        f64::from_bits(C[19].1),
484        f64::from_bits(C[20].1),
485        f64::from_bits(C[21].1),
486        f64::from_bits(C[22].1),
487        f64::from_bits(C[23].1),
488        f64::from_bits(C[24].1),
489        f64::from_bits(C[25].1),
490        f64::from_bits(C[26].1),
491        f64::from_bits(C[27].1),
492    );
493
494    let mut p_e = DoubleDouble::mul_f64_add(r, p0, DoubleDouble::from_bit_pair(C[14]));
495    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[13]));
496    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[12]));
497    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[11]));
498    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[10]));
499    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[9]));
500    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[8]));
501    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[7]));
502    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[6]));
503    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[5]));
504    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[4]));
505    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[3]));
506    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[2]));
507    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[1]));
508    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[0]));
509
510    let p = DoubleDouble::from_exact_add(p_e.hi, p_e.lo);
511    let err = f_fmla(
512        p.hi,
513        f64::from_bits(0x3c10000000000000), // 2^-62
514        f64::from_bits(0x3a30000000000000), // 2^-91
515    );
516    let ub = p.hi + (p.lo + err);
517    let lb = p.hi + (p.lo - err);
518    if ub != lb {
519        return y0_transient_area_hard(x);
520    }
521    p.to_f64()
522}
523
524#[cold]
525#[inline(never)]
526fn y0_transient_area_hard(x: f64) -> f64 {
527    const ZERO: DyadicFloat128 = DyadicFloat128 {
528        sign: DyadicSign::Pos,
529        exponent: -126,
530        mantissa: 0x8c9df6a6_ff921721_70d796f3_2017e155_u128,
531    };
532    let r = DyadicFloat128::new_from_f64(x) - ZERO;
533    /*
534    Polynomial generated by Wolfram:
535    ```text
536    <<FunctionApproximations`
537    ClearAll["Global`*"]
538    f[x_]:= BesselY[0,x + 2.1971413260310170351490335626990]
539    {approx,error} = MiniMaxApproximation[f[x],{x,{1.35 - 2.1971413260310170351490335626990, 2- 2.1971413260310170351490335626990 },27,0},WorkingPrecision->120]
540    poly=error[[1]];
541    coeffs=CoefficientList[poly,x];
542    TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
543    ```
544    */
545    static C: [DyadicFloat128; 28] = [
546        DyadicFloat128 {
547            sign: DyadicSign::Pos,
548            exponent: -128,
549            mantissa: 0x85524221_780a573b_0f774c55_e5a946dc_u128,
550        },
551        DyadicFloat128 {
552            sign: DyadicSign::Pos,
553            exponent: -178,
554            mantissa: 0x8c03783d_6759e3ff_622ac5d1_8df27811_u128,
555        },
556        DyadicFloat128 {
557            sign: DyadicSign::Neg,
558            exponent: -129,
559            mantissa: 0x85524221_78018115_6edff565_13cf55ab_u128,
560        },
561        DyadicFloat128 {
562            sign: DyadicSign::Pos,
563            exponent: -132,
564            mantissa: 0xa1cfd60b_2ed96743_9200508c_125cf3d8_u128,
565        },
566        DyadicFloat128 {
567            sign: DyadicSign::Pos,
568            exponent: -134,
569            mantissa: 0x86957a75_e47e21f9_463c2023_663178df_u128,
570        },
571        DyadicFloat128 {
572            sign: DyadicSign::Pos,
573            exponent: -138,
574            mantissa: 0xfb8b2445_75ecdc3e_c35198bd_b9330775_u128,
575        },
576        DyadicFloat128 {
577            sign: DyadicSign::Neg,
578            exponent: -137,
579            mantissa: 0xa225e953_f312a24c_343e4abb_be73338f_u128,
580        },
581        DyadicFloat128 {
582            sign: DyadicSign::Pos,
583            exponent: -139,
584            mantissa: 0xc2618074_33a6c29e_a86a89e3_fcde0075_u128,
585        },
586        DyadicFloat128 {
587            sign: DyadicSign::Neg,
588            exponent: -140,
589            mantissa: 0x8bd07e41_d7a0f05b_be3657f8_126b9715_u128,
590        },
591        DyadicFloat128 {
592            sign: DyadicSign::Pos,
593            exponent: -142,
594            mantissa: 0xede33202_4724aa57_d04ae75d_31ad69cd_u128,
595        },
596        DyadicFloat128 {
597            sign: DyadicSign::Neg,
598            exponent: -143,
599            mantissa: 0xc2917c1f_251f1a85_62ac8d90_be4550ff_u128,
600        },
601        DyadicFloat128 {
602            sign: DyadicSign::Pos,
603            exponent: -144,
604            mantissa: 0xbbae0086_720a31a1_27725ba4_605e0479_u128,
605        },
606        DyadicFloat128 {
607            sign: DyadicSign::Pos,
608            exponent: -151,
609            mantissa: 0xe8ebe7a0_7cb4b852_80bc2ca8_63864ee0_u128,
610        },
611        DyadicFloat128 {
612            sign: DyadicSign::Pos,
613            exponent: -144,
614            mantissa: 0xd4e11361_0b896bf9_e347a4e4_6d642f8e_u128,
615        },
616        DyadicFloat128 {
617            sign: DyadicSign::Pos,
618            exponent: -143,
619            mantissa: 0xc791bc18_f63a577a_62afaf0b_002753e2_u128,
620        },
621        DyadicFloat128 {
622            sign: DyadicSign::Pos,
623            exponent: -142,
624            mantissa: 0xc75ff51f_5234f34d_8484c248_be6beef2_u128,
625        },
626        DyadicFloat128 {
627            sign: DyadicSign::Pos,
628            exponent: -141,
629            mantissa: 0xa3a0115c_865ed2bf_437188b0_f87883dc_u128,
630        },
631        DyadicFloat128 {
632            sign: DyadicSign::Pos,
633            exponent: -141,
634            mantissa: 0xe8a9ee02_628e0019_545db8d0_98d12eff_u128,
635        },
636        DyadicFloat128 {
637            sign: DyadicSign::Pos,
638            exponent: -140,
639            mantissa: 0x8cc522bc_65b652d1_d56ca41a_43537f6a_u128,
640        },
641        DyadicFloat128 {
642            sign: DyadicSign::Pos,
643            exponent: -140,
644            mantissa: 0x9097f5fb_e766e5b5_5196876c_6ea798ce_u128,
645        },
646        DyadicFloat128 {
647            sign: DyadicSign::Pos,
648            exponent: -141,
649            mantissa: 0xf98a6f76_0b824b1b_a04e4467_3ea487e1_u128,
650        },
651        DyadicFloat128 {
652            sign: DyadicSign::Pos,
653            exponent: -141,
654            mantissa: 0xb2be828b_4c84436d_df1e0964_d44f420f_u128,
655        },
656        DyadicFloat128 {
657            sign: DyadicSign::Pos,
658            exponent: -142,
659            mantissa: 0xd0d2116d_7e77b0bc_119fdfa8_0ee2dfee_u128,
660        },
661        DyadicFloat128 {
662            sign: DyadicSign::Pos,
663            exponent: -143,
664            mantissa: 0xc211e681_0f8ee764_67f63971_21f1a199_u128,
665        },
666        DyadicFloat128 {
667            sign: DyadicSign::Pos,
668            exponent: -144,
669            mantissa: 0x8a2e571b_d7f0d9e2_a0d7009d_70969fa2_u128,
670        },
671        DyadicFloat128 {
672            sign: DyadicSign::Pos,
673            exponent: -146,
674            mantissa: 0x8de3a1c0_7b6bdb5e_435d5749_cadf3edd_u128,
675        },
676        DyadicFloat128 {
677            sign: DyadicSign::Pos,
678            exponent: -149,
679            mantissa: 0xbb9a0fe0_866e3d3f_008db7b3_5029fe59_u128,
680        },
681        DyadicFloat128 {
682            sign: DyadicSign::Pos,
683            exponent: -153,
684            mantissa: 0xf1810c5e_0ff717a2_e1b71dfc_26babb9f_u128,
685        },
686    ];
687    let mut z = C[27];
688    for i in (0..27).rev() {
689        z = r * z + C[i];
690    }
691    z.fast_as_f64()
692}
693
694/// This method on small range searches for nearest zero or extremum.
695/// Then picks stored series expansion at the point end evaluates the poly at the point.
696#[inline]
697pub(crate) fn y0_small_argument_fast(x: f64) -> f64 {
698    let x_abs = x;
699
700    // let avg_step = 74.607799 / 47.0;
701    // let inv_step = 1.0 / avg_step;
702
703    const INV_STEP: f64 = 0.6299609508652038;
704
705    let fx = x_abs * INV_STEP;
706    const Y0_ZEROS_COUNT: f64 = (Y0_ZEROS.len() - 1) as f64;
707    let idx0 = unsafe { fx.min(Y0_ZEROS_COUNT).to_int_unchecked::<usize>() };
708    let idx1 = unsafe {
709        fx.cpu_ceil()
710            .min(Y0_ZEROS_COUNT)
711            .to_int_unchecked::<usize>()
712    };
713
714    let found_zero0 = DoubleDouble::from_bit_pair(Y0_ZEROS[idx0]);
715    let found_zero1 = DoubleDouble::from_bit_pair(Y0_ZEROS[idx1]);
716
717    let dist0 = (found_zero0.hi - x_abs).abs();
718    let dist1 = (found_zero1.hi - x_abs).abs();
719
720    let (found_zero, idx, dist) = if dist0 < dist1 {
721        (found_zero0, idx0, dist0)
722    } else {
723        (found_zero1, idx1, dist1)
724    };
725
726    if idx == 0 {
727        return j0_maclaurin_series(x);
728    }
729
730    let is_too_close_to_zero = dist.abs() < 1e-3;
731
732    let c = if is_too_close_to_zero {
733        &Y0_COEFFS_TAYLOR[idx - 1]
734    } else {
735        &Y0_COEFFS[idx - 1]
736    };
737
738    let r = DoubleDouble::full_add_f64(-found_zero, x);
739
740    // We hit exact zero, value, better to return it directly
741    if dist == 0. {
742        return f64::from_bits(Y0_ZEROS_VALUES[idx]);
743    }
744
745    let p = f_polyeval22(
746        r.hi,
747        f64::from_bits(c[6].1),
748        f64::from_bits(c[7].1),
749        f64::from_bits(c[8].1),
750        f64::from_bits(c[9].1),
751        f64::from_bits(c[10].1),
752        f64::from_bits(c[11].1),
753        f64::from_bits(c[12].1),
754        f64::from_bits(c[13].1),
755        f64::from_bits(c[14].1),
756        f64::from_bits(c[15].1),
757        f64::from_bits(c[16].1),
758        f64::from_bits(c[17].1),
759        f64::from_bits(c[18].1),
760        f64::from_bits(c[19].1),
761        f64::from_bits(c[20].1),
762        f64::from_bits(c[21].1),
763        f64::from_bits(c[22].1),
764        f64::from_bits(c[23].1),
765        f64::from_bits(c[24].1),
766        f64::from_bits(c[25].1),
767        f64::from_bits(c[26].1),
768        f64::from_bits(c[27].1),
769    );
770
771    let mut z = DoubleDouble::mul_f64_add(r, p, DoubleDouble::from_bit_pair(c[5]));
772    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[4]));
773    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[3]));
774    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[2]));
775    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[1]));
776    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[0]));
777    let p = z;
778    let err = f_fmla(
779        p.hi,
780        f64::from_bits(0x3c60000000000000), // 2^-57
781        f64::from_bits(0x3c20000000000000), // 2^-61
782    );
783    let ub = p.hi + (p.lo + err);
784    let lb = p.hi + (p.lo - err);
785    if ub != lb {
786        return y0_small_argument_moderate(r, c);
787    }
788    z.to_f64()
789}
790
791/// This method on small range searches for nearest zero or extremum.
792/// Then picks stored series expansion at the point end evaluates the poly at the point.
793fn y0_small_argument_moderate(r: DoubleDouble, c: &[(u64, u64); 28]) -> f64 {
794    let c0 = &c[15..];
795
796    let p0 = f_polyeval13(
797        r.to_f64(),
798        f64::from_bits(c0[0].1),
799        f64::from_bits(c0[1].1),
800        f64::from_bits(c0[2].1),
801        f64::from_bits(c0[3].1),
802        f64::from_bits(c0[4].1),
803        f64::from_bits(c0[5].1),
804        f64::from_bits(c0[6].1),
805        f64::from_bits(c0[7].1),
806        f64::from_bits(c0[8].1),
807        f64::from_bits(c0[9].1),
808        f64::from_bits(c0[10].1),
809        f64::from_bits(c0[11].1),
810        f64::from_bits(c0[12].1),
811    );
812
813    let mut p_e = DoubleDouble::mul_f64_add(r, p0, DoubleDouble::from_bit_pair(c[14]));
814    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[13]));
815    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[12]));
816    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[11]));
817    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[10]));
818    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[9]));
819    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[8]));
820    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[7]));
821    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[6]));
822    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[5]));
823    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[4]));
824    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[3]));
825    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[2]));
826    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[1]));
827    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[0]));
828
829    let p = DoubleDouble::from_exact_add(p_e.hi, p_e.lo);
830    let err = f_fmla(
831        p.hi,
832        f64::from_bits(0x3c30000000000000), // 2^-60
833        f64::from_bits(0x3a30000000000000), // 2^-91
834    );
835    let ub = p.hi + (p.lo + err);
836    let lb = p.hi + (p.lo - err);
837    if ub != lb {
838        return y0_small_argument_hard(r, c);
839    }
840    p.to_f64()
841}
842
843#[cold]
844#[inline(never)]
845fn y0_small_argument_hard(r: DoubleDouble, c: &[(u64, u64); 28]) -> f64 {
846    let mut p = DoubleDouble::from_bit_pair(c[27]);
847    for i in (0..27).rev() {
848        p = DoubleDouble::mul_add(r, p, DoubleDouble::from_bit_pair(c[i]));
849        p = DoubleDouble::from_exact_add(p.hi, p.lo);
850    }
851    p.to_f64()
852}
853
854/*
855   Evaluates:
856   Y0 = sqrt(2/(PI*x)) * beta(x) * sin(x - PI/4 - alpha(x))
857*/
858#[inline]
859pub(crate) fn y0_asympt_fast(x: f64) -> f64 {
860    const SQRT_2_OVER_PI: DoubleDouble = DoubleDouble::new(
861        f64::from_bits(0xbc8cbc0d30ebfd15),
862        f64::from_bits(0x3fe9884533d43651),
863    );
864    const MPI_OVER_4: DoubleDouble = DoubleDouble::new(
865        f64::from_bits(0xbc81a62633145c07),
866        f64::from_bits(0xbfe921fb54442d18),
867    );
868
869    let recip = if x.to_bits() > 0x7fd000000000000u64 {
870        DoubleDouble::quick_mult_f64(DoubleDouble::from_exact_safe_div(4.0, x), 0.25)
871    } else {
872        DoubleDouble::from_recip(x)
873    };
874
875    let alpha = bessel_0_asympt_alpha_fast(recip);
876    let beta = bessel_0_asympt_beta_fast(recip);
877
878    let AngleReduced { angle } = rem2pi_any(x);
879
880    // Without full subtraction cancellation happens sometimes
881    let x0pi34 = DoubleDouble::full_dd_sub(MPI_OVER_4, alpha);
882    let r0 = DoubleDouble::full_dd_add(angle, x0pi34);
883
884    let m_cos = sin_dd_small_fast(r0);
885    let z0 = DoubleDouble::quick_mult(beta, m_cos);
886    let r_sqrt = DoubleDouble::from_rsqrt_fast(x);
887    let scale = DoubleDouble::quick_mult(SQRT_2_OVER_PI, r_sqrt);
888    let r = DoubleDouble::quick_mult(scale, z0);
889    let p = DoubleDouble::from_exact_add(r.hi, r.lo);
890    let err = f_fmla(
891        p.hi,
892        f64::from_bits(0x3c40000000000000), // 2^-59
893        f64::from_bits(0x3c10000000000000), // 2^-62
894    );
895    let ub = p.hi + (p.lo + err);
896    let lb = p.hi + (p.lo - err);
897
898    if ub == lb {
899        return p.to_f64();
900    }
901    y0_asympt(x, recip, r_sqrt, angle)
902}
903
904/*
905   Evaluates:
906   Y0 = sqrt(2/(PI*x)) * beta(x) * sin(x - PI/4 - alpha(x))
907*/
908fn y0_asympt(x: f64, recip: DoubleDouble, r_sqrt: DoubleDouble, angle: DoubleDouble) -> f64 {
909    const SQRT_2_OVER_PI: DoubleDouble = DoubleDouble::new(
910        f64::from_bits(0xbc8cbc0d30ebfd15),
911        f64::from_bits(0x3fe9884533d43651),
912    );
913    const MPI_OVER_4: DoubleDouble = DoubleDouble::new(
914        f64::from_bits(0xbc81a62633145c07),
915        f64::from_bits(0xbfe921fb54442d18),
916    );
917
918    let alpha = bessel_0_asympt_alpha(recip);
919    let beta = bessel_0_asympt_beta(recip);
920
921    // Without full subtraction cancellation happens sometimes
922    let x0pi34 = DoubleDouble::full_dd_sub(MPI_OVER_4, alpha);
923    let r0 = DoubleDouble::full_dd_add(angle, x0pi34);
924
925    let m_cos = sin_dd_small(r0);
926    let z0 = DoubleDouble::quick_mult(beta, m_cos);
927    let scale = DoubleDouble::quick_mult(SQRT_2_OVER_PI, r_sqrt);
928    let r = DoubleDouble::quick_mult(scale, z0);
929    let p = DoubleDouble::from_exact_add(r.hi, r.lo);
930    let err = f_fmla(
931        p.hi,
932        f64::from_bits(0x3c30000000000000), // 2^-60
933        f64::from_bits(0x3a20000000000000), // 2^-93
934    );
935    let ub = p.hi + (p.lo + err);
936    let lb = p.hi + (p.lo - err);
937
938    if ub == lb {
939        return p.to_f64();
940    }
941    y0_asympt_hard(x)
942}
943
944/*
945   Evaluates:
946   Y0 = sqrt(2/(PI*x)) * beta(x) * sin(x - PI/4 - alpha(x))
947*/
948#[cold]
949#[inline(never)]
950fn y0_asympt_hard(x: f64) -> f64 {
951    const SQRT_2_OVER_PI: DyadicFloat128 = DyadicFloat128 {
952        sign: DyadicSign::Pos,
953        exponent: -128,
954        mantissa: 0xcc42299e_a1b28468_7e59e280_5d5c7180_u128,
955    };
956
957    const MPI_OVER_4: DyadicFloat128 = DyadicFloat128 {
958        sign: DyadicSign::Neg,
959        exponent: -128,
960        mantissa: 0xc90fdaa2_2168c234_c4c6628b_80dc1cd1_u128,
961    };
962
963    let x_dyadic = DyadicFloat128::new_from_f64(x);
964    let recip = DyadicFloat128::accurate_reciprocal(x);
965
966    let alpha = bessel_0_asympt_alpha_hard(recip);
967    let beta = bessel_0_asympt_beta_hard(recip);
968
969    let angle = rem2pi_f128(x_dyadic);
970
971    let x0pi34 = MPI_OVER_4 - alpha;
972    let r0 = angle + x0pi34;
973
974    let m_sin = sin_f128_small(r0);
975
976    let z0 = beta * m_sin;
977    let r_sqrt = bessel_rsqrt_hard(x, recip);
978    let scale = SQRT_2_OVER_PI * r_sqrt;
979    let p = scale * z0;
980    p.fast_as_f64()
981}
982
983#[cfg(test)]
984mod tests {
985    use super::*;
986
987    #[test]
988    fn test_y0() {
989        //ULP should be less than 0.5, but it was 1017.1969361449036, on 3 result 0.37685001001284685, using f_y0 and MPFR 0.3768500100127904
990        assert_eq!(f_y0(3.), 0.3768500100127904);
991        assert_eq!(f_y0(0.906009703874588), 0.01085796448629276);
992        assert_eq!(f_y0(80.), -0.05562033908977);
993        assert_eq!(f_y0(5.), -0.30851762524903376);
994        assert_eq!(
995            f_y0(f64::from_bits(0x3fec982eb8d417ea)),
996            -0.000000000000000023389279284062102
997        );
998        assert!(f_y0(f64::NAN).is_nan());
999        assert_eq!(f_y0(f64::INFINITY), 0.);
1000        assert!(f_y0(f64::NEG_INFINITY).is_nan());
1001    }
1002
1003    #[test]
1004    fn test_y0_edge_values() {
1005        assert_eq!(f_y0(0.8900000000138676), -0.0031519646708080126);
1006        assert_eq!(f_y0(0.8900000000409116), -0.0031519646469294936);
1007        assert_eq!(f_y0(98.1760435789366), 0.0000000000000056889416242533015);
1008        assert_eq!(
1009            f_y0(91.8929453121571802176),
1010            -0.00000000000000007281665706677893
1011        );
1012        assert_eq!(
1013            f_y0(f64::from_bits(0x6e7c1d741dc52512u64)),
1014            f64::from_bits(0x2696f860815bc669)
1015        );
1016        assert_eq!(f_y0(f64::from_bits(0x3e04cdee58a47edd)), -13.58605001628649);
1017        assert_eq!(
1018            f_y0(0.89357696627916749),
1019            -0.000000000000000023389279284062102
1020        );
1021    }
1022}
pxfm/bessel/y0.rs

pxfm/bessel/
y0.rs
