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

pxfm/bessel/
y0.rs
