pxfm/bessel/
y1.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::alpha1::{
30    bessel_1_asympt_alpha, bessel_1_asympt_alpha_fast, bessel_1_asympt_alpha_hard,
31};
32use crate::bessel::beta1::{
33    bessel_1_asympt_beta, bessel_1_asympt_beta_fast, bessel_1_asympt_beta_hard,
34};
35use crate::bessel::i0::bessel_rsqrt_hard;
36use crate::bessel::y1_coeffs::Y1_COEFFS_REMEZ;
37use crate::bessel::y1_coeffs_taylor::Y1_COEFFS;
38use crate::bessel::y1f_coeffs::{Y1_ZEROS, Y1_ZEROS_VALUES};
39use crate::common::f_fmla;
40use crate::double_double::DoubleDouble;
41use crate::dyadic_float::{DyadicFloat128, DyadicSign};
42use crate::logs::log_dd_fast;
43use crate::polyeval::{f_polyeval12, f_polyeval13, f_polyeval15, f_polyeval22, f_polyeval24};
44use crate::sin_helper::{cos_dd_small, cos_dd_small_fast, cos_f128_small};
45use crate::sincos_reduce::{AngleReduced, rem2pi_any, rem2pi_f128};
46
47/// Bessel of the second kind of order 1 ( Y1 )
48pub fn f_y1(x: f64) -> f64 {
49    let ix = x.to_bits();
50
51    if ix >= 0x7ffu64 << 52 || ix == 0 {
52        // |x| == NaN, x == inf, |x| == 0, x < 0
53        if ix.wrapping_shl(1) == 0 {
54            // |x| == 0
55            return f64::NEG_INFINITY;
56        }
57
58        if x.is_infinite() {
59            if x.is_sign_negative() {
60                return f64::NAN;
61            }
62            return 0.;
63        }
64
65        return x + f64::NAN; // x == NaN
66    }
67
68    let xb = x.to_bits();
69
70    if xb <= 0x4049c00000000000u64 {
71        // x <= 51.5
72        if xb <= 0x4000000000000000u64 {
73            // x <= 2
74            if xb <= 0x3ff75c28f5c28f5cu64 {
75                // x <= 1.46
76                return y1_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 y1_transient_zone_fast(x);
81        }
82
83        return y1_small_argument_fast(x);
84    }
85
86    y1_asympt_fast(x)
87}
88
89/**
90Generated by SageMath:
91Evaluates:
92y2 = -J1(x)*log(x) + 1/x * (1 - sum((-1)^m*(H(m)+H(m-1))/(2^m*m!*(m-1)!)*x^(2*m))
93Y1(x) = 2/pi*(-y2(x)+(euler_gamma - log(2))*J1(x))
94expressed as:
95Y1(x)=log(x)*W1(x) - Z1(x) - 2/(pi*x)
96```python
97from sage.all import *
98
99R = LaurentSeriesRing(RealField(300), 'x', default_prec=300)
100x = R.gen()
101N = 16  # Number of terms (adjust as needed)
102gamma = RealField(300)(euler_gamma)
103d2 = RealField(300)(2)
104pi = RealField(300).pi()
105log2 = RealField(300)(2).log()
106
107def j_series(n, x):
108    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)])
109
110J1_series = j_series(1, x)
111
112def harmony(m):
113    return sum(RealField(300)(1)/RealField(300)(k) for k in range(1, m+1))
114
115def z_series(x):
116    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)])
117
118W1 = d2/pi * J1_series
119Z1 = -(d2/(x*pi) * z_series(x) + d2/pi * gamma * J1_series(x) - d2/pi * log2 * J1_series(x))
120
121def y1_full(x):
122    return d2/pi * (J1_series(x) * x.log() - 1/x * ( 1 - z_series(x)) + (gamma - log2) * J1_series(x))
123
124# see the series
125print(W0)
126print(Z0)
127```
128See ./notes/bessel_y1_taylor.ipynb for generation
129**/
130#[inline]
131fn y1_near_zero_fast(x: f64) -> f64 {
132    const W: [(u64, u64); 15] = [
133        (0xbc76b01ec5417056, 0x3fd45f306dc9c883),
134        (0x3c46b01ec5417056, 0xbfa45f306dc9c883),
135        (0xbbfe40290701eb1e, 0x3f5b2995e7b7b604),
136        (0xbba67fe4a5feb897, 0xbf021bb945252402),
137        (0xbb0334914cdd2ba9, 0x3e9cf9286ea1d337),
138        (0x3aab4f3c6d42c1f4, 0xbe2ee7a29824147f),
139        (0xba407fb57ef4dc2c, 0x3db78be9987d036d),
140        (0x39d2921e91b07dd0, 0xbd3ae90af76a4d0f),
141        (0x395a28c8620dc90e, 0x3cb7eb97f85e7d62),
142        (0xb8cf83f52abe45c5, 0xbc31028e3376648a),
143        (0xb8441050c68ca435, 0x3ba3cb1e7d0c17e7),
144        (0xb7ab072548a1aa43, 0xbb133191ed9f1eef),
145        (0xb6f05192c2d9b6ee, 0x3a7f7f4b5e8ef7b0),
146        (0x367ad65afe306d57, 0xb9e626e36cb3515d),
147        (0xb5ea1c4136f8f230, 0x394b01153dce6810),
148    ];
149    let x2 = DoubleDouble::from_exact_mult(x, x);
150    let w0 = f_polyeval12(
151        x2.hi,
152        f64::from_bits(W[3].1),
153        f64::from_bits(W[4].1),
154        f64::from_bits(W[5].1),
155        f64::from_bits(W[6].1),
156        f64::from_bits(W[7].1),
157        f64::from_bits(W[8].1),
158        f64::from_bits(W[9].1),
159        f64::from_bits(W[10].1),
160        f64::from_bits(W[11].1),
161        f64::from_bits(W[12].1),
162        f64::from_bits(W[13].1),
163        f64::from_bits(W[14].1),
164    );
165
166    let mut w = DoubleDouble::mul_f64_add(x2, w0, DoubleDouble::from_bit_pair(W[2]));
167    w = DoubleDouble::mul_add(x2, w, DoubleDouble::from_bit_pair(W[1]));
168    w = DoubleDouble::mul_add(x2, w, DoubleDouble::from_bit_pair(W[0]));
169    w = DoubleDouble::quick_mult_f64(w, x);
170
171    const Z: [(u64, u64); 15] = [
172        (0x3c61d7eb2e54cda1, 0x3fc91866143cbc8a),
173        (0xbc2f9f7a0ce54a40, 0xbfabd3975c75b4a7),
174        (0xbbf7659313f45e8c, 0x3f6835b97894be5b),
175        (0x3b9cbcd40f1be7b9, 0xbf12c7dbffcde97d),
176        (0xbb495d78778645b4, 0x3eb0a780ac776eac),
177        (0xbae15be86455c1ab, 0xbe432e5a4ddeea30),
178        (0xba5ad966c12f1e3c, 0x3dcf0ce34d2066a6),
179        (0x39e9717155dc7521, 0xbd52a4e1aea45c18),
180        (0x394f447fe5de1290, 0x3cd1474ade9154ac),
181        (0xb8e1699d9009a7fc, 0xbc4978ba84f218c0),
182        (0xb8505502096ead17, 0x3bbe9598c016378b),
183        (0x37942b6c36b2c5f1, 0xbb2e7e5fcfc4b7b1),
184        (0x37210853b78bd08a, 0x3a99a6c1266c116d),
185        (0xb686c9639c9d976e, 0xba02738998fe7337),
186        (0xb603b739ee04b9fe, 0x3966f58cd41b6d08),
187    ];
188    let z0 = f_polyeval12(
189        x2.hi,
190        f64::from_bits(Z[3].1),
191        f64::from_bits(Z[4].1),
192        f64::from_bits(Z[5].1),
193        f64::from_bits(Z[6].1),
194        f64::from_bits(Z[7].1),
195        f64::from_bits(Z[8].1),
196        f64::from_bits(Z[9].1),
197        f64::from_bits(Z[10].1),
198        f64::from_bits(Z[11].1),
199        f64::from_bits(Z[12].1),
200        f64::from_bits(Z[13].1),
201        f64::from_bits(Z[14].1),
202    );
203
204    let mut z = DoubleDouble::mul_f64_add(x2, z0, DoubleDouble::from_bit_pair(Z[2]));
205    z = DoubleDouble::mul_add(x2, z, DoubleDouble::from_bit_pair(Z[1]));
206    z = DoubleDouble::mul_add(x2, z, DoubleDouble::from_bit_pair(Z[0]));
207    z = DoubleDouble::quick_mult_f64(z, x);
208
209    let w_log = log_dd_fast(x);
210
211    const MINUS_TWO_OVER_PI: DoubleDouble =
212        DoubleDouble::from_bit_pair((0x3c86b01ec5417056, 0xbfe45f306dc9c883));
213
214    let m_two_over_pi_div_x: DoubleDouble;
215    #[cfg(any(
216        all(
217            any(target_arch = "x86", target_arch = "x86_64"),
218            target_feature = "fma"
219        ),
220        all(target_arch = "aarch64", target_feature = "neon")
221    ))]
222    {
223        m_two_over_pi_div_x = DoubleDouble::div_dd_f64(MINUS_TWO_OVER_PI, x)
224    }
225    #[cfg(not(any(
226        all(
227            any(target_arch = "x86", target_arch = "x86_64"),
228            target_feature = "fma"
229        ),
230        all(target_arch = "aarch64", target_feature = "neon")
231    )))]
232    {
233        use crate::double_double::two_product_compatible;
234        m_two_over_pi_div_x = if two_product_compatible(x) {
235            DoubleDouble::div_dd_f64(MINUS_TWO_OVER_PI, x)
236        } else {
237            DoubleDouble::div_safe_dd_f64(MINUS_TWO_OVER_PI, x)
238        };
239    }
240    if m_two_over_pi_div_x.hi.is_infinite() {
241        return f64::NEG_INFINITY;
242    }
243
244    let zvp = DoubleDouble::mul_add(w, w_log, -z);
245    let p = DoubleDouble::add(m_two_over_pi_div_x, zvp);
246    let err = f_fmla(
247        p.hi,
248        f64::from_bits(0x3c30000000000000), // 2^-60
249        f64::from_bits(0x3be0000000000000), // 2^-65
250    );
251    let ub = p.hi + (p.lo + err);
252    let lb = p.hi + (p.lo - err);
253    if ub == lb {
254        return p.to_f64();
255    }
256    y1_near_zero(x, w_log)
257}
258
259/**
260Generated by SageMath:
261Evaluates:
262y2 = -J1(x)*log(x) + 1/x * (1 - sum((-1)^m*(H(m)+H(m-1))/(2^m*m!*(m-1)!)*x^(2*m))
263Y1(x) = 2/pi*(-y2(x)+(euler_gamma - log(2))*J1(x))
264expressed as:
265Y1(x)=log(x)*W1(x) - Z1(x) - 2/(pi*x)
266```python
267from sage.all import *
268
269R = LaurentSeriesRing(RealField(300), 'x', default_prec=300)
270x = R.gen()
271N = 16  # Number of terms (adjust as needed)
272gamma = RealField(300)(euler_gamma)
273d2 = RealField(300)(2)
274pi = RealField(300).pi()
275log2 = RealField(300)(2).log()
276
277def j_series(n, x):
278    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)])
279
280J1_series = j_series(1, x)
281
282def harmony(m):
283    return sum(RealField(300)(1)/RealField(300)(k) for k in range(1, m+1))
284
285def z_series(x):
286    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)])
287
288W1 = d2/pi * J1_series
289Z1 = -(d2/(x*pi) * z_series(x) + d2/pi * gamma * J1_series(x) - d2/pi * log2 * J1_series(x))
290
291def y1_full(x):
292    return d2/pi * (J1_series(x) * x.log() - 1/x * ( 1 - z_series(x)) + (gamma - log2) * J1_series(x))
293
294# see the series
295print(W0)
296print(Z0)
297```
298See ./notes/bessel_y1_taylor.ipynb for generation
299**/
300#[cold]
301#[inline(never)]
302fn y1_near_zero(x: f64, w_log: DoubleDouble) -> f64 {
303    const W: [(u64, u64); 15] = [
304        (0xbc76b01ec5417056, 0x3fd45f306dc9c883),
305        (0x3c46b01ec5417056, 0xbfa45f306dc9c883),
306        (0xbbfe40290701eb1e, 0x3f5b2995e7b7b604),
307        (0xbba67fe4a5feb897, 0xbf021bb945252402),
308        (0xbb0334914cdd2ba9, 0x3e9cf9286ea1d337),
309        (0x3aab4f3c6d42c1f4, 0xbe2ee7a29824147f),
310        (0xba407fb57ef4dc2c, 0x3db78be9987d036d),
311        (0x39d2921e91b07dd0, 0xbd3ae90af76a4d0f),
312        (0x395a28c8620dc90e, 0x3cb7eb97f85e7d62),
313        (0xb8cf83f52abe45c5, 0xbc31028e3376648a),
314        (0xb8441050c68ca435, 0x3ba3cb1e7d0c17e7),
315        (0xb7ab072548a1aa43, 0xbb133191ed9f1eef),
316        (0xb6f05192c2d9b6ee, 0x3a7f7f4b5e8ef7b0),
317        (0x367ad65afe306d57, 0xb9e626e36cb3515d),
318        (0xb5ea1c4136f8f230, 0x394b01153dce6810),
319    ];
320    let x2 = DoubleDouble::from_exact_mult(x, x);
321    let mut w = f_polyeval15(
322        x2,
323        DoubleDouble::from_bit_pair(W[0]),
324        DoubleDouble::from_bit_pair(W[1]),
325        DoubleDouble::from_bit_pair(W[2]),
326        DoubleDouble::from_bit_pair(W[3]),
327        DoubleDouble::from_bit_pair(W[4]),
328        DoubleDouble::from_bit_pair(W[5]),
329        DoubleDouble::from_bit_pair(W[6]),
330        DoubleDouble::from_bit_pair(W[7]),
331        DoubleDouble::from_bit_pair(W[8]),
332        DoubleDouble::from_bit_pair(W[9]),
333        DoubleDouble::from_bit_pair(W[10]),
334        DoubleDouble::from_bit_pair(W[11]),
335        DoubleDouble::from_bit_pair(W[12]),
336        DoubleDouble::from_bit_pair(W[13]),
337        DoubleDouble::from_bit_pair(W[14]),
338    );
339    w = DoubleDouble::quick_mult_f64(w, x);
340
341    const Z: [(u64, u64); 15] = [
342        (0x3c61d7eb2e54cda1, 0x3fc91866143cbc8a),
343        (0xbc2f9f7a0ce54a40, 0xbfabd3975c75b4a7),
344        (0xbbf7659313f45e8c, 0x3f6835b97894be5b),
345        (0x3b9cbcd40f1be7b9, 0xbf12c7dbffcde97d),
346        (0xbb495d78778645b4, 0x3eb0a780ac776eac),
347        (0xbae15be86455c1ab, 0xbe432e5a4ddeea30),
348        (0xba5ad966c12f1e3c, 0x3dcf0ce34d2066a6),
349        (0x39e9717155dc7521, 0xbd52a4e1aea45c18),
350        (0x394f447fe5de1290, 0x3cd1474ade9154ac),
351        (0xb8e1699d9009a7fc, 0xbc4978ba84f218c0),
352        (0xb8505502096ead17, 0x3bbe9598c016378b),
353        (0x37942b6c36b2c5f1, 0xbb2e7e5fcfc4b7b1),
354        (0x37210853b78bd08a, 0x3a99a6c1266c116d),
355        (0xb686c9639c9d976e, 0xba02738998fe7337),
356        (0xb603b739ee04b9fe, 0x3966f58cd41b6d08),
357    ];
358    let mut z = f_polyeval15(
359        x2,
360        DoubleDouble::from_bit_pair(Z[0]),
361        DoubleDouble::from_bit_pair(Z[1]),
362        DoubleDouble::from_bit_pair(Z[2]),
363        DoubleDouble::from_bit_pair(Z[3]),
364        DoubleDouble::from_bit_pair(Z[4]),
365        DoubleDouble::from_bit_pair(Z[5]),
366        DoubleDouble::from_bit_pair(Z[6]),
367        DoubleDouble::from_bit_pair(Z[7]),
368        DoubleDouble::from_bit_pair(Z[8]),
369        DoubleDouble::from_bit_pair(Z[9]),
370        DoubleDouble::from_bit_pair(Z[10]),
371        DoubleDouble::from_bit_pair(Z[11]),
372        DoubleDouble::from_bit_pair(Z[12]),
373        DoubleDouble::from_bit_pair(Z[13]),
374        DoubleDouble::from_bit_pair(Z[14]),
375    );
376    z = DoubleDouble::quick_mult_f64(z, x);
377
378    const MINUS_TWO_OVER_PI: DoubleDouble =
379        DoubleDouble::from_bit_pair((0x3c86b01ec5417056, 0xbfe45f306dc9c883));
380
381    let m_two_over_pi_div_x: DoubleDouble;
382    #[cfg(any(
383        all(
384            any(target_arch = "x86", target_arch = "x86_64"),
385            target_feature = "fma"
386        ),
387        all(target_arch = "aarch64", target_feature = "neon")
388    ))]
389    {
390        m_two_over_pi_div_x = DoubleDouble::div_dd_f64(MINUS_TWO_OVER_PI, x)
391    }
392    #[cfg(not(any(
393        all(
394            any(target_arch = "x86", target_arch = "x86_64"),
395            target_feature = "fma"
396        ),
397        all(target_arch = "aarch64", target_feature = "neon")
398    )))]
399    {
400        use crate::double_double::two_product_compatible;
401        m_two_over_pi_div_x = if two_product_compatible(x) {
402            DoubleDouble::div_dd_f64(MINUS_TWO_OVER_PI, x)
403        } else {
404            DoubleDouble::div_safe_dd_f64(MINUS_TWO_OVER_PI, x)
405        };
406    }
407    if m_two_over_pi_div_x.hi.is_infinite() {
408        return f64::NEG_INFINITY;
409    }
410
411    let zvp = DoubleDouble::mul_add(w, w_log, -z);
412    DoubleDouble::full_dd_add(m_two_over_pi_div_x, zvp).to_f64()
413}
414
415#[inline]
416fn y1_transient_zone_fast(x: f64) -> f64 {
417    /*
418    <<FunctionApproximations`
419    ClearAll["Global`*"]
420    f[x_]:= BesselY[1,x + 2.1971413260310170351490335626990]
421    {approx,error} =MiniMaxApproximation[f[x],{x,{1.42 -  2.1971413260310170351490335626990, 2-  2.1971413260310170351490335626990 },27,0},WorkingPrecision->120]
422    poly=error[[1]];
423    coeffs=CoefficientList[poly,x];
424    TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
425     */
426    const C: [(u64, u64); 28] = [
427        (0x38d23216c8eac2ee, 0x3c631de77e55e9b0),
428        (0x3c77f3a6718f8e03, 0x3fe0aa48442f0150),
429        (0xbc52b77d09be8679, 0xbfbe56f82217b33a),
430        (0x3c237af39fc9d759, 0xbfa0d2af4e922afe),
431        (0xbc08241e95389bc3, 0xbf73a6dec2f9f739),
432        (0x3c18eac6a35acd63, 0x3f7e671c82a1c956),
433        (0x3be0a0f3c8908083, 0xbf65429d5dbc3bb0),
434        (0x3be217fa58861600, 0x3f517abad71c26c0),
435        (0x3bebc327d6c65514, 0xbf40b28b4ef33a56),
436        (0xbbcb597d6bdd5992, 0x3f2ef0d150ec9934),
437        (0x3bb1a480de696e07, 0xbf1c0758a86844be),
438        (0xbb9c6ea6352ab84e, 0x3f0b7c88b58d9ef3),
439        (0x3b7dbb78dfd868a9, 0xbee97ef9519bdcfd),
440        (0xbb8d4519030f499d, 0x3f04815f08e7ad5e),
441        (0x3bbbd70e1480e260, 0x3f10f348483b57bc),
442        (0x3b7117cf4d2b6f3c, 0x3f231914389bb1bb),
443        (0x3b8ca48beaf6a58d, 0x3f30b29e838345b4),
444        (0x3bccfac65ce17cf9, 0x3f39a69e98f61897),
445        (0xbbeae7e3065b09c9, 0x3f40b9511666fcf0),
446        (0xbbe5cbddf691e7e6, 0x3f428cd8388e634b),
447        (0xbbd91372412d1e1b, 0x3f414ba048d9e1d5),
448        (0xbbb0781a70c6f715, 0x3f3acdcf66f1de95),
449        (0xbba3ae83fd425494, 0x3f30f44ae6620bba),
450        (0x3bc001d75da77b74, 0x3f21154a0a1f2161),
451        (0xbb91c9afb1a1b874, 0x3f0a687b664cbac6),
452        (0x3b8e0a06b9444963, 0x3eed7c3cbb4ba5d8),
453        (0x3b4f0e9dfc915934, 0x3ec53ca23fdd0999),
454        (0xbb0409258f8ffca8, 0x3e8de620acb51b2d),
455    ];
456
457    // this poly relative to first zero
458    const ZERO: DoubleDouble =
459        DoubleDouble::from_bit_pair((0xbc8bd1e50d219bfd, 0x400193bed4dff243));
460
461    let mut r = DoubleDouble::full_add_f64(-ZERO, x);
462    r = DoubleDouble::from_exact_add(r.hi, r.lo);
463
464    let p0 = f_polyeval24(
465        r.to_f64(),
466        f64::from_bits(C[4].1),
467        f64::from_bits(C[5].1),
468        f64::from_bits(C[6].1),
469        f64::from_bits(C[7].1),
470        f64::from_bits(C[8].1),
471        f64::from_bits(C[9].1),
472        f64::from_bits(C[10].1),
473        f64::from_bits(C[11].1),
474        f64::from_bits(C[12].1),
475        f64::from_bits(C[13].1),
476        f64::from_bits(C[14].1),
477        f64::from_bits(C[15].1),
478        f64::from_bits(C[16].1),
479        f64::from_bits(C[17].1),
480        f64::from_bits(C[18].1),
481        f64::from_bits(C[19].1),
482        f64::from_bits(C[20].1),
483        f64::from_bits(C[21].1),
484        f64::from_bits(C[22].1),
485        f64::from_bits(C[23].1),
486        f64::from_bits(C[24].1),
487        f64::from_bits(C[25].1),
488        f64::from_bits(C[26].1),
489        f64::from_bits(C[27].1),
490    );
491
492    let mut p = DoubleDouble::mul_f64_add(r, p0, DoubleDouble::from_bit_pair(C[3]));
493    p = DoubleDouble::mul_add(p, r, DoubleDouble::from_bit_pair(C[2]));
494    p = DoubleDouble::mul_add(p, r, DoubleDouble::from_bit_pair(C[1]));
495    p = DoubleDouble::mul_add(p, r, DoubleDouble::from_bit_pair(C[0]));
496
497    let err = f_fmla(
498        p.hi,
499        f64::from_bits(0x3c50000000000000), // 2^-58
500        f64::from_bits(0x3b90000000000000), // 2^-70
501    );
502    let ub = p.hi + (p.lo + err);
503    let lb = p.hi + (p.lo - err);
504    if ub == lb {
505        return p.to_f64();
506    }
507    y1_transient_zone(x)
508}
509
510fn y1_transient_zone(x: f64) -> f64 {
511    /*
512    <<FunctionApproximations`
513    ClearAll["Global`*"]
514    f[x_]:= BesselY[1,x + 2.1971413260310170351490335626990]
515    {approx,error} =MiniMaxApproximation[f[x],{x,{1.42 -  2.1971413260310170351490335626990, 2-  2.1971413260310170351490335626990 },27,0},WorkingPrecision->120]
516    poly=error[[1]];
517    coeffs=CoefficientList[poly,x];
518    TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
519     */
520    const C: [(u64, u64); 28] = [
521        (0x38d23216c8eac2ee, 0x3c631de77e55e9b0),
522        (0x3c77f3a6718f8e03, 0x3fe0aa48442f0150),
523        (0xbc52b77d09be8679, 0xbfbe56f82217b33a),
524        (0x3c237af39fc9d759, 0xbfa0d2af4e922afe),
525        (0xbc08241e95389bc3, 0xbf73a6dec2f9f739),
526        (0x3c18eac6a35acd63, 0x3f7e671c82a1c956),
527        (0x3be0a0f3c8908083, 0xbf65429d5dbc3bb0),
528        (0x3be217fa58861600, 0x3f517abad71c26c0),
529        (0x3bebc327d6c65514, 0xbf40b28b4ef33a56),
530        (0xbbcb597d6bdd5992, 0x3f2ef0d150ec9934),
531        (0x3bb1a480de696e07, 0xbf1c0758a86844be),
532        (0xbb9c6ea6352ab84e, 0x3f0b7c88b58d9ef3),
533        (0x3b7dbb78dfd868a9, 0xbee97ef9519bdcfd),
534        (0xbb8d4519030f499d, 0x3f04815f08e7ad5e),
535        (0x3bbbd70e1480e260, 0x3f10f348483b57bc),
536        (0x3b7117cf4d2b6f3c, 0x3f231914389bb1bb),
537        (0x3b8ca48beaf6a58d, 0x3f30b29e838345b4),
538        (0x3bccfac65ce17cf9, 0x3f39a69e98f61897),
539        (0xbbeae7e3065b09c9, 0x3f40b9511666fcf0),
540        (0xbbe5cbddf691e7e6, 0x3f428cd8388e634b),
541        (0xbbd91372412d1e1b, 0x3f414ba048d9e1d5),
542        (0xbbb0781a70c6f715, 0x3f3acdcf66f1de95),
543        (0xbba3ae83fd425494, 0x3f30f44ae6620bba),
544        (0x3bc001d75da77b74, 0x3f21154a0a1f2161),
545        (0xbb91c9afb1a1b874, 0x3f0a687b664cbac6),
546        (0x3b8e0a06b9444963, 0x3eed7c3cbb4ba5d8),
547        (0x3b4f0e9dfc915934, 0x3ec53ca23fdd0999),
548        (0xbb0409258f8ffca8, 0x3e8de620acb51b2d),
549    ];
550
551    // this poly relative to first zero
552    const ZERO: DoubleDouble =
553        DoubleDouble::from_bit_pair((0xbc8bd1e50d219bfd, 0x400193bed4dff243));
554
555    let r = DoubleDouble::full_add_f64(-ZERO, x);
556
557    let p0 = f_polyeval13(
558        r.to_f64(),
559        f64::from_bits(C[15].1),
560        f64::from_bits(C[16].1),
561        f64::from_bits(C[17].1),
562        f64::from_bits(C[18].1),
563        f64::from_bits(C[19].1),
564        f64::from_bits(C[20].1),
565        f64::from_bits(C[21].1),
566        f64::from_bits(C[22].1),
567        f64::from_bits(C[23].1),
568        f64::from_bits(C[24].1),
569        f64::from_bits(C[25].1),
570        f64::from_bits(C[26].1),
571        f64::from_bits(C[27].1),
572    );
573
574    let mut p_e = DoubleDouble::mul_f64_add(r, p0, DoubleDouble::from_bit_pair(C[14]));
575    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[13]));
576    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[12]));
577    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[11]));
578    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[10]));
579    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[9]));
580    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[8]));
581    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[7]));
582    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[6]));
583    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[5]));
584    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[4]));
585    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[3]));
586    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[2]));
587    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[1]));
588    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[0]));
589
590    let p = DoubleDouble::from_exact_add(p_e.hi, p_e.lo);
591    let err = f_fmla(
592        p.hi,
593        f64::from_bits(0x3c20000000000000), // 2^-61
594        f64::from_bits(0x3a90000000000000), // 2^-86
595    );
596    let ub = p.hi + (p.lo + err);
597    let lb = p.hi + (p.lo - err);
598    if ub == lb {
599        return p.to_f64();
600    }
601    y1_transient_hard(x)
602}
603
604#[cold]
605#[inline(never)]
606fn y1_transient_hard(x: f64) -> f64 {
607    /*
608    <<FunctionApproximations`
609    ClearAll["Global`*"]
610    f[x_]:= BesselY[1,x + 2.1971413260310170351490335626990]
611    {approx,error} =MiniMaxApproximation[f[x],{x,{1.42 -  2.1971413260310170351490335626990, 2-  2.1971413260310170351490335626990 },27,0},WorkingPrecision->120]
612    poly=error[[1]];
613    coeffs=CoefficientList[poly,x];
614    TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
615     */
616    pub(crate) static C: [DyadicFloat128; 28] = [
617        DyadicFloat128 {
618            sign: DyadicSign::Pos,
619            exponent: -184,
620            mantissa: 0x98ef3bf2_af4d8048_c85b23ab_0bb72488_u128,
621        },
622        DyadicFloat128 {
623            sign: DyadicSign::Pos,
624            exponent: -128,
625            mantissa: 0x85524221_780a817f_3a6718f8_e03513ec_u128,
626        },
627        DyadicFloat128 {
628            sign: DyadicSign::Neg,
629            exponent: -131,
630            mantissa: 0xf2b7c110_bd99d256_efa137d0_cf1f7988_u128,
631        },
632        DyadicFloat128 {
633            sign: DyadicSign::Neg,
634            exponent: -132,
635            mantissa: 0x86957a74_9157ef64_286301b1_453792e2_u128,
636        },
637        DyadicFloat128 {
638            sign: DyadicSign::Neg,
639            exponent: -135,
640            mantissa: 0x9d36f617_cfb9c982_41e95389_bc32e8c2_u128,
641        },
642        DyadicFloat128 {
643            sign: DyadicSign::Pos,
644            exponent: -135,
645            mantissa: 0xf338e415_0e4ab31d_58d46b59_ac6792eb_u128,
646        },
647        DyadicFloat128 {
648            sign: DyadicSign::Neg,
649            exponent: -136,
650            mantissa: 0xaa14eaed_e1dd7f7a_f861bb7b_fbe6acb5_u128,
651        },
652        DyadicFloat128 {
653            sign: DyadicSign::Pos,
654            exponent: -137,
655            mantissa: 0x8bd5d6b8_e1360121_7fa58861_5ff9b614_u128,
656        },
657        DyadicFloat128 {
658            sign: DyadicSign::Neg,
659            exponent: -138,
660            mantissa: 0x85945a77_99d2ac87_9b052735_5d7f5f59_u128,
661        },
662        DyadicFloat128 {
663            sign: DyadicSign::Pos,
664            exponent: -140,
665            mantissa: 0xf7868a87_64c99c94_d0528454_cdccd44c_u128,
666        },
667        DyadicFloat128 {
668            sign: DyadicSign::Neg,
669            exponent: -141,
670            mantissa: 0xe03ac543_4225edcb_6fe432d2_3f11a8b0_u128,
671        },
672        DyadicFloat128 {
673            sign: DyadicSign::Pos,
674            exponent: -142,
675            mantissa: 0xdbe445ac_6cf79639_159cad54_7b1eda7c_u128,
676        },
677        DyadicFloat128 {
678            sign: DyadicSign::Neg,
679            exponent: -144,
680            mantissa: 0xcbf7ca8c_dee7e624_48720279_75757a17_u128,
681        },
682        DyadicFloat128 {
683            sign: DyadicSign::Pos,
684            exponent: -142,
685            mantissa: 0xa40af847_3d6aef15_d737e785_b3163322_u128,
686        },
687        DyadicFloat128 {
688            sign: DyadicSign::Pos,
689            exponent: -141,
690            mantissa: 0x879a4241_dabde37a_e1c2901c_4c06b1dc_u128,
691        },
692        DyadicFloat128 {
693            sign: DyadicSign::Pos,
694            exponent: -140,
695            mantissa: 0x98c8a1c4_dd8dd811_17cf4d2b_6f3c3910_u128,
696        },
697        DyadicFloat128 {
698            sign: DyadicSign::Pos,
699            exponent: -139,
700            mantissa: 0x8594f41c_1a2da01c_a48beaf6_a58cef9f_u128,
701        },
702        DyadicFloat128 {
703            sign: DyadicSign::Pos,
704            exponent: -139,
705            mantissa: 0xcd34f4c7_b0c4b9cf_ac65ce17_cf940c9c_u128,
706        },
707        DyadicFloat128 {
708            sign: DyadicSign::Pos,
709            exponent: -138,
710            mantissa: 0x85ca88b3_37e77ca3_039f349e_c6ddb06d_u128,
711        },
712        DyadicFloat128 {
713            sign: DyadicSign::Pos,
714            exponent: -138,
715            mantissa: 0x9466c1c4_731a5546_84412dc3_033a8ee7_u128,
716        },
717        DyadicFloat128 {
718            sign: DyadicSign::Pos,
719            exponent: -138,
720            mantissa: 0x8a5d0246_cf0ea66e_c8dbed2e_1e50b14f_u128,
721        },
722        DyadicFloat128 {
723            sign: DyadicSign::Pos,
724            exponent: -139,
725            mantissa: 0xd66e7b37_8ef4a77c_3f2c79c8_4757a40d_u128,
726        },
727        DyadicFloat128 {
728            sign: DyadicSign::Pos,
729            exponent: -139,
730            mantissa: 0x87a25733_105dcfb1_45f00af6_adb11b5f_u128,
731        },
732        DyadicFloat128 {
733            sign: DyadicSign::Pos,
734            exponent: -140,
735            mantissa: 0x88aa5050_f90b0a00_3aebb4ef_6e7e9ce1_u128,
736        },
737        DyadicFloat128 {
738            sign: DyadicSign::Pos,
739            exponent: -142,
740            mantissa: 0xd343db32_65d62ee3_6504e5e4_78c55965_u128,
741        },
742        DyadicFloat128 {
743            sign: DyadicSign::Pos,
744            exponent: -144,
745            mantissa: 0xebe1e5da_5d2ec3c1_40d72889_2c57e483_u128,
746        },
747        DyadicFloat128 {
748            sign: DyadicSign::Pos,
749            exponent: -146,
750            mantissa: 0xa9e511fe_e84cc8f8_74efe48a_c9a09ebc_u128,
751        },
752        DyadicFloat128 {
753            sign: DyadicSign::Pos,
754            exponent: -150,
755            mantissa: 0xef310565_a8d9675f_b6d38380_1abf92a6_u128,
756        },
757    ];
758    const ZERO: DyadicFloat128 = DyadicFloat128 {
759        sign: DyadicSign::Pos,
760        exponent: -126,
761        mantissa: 0x8c9df6a6_ff921721_70d796f3_2017e155_u128,
762    };
763    let r = DyadicFloat128::new_from_f64(x) - ZERO;
764
765    let mut p = C[27];
766    for i in (0..27).rev() {
767        p = r * p + C[i];
768    }
769    p.fast_as_f64()
770}
771
772/// This method on small range searches for nearest zero or extremum.
773/// Then picks stored series expansion at the point end evaluates the poly at the point.
774#[inline]
775pub(crate) fn y1_small_argument_fast(x: f64) -> f64 {
776    // let avg_step = 51.03 / 33.0;
777    // let inv_step = 1.0 / avg_step;
778    //
779    // println!("inv_step {}", inv_step);
780
781    const INV_STEP: f64 = 0.6466784244562023;
782
783    let fx = x * INV_STEP;
784    const Y1_ZEROS_COUNT: f64 = (Y1_ZEROS.len() - 1) as f64;
785    let idx0 = unsafe { fx.min(Y1_ZEROS_COUNT).to_int_unchecked::<usize>() };
786    let idx1 = unsafe { fx.ceil().min(Y1_ZEROS_COUNT).to_int_unchecked::<usize>() };
787
788    let found_zero0 = DoubleDouble::from_bit_pair(Y1_ZEROS[idx0]);
789    let found_zero1 = DoubleDouble::from_bit_pair(Y1_ZEROS[idx1]);
790
791    let dist0 = (found_zero0.hi - x).abs();
792    let dist1 = (found_zero1.hi - x).abs();
793
794    let (found_zero, idx, dist) = if dist0 < dist1 {
795        (found_zero0, idx0, dist0)
796    } else {
797        (found_zero1, idx1, dist1)
798    };
799
800    if idx == 0 {
801        return y1_near_zero_fast(x);
802    }
803
804    let close_to_zero = dist.abs() < 1e-3;
805
806    let c = if close_to_zero {
807        &Y1_COEFFS[idx - 1]
808    } else {
809        &Y1_COEFFS_REMEZ[idx - 1]
810    };
811
812    let r = DoubleDouble::full_add_f64(-found_zero, x);
813
814    // We hit exact zero, value, better to return it directly
815    if dist == 0. {
816        return f64::from_bits(Y1_ZEROS_VALUES[idx]);
817    }
818
819    let p = f_polyeval22(
820        r.hi,
821        f64::from_bits(c[6].1),
822        f64::from_bits(c[7].1),
823        f64::from_bits(c[8].1),
824        f64::from_bits(c[9].1),
825        f64::from_bits(c[10].1),
826        f64::from_bits(c[11].1),
827        f64::from_bits(c[12].1),
828        f64::from_bits(c[13].1),
829        f64::from_bits(c[14].1),
830        f64::from_bits(c[15].1),
831        f64::from_bits(c[16].1),
832        f64::from_bits(c[17].1),
833        f64::from_bits(c[18].1),
834        f64::from_bits(c[19].1),
835        f64::from_bits(c[20].1),
836        f64::from_bits(c[21].1),
837        f64::from_bits(c[22].1),
838        f64::from_bits(c[23].1),
839        f64::from_bits(c[24].1),
840        f64::from_bits(c[25].1),
841        f64::from_bits(c[26].1),
842        f64::from_bits(c[27].1),
843    );
844
845    let mut z = DoubleDouble::mul_f64_add(r, p, DoubleDouble::from_bit_pair(c[5]));
846    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[4]));
847    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[3]));
848    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[2]));
849    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[1]));
850    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[0]));
851    let p = z;
852    let err = f_fmla(
853        p.hi,
854        f64::from_bits(0x3c60000000000000), // 2^-57
855        f64::from_bits(0x3c20000000000000), // 2^-61
856    );
857    let ub = p.hi + (p.lo + err);
858    let lb = p.hi + (p.lo - err);
859    if ub == lb {
860        return p.to_f64();
861    }
862    y0_small_argument_moderate(r, c)
863}
864
865fn y0_small_argument_moderate(r: DoubleDouble, c0: &[(u64, u64); 28]) -> f64 {
866    let c = &c0[15..];
867
868    let p0 = f_polyeval13(
869        r.to_f64(),
870        f64::from_bits(c[0].1),
871        f64::from_bits(c[1].1),
872        f64::from_bits(c[2].1),
873        f64::from_bits(c[3].1),
874        f64::from_bits(c[4].1),
875        f64::from_bits(c[5].1),
876        f64::from_bits(c[6].1),
877        f64::from_bits(c[7].1),
878        f64::from_bits(c[8].1),
879        f64::from_bits(c[9].1),
880        f64::from_bits(c[10].1),
881        f64::from_bits(c[11].1),
882        f64::from_bits(c[12].1),
883    );
884
885    let c = c0;
886
887    let mut p_e = DoubleDouble::mul_f64_add(r, p0, DoubleDouble::from_bit_pair(c[14]));
888    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[13]));
889    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[12]));
890    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[11]));
891    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[10]));
892    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[9]));
893    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[8]));
894    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[7]));
895    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[6]));
896    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[5]));
897    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[4]));
898    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[3]));
899    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[2]));
900    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[1]));
901    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[0]));
902
903    let p = DoubleDouble::from_exact_add(p_e.hi, p_e.lo);
904    let err = f_fmla(
905        p.hi,
906        f64::from_bits(0x3c30000000000000), // 2^-60
907        f64::from_bits(0x3bf0000000000000), // 2^-64
908    );
909    let ub = p.hi + (p.lo + err);
910    let lb = p.hi + (p.lo - err);
911    if ub == lb {
912        return p.to_f64();
913    }
914    y1_small_argument_hard(r, c)
915}
916
917#[cold]
918#[inline(never)]
919fn y1_small_argument_hard(r: DoubleDouble, c: &[(u64, u64); 28]) -> f64 {
920    // if we're too close to zero taylor will converge faster and more accurate,
921    // since remez, minimax and other almost cannot polynomial optimize near zero
922    let mut p = DoubleDouble::from_bit_pair(c[27]);
923    for i in (0..27).rev() {
924        p = DoubleDouble::mul_add(r, p, DoubleDouble::from_bit_pair(c[i]));
925        p = DoubleDouble::from_exact_add(p.hi, p.lo);
926    }
927    p.to_f64()
928}
929
930/*
931   Evaluates:
932   Y1 = sqrt(2/(PI*x)) * beta(x) * sin(x - 3*PI/4 - alpha(x))
933
934   Discarding 1/2*PI gives:
935   Y1 = sqrt(2/(PI*x)) * beta(x) * (-cos(x - PI/4 - alpha(x)))
936*/
937#[inline]
938pub(crate) fn y1_asympt_fast(x: f64) -> f64 {
939    const SQRT_2_OVER_PI: DoubleDouble = DoubleDouble::new(
940        f64::from_bits(0xbc8cbc0d30ebfd15),
941        f64::from_bits(0x3fe9884533d43651),
942    );
943    const MPI_OVER_4: DoubleDouble = DoubleDouble::new(
944        f64::from_bits(0xbc81a62633145c07),
945        f64::from_bits(0xbfe921fb54442d18),
946    );
947
948    let recip = if x.to_bits() > 0x7fd000000000000u64 {
949        DoubleDouble::quick_mult_f64(DoubleDouble::from_exact_safe_div(4.0, x), 0.25)
950    } else {
951        DoubleDouble::from_recip(x)
952    };
953
954    let alpha = bessel_1_asympt_alpha_fast(recip);
955    let beta = bessel_1_asympt_beta_fast(recip);
956
957    let AngleReduced { angle } = rem2pi_any(x);
958
959    // Without full subtraction cancellation happens sometimes
960    let x0pi34 = DoubleDouble::full_dd_sub(MPI_OVER_4, alpha);
961    let r0 = DoubleDouble::full_dd_add(angle, x0pi34);
962
963    let m_cos = -cos_dd_small_fast(r0);
964    let z0 = DoubleDouble::quick_mult(beta, m_cos);
965    let r_sqrt = DoubleDouble::from_rsqrt_fast(x);
966    let scale = DoubleDouble::quick_mult(SQRT_2_OVER_PI, r_sqrt);
967    let r = DoubleDouble::quick_mult(scale, z0);
968    let p = DoubleDouble::from_exact_add(r.hi, r.lo);
969    let err = f_fmla(
970        p.hi,
971        f64::from_bits(0x3c40000000000000), // 2^-60
972        f64::from_bits(0x3bf0000000000000), // 2^-87
973    );
974    let ub = p.hi + (p.lo + err);
975    let lb = p.hi + (p.lo - err);
976    if ub == lb {
977        return p.to_f64();
978    }
979    y1_asympt(x, recip, r_sqrt, angle)
980}
981
982/*
983   Evaluates:
984   Y1 = sqrt(2/(PI*x)) * beta(x) * sin(x - 3*PI/4 - alpha(x))
985
986   Discarding 1/2*PI gives:
987   Y1 = sqrt(2/(PI*x)) * beta(x) * (-cos(x - PI/4 - alpha(x)))
988*/
989fn y1_asympt(x: f64, recip: DoubleDouble, r_sqrt: DoubleDouble, angle: DoubleDouble) -> f64 {
990    const SQRT_2_OVER_PI: DoubleDouble = DoubleDouble::new(
991        f64::from_bits(0xbc8cbc0d30ebfd15),
992        f64::from_bits(0x3fe9884533d43651),
993    );
994    const MPI_OVER_4: DoubleDouble = DoubleDouble::new(
995        f64::from_bits(0xbc81a62633145c07),
996        f64::from_bits(0xbfe921fb54442d18),
997    );
998
999    let alpha = bessel_1_asympt_alpha(recip);
1000    let beta = bessel_1_asympt_beta(recip);
1001
1002    // Without full subtraction cancellation happens sometimes
1003    let x0pi34 = DoubleDouble::full_dd_sub(MPI_OVER_4, alpha);
1004    let r0 = DoubleDouble::full_dd_add(angle, x0pi34);
1005
1006    let m_cos = -cos_dd_small(r0);
1007    let z0 = DoubleDouble::quick_mult(beta, m_cos);
1008    let scale = DoubleDouble::quick_mult(SQRT_2_OVER_PI, r_sqrt);
1009    let r = DoubleDouble::quick_mult(scale, z0);
1010    let p = DoubleDouble::from_exact_add(r.hi, r.lo);
1011    let err = f_fmla(
1012        p.hi,
1013        f64::from_bits(0x3c30000000000000), // 2^-60
1014        f64::from_bits(0x3a80000000000000), // 2^-87
1015    );
1016    let ub = p.hi + (p.lo + err);
1017    let lb = p.hi + (p.lo - err);
1018    if ub == lb {
1019        return p.to_f64();
1020    }
1021    y1_asympt_hard(x)
1022}
1023
1024/*
1025   Evaluates:
1026   Y1 = sqrt(2/(PI*x)) * beta(x) * sin(x - 3*PI/4 - alpha(x))
1027
1028   Discarding 1/2*PI gives:
1029   Y1 = sqrt(2/(PI*x)) * beta(x) * (-cos(x - PI/4 - alpha(x)))
1030*/
1031#[cold]
1032#[inline(never)]
1033fn y1_asympt_hard(x: f64) -> f64 {
1034    const SQRT_2_OVER_PI: DyadicFloat128 = DyadicFloat128 {
1035        sign: DyadicSign::Pos,
1036        exponent: -128,
1037        mantissa: 0xcc42299e_a1b28468_7e59e280_5d5c7180_u128,
1038    };
1039
1040    const MPI_OVER_4: DyadicFloat128 = DyadicFloat128 {
1041        sign: DyadicSign::Neg,
1042        exponent: -128,
1043        mantissa: 0xc90fdaa2_2168c234_c4c6628b_80dc1cd1_u128,
1044    };
1045
1046    let x_dyadic = DyadicFloat128::new_from_f64(x);
1047    let recip = DyadicFloat128::accurate_reciprocal(x);
1048
1049    let alpha = bessel_1_asympt_alpha_hard(recip);
1050    let beta = bessel_1_asympt_beta_hard(recip);
1051
1052    let angle = rem2pi_f128(x_dyadic);
1053
1054    let x0pi34 = MPI_OVER_4 - alpha;
1055    let r0 = angle + x0pi34;
1056
1057    let m_cos = cos_f128_small(r0).negated();
1058
1059    let z0 = beta * m_cos;
1060    let r_sqrt = bessel_rsqrt_hard(x, recip);
1061    let scale = SQRT_2_OVER_PI * r_sqrt;
1062    let p = scale * z0;
1063    p.fast_as_f64()
1064}
1065
1066#[cfg(test)]
1067mod tests {
1068    use super::*;
1069
1070    #[test]
1071    fn test_y1() {
1072        // ULP should be less than 0.500001, but it was 0.5089558379720955, on 2.1957931471395398 result -0.0007023285780874727, using f_y1 and MPFR -0.0007023285780874729
1073        assert_eq!(f_y1(0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000007291282546733975),
1074                   -873124540555277200000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.);
1075        assert_eq!(f_y1(2.1957931471395398), -0.0007023285780874727);
1076        assert_eq!(
1077            f_y1(f64::from_bits(0x571a31ffe2ff7e9fu64)),
1078            f64::from_bits(0x32e58532f95056ffu64)
1079        );
1080        assert_eq!(
1081            f_y1(f64::from_bits(0x400193bed4dff243)),
1082            0.00000000000000002513306678922122
1083        );
1084        assert_eq!(
1085            f_y1(f64::from_bits(0x3ffc513c569fe78e)),
1086            -0.24189760895998239
1087        );
1088        assert_eq!(
1089            f_y1(f64::from_bits(0x4192391e4c8faa60)),
1090            -0.000000000000000002572292246748134
1091        );
1092        assert_eq!(
1093            f_y1(f64::from_bits(0x403e9e480605283c)),
1094            -0.00000000000000001524456280251315
1095        );
1096        assert_eq!(
1097            f_y1(f64::from_bits(0x40277f9138d43206)),
1098            0.000000000000000006849807120770496
1099        );
1100        assert_eq!(f_y1(f64::INFINITY), 0.);
1101        assert!(f_y1(f64::NEG_INFINITY).is_nan());
1102        assert!(f_y1(f64::NAN).is_nan());
1103    }
1104
1105    #[test]
1106    fn test_y1_edge_cases() {
1107        assert_eq!(f_y1(2.1904620854463985), -0.0034837351616785234);
1108        assert_eq!(f_y1(2.197142201034536), 4.5568985277260593e-7);
1109        assert_eq!(f_y1(1.4000000000000004), -0.4791469742327998);
1110        assert_eq!(f_y1(2.0002288794493848), -0.1069033735586767);
1111    }
1112}
pxfm/bessel/y1.rs

pxfm/bessel/
y1.rs
