pxfm 0.1.29

/*
 * // Copyright (c) Radzivon Bartoshyk 7/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
use crate::bessel::alpha1::{
    bessel_1_asympt_alpha, bessel_1_asympt_alpha_fast, bessel_1_asympt_alpha_hard,
};
use crate::bessel::beta1::{
    bessel_1_asympt_beta, bessel_1_asympt_beta_fast, bessel_1_asympt_beta_hard,
};
use crate::bessel::i0::bessel_rsqrt_hard;
use crate::bessel::y1_coeffs::Y1_COEFFS_REMEZ;
use crate::bessel::y1_coeffs_taylor::Y1_COEFFS;
use crate::bessel::y1f_coeffs::{Y1_ZEROS, Y1_ZEROS_VALUES};
use crate::common::f_fmla;
use crate::double_double::DoubleDouble;
use crate::dyadic_float::{DyadicFloat128, DyadicSign};
use crate::logs::log_dd_fast;
use crate::polyeval::{f_polyeval12, f_polyeval13, f_polyeval15, f_polyeval22, f_polyeval24};
use crate::rounding::CpuCeil;
use crate::sin_helper::{cos_dd_small, cos_dd_small_fast, cos_f128_small};
use crate::sincos_reduce::{AngleReduced, rem2pi_any, rem2pi_f128};

/// Bessel of the second kind of order 1 ( Y1 )
pub fn f_y1(x: f64) -> f64 {
    let ix = x.to_bits();

    if ix >= 0x7ffu64 << 52 || ix == 0 {
        // |x| == NaN, x == inf, |x| == 0, x < 0
        if ix.wrapping_shl(1) == 0 {
            // |x| == 0
            return f64::NEG_INFINITY;
        }

        if x.is_infinite() {
            if x.is_sign_negative() {
                return f64::NAN;
            }
            return 0.;
        }

        return x + f64::NAN; // x == NaN
    }

    let xb = x.to_bits();

    if xb <= 0x4049c00000000000u64 {
        // x <= 51.5
        if xb <= 0x4000000000000000u64 {
            // x <= 2
            if xb <= 0x3ff75c28f5c28f5cu64 {
                // x <= 1.46
                return y1_near_zero_fast(x);
            }
            // transient zone from 1.46 to 2 have bad behavior for log poly already,
            // and not yet good to be easily covered, thus it use its own poly
            return y1_transient_zone_fast(x);
        }

        return y1_small_argument_fast(x);
    }

    y1_asympt_fast(x)
}

/**
Generated by SageMath:
Evaluates:
y2 = -J1(x)*log(x) + 1/x * (1 - sum((-1)^m*(H(m)+H(m-1))/(2^m*m!*(m-1)!)*x^(2*m))
Y1(x) = 2/pi*(-y2(x)+(euler_gamma - log(2))*J1(x))
expressed as:
Y1(x)=log(x)*W1(x) - Z1(x) - 2/(pi*x)
```python
from sage.all import *

R = LaurentSeriesRing(RealField(300), 'x', default_prec=300)
x = R.gen()
N = 16  # Number of terms (adjust as needed)
gamma = RealField(300)(euler_gamma)
d2 = RealField(300)(2)
pi = RealField(300).pi()
log2 = RealField(300)(2).log()

def j_series(n, x):
    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)])

J1_series = j_series(1, x)

def harmony(m):
    return sum(RealField(300)(1)/RealField(300)(k) for k in range(1, m+1))

def z_series(x):
    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)])

W1 = d2/pi * J1_series
Z1 = -(d2/(x*pi) * z_series(x) + d2/pi * gamma * J1_series(x) - d2/pi * log2 * J1_series(x))

def y1_full(x):
    return d2/pi * (J1_series(x) * x.log() - 1/x * ( 1 - z_series(x)) + (gamma - log2) * J1_series(x))

# see the series
print(W0)
print(Z0)
```
See ./notes/bessel_y1_taylor.ipynb for generation
**/
#[inline]
fn y1_near_zero_fast(x: f64) -> f64 {
    const W: [(u64, u64); 15] = [
        (0xbc76b01ec5417056, 0x3fd45f306dc9c883),
        (0x3c46b01ec5417056, 0xbfa45f306dc9c883),
        (0xbbfe40290701eb1e, 0x3f5b2995e7b7b604),
        (0xbba67fe4a5feb897, 0xbf021bb945252402),
        (0xbb0334914cdd2ba9, 0x3e9cf9286ea1d337),
        (0x3aab4f3c6d42c1f4, 0xbe2ee7a29824147f),
        (0xba407fb57ef4dc2c, 0x3db78be9987d036d),
        (0x39d2921e91b07dd0, 0xbd3ae90af76a4d0f),
        (0x395a28c8620dc90e, 0x3cb7eb97f85e7d62),
        (0xb8cf83f52abe45c5, 0xbc31028e3376648a),
        (0xb8441050c68ca435, 0x3ba3cb1e7d0c17e7),
        (0xb7ab072548a1aa43, 0xbb133191ed9f1eef),
        (0xb6f05192c2d9b6ee, 0x3a7f7f4b5e8ef7b0),
        (0x367ad65afe306d57, 0xb9e626e36cb3515d),
        (0xb5ea1c4136f8f230, 0x394b01153dce6810),
    ];
    let x2 = DoubleDouble::from_exact_mult(x, x);
    let w0 = f_polyeval12(
        x2.hi,
        f64::from_bits(W[3].1),
        f64::from_bits(W[4].1),
        f64::from_bits(W[5].1),
        f64::from_bits(W[6].1),
        f64::from_bits(W[7].1),
        f64::from_bits(W[8].1),
        f64::from_bits(W[9].1),
        f64::from_bits(W[10].1),
        f64::from_bits(W[11].1),
        f64::from_bits(W[12].1),
        f64::from_bits(W[13].1),
        f64::from_bits(W[14].1),
    );

    let mut w = DoubleDouble::mul_f64_add(x2, w0, DoubleDouble::from_bit_pair(W[2]));
    w = DoubleDouble::mul_add(x2, w, DoubleDouble::from_bit_pair(W[1]));
    w = DoubleDouble::mul_add(x2, w, DoubleDouble::from_bit_pair(W[0]));
    w = DoubleDouble::quick_mult_f64(w, x);

    const Z: [(u64, u64); 15] = [
        (0x3c61d7eb2e54cda1, 0x3fc91866143cbc8a),
        (0xbc2f9f7a0ce54a40, 0xbfabd3975c75b4a7),
        (0xbbf7659313f45e8c, 0x3f6835b97894be5b),
        (0x3b9cbcd40f1be7b9, 0xbf12c7dbffcde97d),
        (0xbb495d78778645b4, 0x3eb0a780ac776eac),
        (0xbae15be86455c1ab, 0xbe432e5a4ddeea30),
        (0xba5ad966c12f1e3c, 0x3dcf0ce34d2066a6),
        (0x39e9717155dc7521, 0xbd52a4e1aea45c18),
        (0x394f447fe5de1290, 0x3cd1474ade9154ac),
        (0xb8e1699d9009a7fc, 0xbc4978ba84f218c0),
        (0xb8505502096ead17, 0x3bbe9598c016378b),
        (0x37942b6c36b2c5f1, 0xbb2e7e5fcfc4b7b1),
        (0x37210853b78bd08a, 0x3a99a6c1266c116d),
        (0xb686c9639c9d976e, 0xba02738998fe7337),
        (0xb603b739ee04b9fe, 0x3966f58cd41b6d08),
    ];
    let z0 = f_polyeval12(
        x2.hi,
        f64::from_bits(Z[3].1),
        f64::from_bits(Z[4].1),
        f64::from_bits(Z[5].1),
        f64::from_bits(Z[6].1),
        f64::from_bits(Z[7].1),
        f64::from_bits(Z[8].1),
        f64::from_bits(Z[9].1),
        f64::from_bits(Z[10].1),
        f64::from_bits(Z[11].1),
        f64::from_bits(Z[12].1),
        f64::from_bits(Z[13].1),
        f64::from_bits(Z[14].1),
    );

    let mut z = DoubleDouble::mul_f64_add(x2, z0, DoubleDouble::from_bit_pair(Z[2]));
    z = DoubleDouble::mul_add(x2, z, DoubleDouble::from_bit_pair(Z[1]));
    z = DoubleDouble::mul_add(x2, z, DoubleDouble::from_bit_pair(Z[0]));
    z = DoubleDouble::quick_mult_f64(z, x);

    let w_log = log_dd_fast(x);

    const MINUS_TWO_OVER_PI: DoubleDouble =
        DoubleDouble::from_bit_pair((0x3c86b01ec5417056, 0xbfe45f306dc9c883));

    let m_two_over_pi_div_x: DoubleDouble;
    #[cfg(any(
        all(
            any(target_arch = "x86", target_arch = "x86_64"),
            target_feature = "fma"
        ),
        target_arch = "aarch64"
    ))]
    {
        m_two_over_pi_div_x = DoubleDouble::div_dd_f64(MINUS_TWO_OVER_PI, x)
    }
    #[cfg(not(any(
        all(
            any(target_arch = "x86", target_arch = "x86_64"),
            target_feature = "fma"
        ),
        target_arch = "aarch64"
    )))]
    {
        use crate::double_double::two_product_compatible;
        m_two_over_pi_div_x = if two_product_compatible(x) {
            DoubleDouble::div_dd_f64(MINUS_TWO_OVER_PI, x)
        } else {
            DoubleDouble::div_safe_dd_f64(MINUS_TWO_OVER_PI, x)
        };
    }
    if m_two_over_pi_div_x.hi.is_infinite() {
        return f64::NEG_INFINITY;
    }

    let zvp = DoubleDouble::mul_add(w, w_log, -z);
    let p = DoubleDouble::add(m_two_over_pi_div_x, zvp);
    let err = f_fmla(
        p.hi,
        f64::from_bits(0x3c30000000000000), // 2^-60
        f64::from_bits(0x3be0000000000000), // 2^-65
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);
    if ub == lb {
        return p.to_f64();
    }
    y1_near_zero(x, w_log)
}

/**
Generated by SageMath:
Evaluates:
y2 = -J1(x)*log(x) + 1/x * (1 - sum((-1)^m*(H(m)+H(m-1))/(2^m*m!*(m-1)!)*x^(2*m))
Y1(x) = 2/pi*(-y2(x)+(euler_gamma - log(2))*J1(x))
expressed as:
Y1(x)=log(x)*W1(x) - Z1(x) - 2/(pi*x)
```python
from sage.all import *

R = LaurentSeriesRing(RealField(300), 'x', default_prec=300)
x = R.gen()
N = 16  # Number of terms (adjust as needed)
gamma = RealField(300)(euler_gamma)
d2 = RealField(300)(2)
pi = RealField(300).pi()
log2 = RealField(300)(2).log()

def j_series(n, x):
    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)])

J1_series = j_series(1, x)

def harmony(m):
    return sum(RealField(300)(1)/RealField(300)(k) for k in range(1, m+1))

def z_series(x):
    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)])

W1 = d2/pi * J1_series
Z1 = -(d2/(x*pi) * z_series(x) + d2/pi * gamma * J1_series(x) - d2/pi * log2 * J1_series(x))

def y1_full(x):
    return d2/pi * (J1_series(x) * x.log() - 1/x * ( 1 - z_series(x)) + (gamma - log2) * J1_series(x))

# see the series
print(W0)
print(Z0)
```
See ./notes/bessel_y1_taylor.ipynb for generation
**/
#[cold]
#[inline(never)]
fn y1_near_zero(x: f64, w_log: DoubleDouble) -> f64 {
    const W: [(u64, u64); 15] = [
        (0xbc76b01ec5417056, 0x3fd45f306dc9c883),
        (0x3c46b01ec5417056, 0xbfa45f306dc9c883),
        (0xbbfe40290701eb1e, 0x3f5b2995e7b7b604),
        (0xbba67fe4a5feb897, 0xbf021bb945252402),
        (0xbb0334914cdd2ba9, 0x3e9cf9286ea1d337),
        (0x3aab4f3c6d42c1f4, 0xbe2ee7a29824147f),
        (0xba407fb57ef4dc2c, 0x3db78be9987d036d),
        (0x39d2921e91b07dd0, 0xbd3ae90af76a4d0f),
        (0x395a28c8620dc90e, 0x3cb7eb97f85e7d62),
        (0xb8cf83f52abe45c5, 0xbc31028e3376648a),
        (0xb8441050c68ca435, 0x3ba3cb1e7d0c17e7),
        (0xb7ab072548a1aa43, 0xbb133191ed9f1eef),
        (0xb6f05192c2d9b6ee, 0x3a7f7f4b5e8ef7b0),
        (0x367ad65afe306d57, 0xb9e626e36cb3515d),
        (0xb5ea1c4136f8f230, 0x394b01153dce6810),
    ];
    let x2 = DoubleDouble::from_exact_mult(x, x);
    let mut w = f_polyeval15(
        x2,
        DoubleDouble::from_bit_pair(W[0]),
        DoubleDouble::from_bit_pair(W[1]),
        DoubleDouble::from_bit_pair(W[2]),
        DoubleDouble::from_bit_pair(W[3]),
        DoubleDouble::from_bit_pair(W[4]),
        DoubleDouble::from_bit_pair(W[5]),
        DoubleDouble::from_bit_pair(W[6]),
        DoubleDouble::from_bit_pair(W[7]),
        DoubleDouble::from_bit_pair(W[8]),
        DoubleDouble::from_bit_pair(W[9]),
        DoubleDouble::from_bit_pair(W[10]),
        DoubleDouble::from_bit_pair(W[11]),
        DoubleDouble::from_bit_pair(W[12]),
        DoubleDouble::from_bit_pair(W[13]),
        DoubleDouble::from_bit_pair(W[14]),
    );
    w = DoubleDouble::quick_mult_f64(w, x);

    const Z: [(u64, u64); 15] = [
        (0x3c61d7eb2e54cda1, 0x3fc91866143cbc8a),
        (0xbc2f9f7a0ce54a40, 0xbfabd3975c75b4a7),
        (0xbbf7659313f45e8c, 0x3f6835b97894be5b),
        (0x3b9cbcd40f1be7b9, 0xbf12c7dbffcde97d),
        (0xbb495d78778645b4, 0x3eb0a780ac776eac),
        (0xbae15be86455c1ab, 0xbe432e5a4ddeea30),
        (0xba5ad966c12f1e3c, 0x3dcf0ce34d2066a6),
        (0x39e9717155dc7521, 0xbd52a4e1aea45c18),
        (0x394f447fe5de1290, 0x3cd1474ade9154ac),
        (0xb8e1699d9009a7fc, 0xbc4978ba84f218c0),
        (0xb8505502096ead17, 0x3bbe9598c016378b),
        (0x37942b6c36b2c5f1, 0xbb2e7e5fcfc4b7b1),
        (0x37210853b78bd08a, 0x3a99a6c1266c116d),
        (0xb686c9639c9d976e, 0xba02738998fe7337),
        (0xb603b739ee04b9fe, 0x3966f58cd41b6d08),
    ];
    let mut z = f_polyeval15(
        x2,
        DoubleDouble::from_bit_pair(Z[0]),
        DoubleDouble::from_bit_pair(Z[1]),
        DoubleDouble::from_bit_pair(Z[2]),
        DoubleDouble::from_bit_pair(Z[3]),
        DoubleDouble::from_bit_pair(Z[4]),
        DoubleDouble::from_bit_pair(Z[5]),
        DoubleDouble::from_bit_pair(Z[6]),
        DoubleDouble::from_bit_pair(Z[7]),
        DoubleDouble::from_bit_pair(Z[8]),
        DoubleDouble::from_bit_pair(Z[9]),
        DoubleDouble::from_bit_pair(Z[10]),
        DoubleDouble::from_bit_pair(Z[11]),
        DoubleDouble::from_bit_pair(Z[12]),
        DoubleDouble::from_bit_pair(Z[13]),
        DoubleDouble::from_bit_pair(Z[14]),
    );
    z = DoubleDouble::quick_mult_f64(z, x);

    const MINUS_TWO_OVER_PI: DoubleDouble =
        DoubleDouble::from_bit_pair((0x3c86b01ec5417056, 0xbfe45f306dc9c883));

    let m_two_over_pi_div_x: DoubleDouble;
    #[cfg(any(
        all(
            any(target_arch = "x86", target_arch = "x86_64"),
            target_feature = "fma"
        ),
        target_arch = "aarch64"
    ))]
    {
        m_two_over_pi_div_x = DoubleDouble::div_dd_f64(MINUS_TWO_OVER_PI, x)
    }
    #[cfg(not(any(
        all(
            any(target_arch = "x86", target_arch = "x86_64"),
            target_feature = "fma"
        ),
        target_arch = "aarch64"
    )))]
    {
        use crate::double_double::two_product_compatible;
        m_two_over_pi_div_x = if two_product_compatible(x) {
            DoubleDouble::div_dd_f64(MINUS_TWO_OVER_PI, x)
        } else {
            DoubleDouble::div_safe_dd_f64(MINUS_TWO_OVER_PI, x)
        };
    }
    if m_two_over_pi_div_x.hi.is_infinite() {
        return f64::NEG_INFINITY;
    }

    let zvp = DoubleDouble::mul_add(w, w_log, -z);
    DoubleDouble::full_dd_add(m_two_over_pi_div_x, zvp).to_f64()
}

#[inline]
fn y1_transient_zone_fast(x: f64) -> f64 {
    /*
    <<FunctionApproximations`
    ClearAll["Global`*"]
    f[x_]:= BesselY[1,x + 2.1971413260310170351490335626990]
    {approx,error} =MiniMaxApproximation[f[x],{x,{1.42 -  2.1971413260310170351490335626990, 2-  2.1971413260310170351490335626990 },27,0},WorkingPrecision->120]
    poly=error[[1]];
    coeffs=CoefficientList[poly,x];
    TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
     */
    const C: [(u64, u64); 28] = [
        (0x38d23216c8eac2ee, 0x3c631de77e55e9b0),
        (0x3c77f3a6718f8e03, 0x3fe0aa48442f0150),
        (0xbc52b77d09be8679, 0xbfbe56f82217b33a),
        (0x3c237af39fc9d759, 0xbfa0d2af4e922afe),
        (0xbc08241e95389bc3, 0xbf73a6dec2f9f739),
        (0x3c18eac6a35acd63, 0x3f7e671c82a1c956),
        (0x3be0a0f3c8908083, 0xbf65429d5dbc3bb0),
        (0x3be217fa58861600, 0x3f517abad71c26c0),
        (0x3bebc327d6c65514, 0xbf40b28b4ef33a56),
        (0xbbcb597d6bdd5992, 0x3f2ef0d150ec9934),
        (0x3bb1a480de696e07, 0xbf1c0758a86844be),
        (0xbb9c6ea6352ab84e, 0x3f0b7c88b58d9ef3),
        (0x3b7dbb78dfd868a9, 0xbee97ef9519bdcfd),
        (0xbb8d4519030f499d, 0x3f04815f08e7ad5e),
        (0x3bbbd70e1480e260, 0x3f10f348483b57bc),
        (0x3b7117cf4d2b6f3c, 0x3f231914389bb1bb),
        (0x3b8ca48beaf6a58d, 0x3f30b29e838345b4),
        (0x3bccfac65ce17cf9, 0x3f39a69e98f61897),
        (0xbbeae7e3065b09c9, 0x3f40b9511666fcf0),
        (0xbbe5cbddf691e7e6, 0x3f428cd8388e634b),
        (0xbbd91372412d1e1b, 0x3f414ba048d9e1d5),
        (0xbbb0781a70c6f715, 0x3f3acdcf66f1de95),
        (0xbba3ae83fd425494, 0x3f30f44ae6620bba),
        (0x3bc001d75da77b74, 0x3f21154a0a1f2161),
        (0xbb91c9afb1a1b874, 0x3f0a687b664cbac6),
        (0x3b8e0a06b9444963, 0x3eed7c3cbb4ba5d8),
        (0x3b4f0e9dfc915934, 0x3ec53ca23fdd0999),
        (0xbb0409258f8ffca8, 0x3e8de620acb51b2d),
    ];

    // this poly relative to first zero
    const ZERO: DoubleDouble =
        DoubleDouble::from_bit_pair((0xbc8bd1e50d219bfd, 0x400193bed4dff243));

    let mut r = DoubleDouble::full_add_f64(-ZERO, x);
    r = DoubleDouble::from_exact_add(r.hi, r.lo);

    let p0 = f_polyeval24(
        r.to_f64(),
        f64::from_bits(C[4].1),
        f64::from_bits(C[5].1),
        f64::from_bits(C[6].1),
        f64::from_bits(C[7].1),
        f64::from_bits(C[8].1),
        f64::from_bits(C[9].1),
        f64::from_bits(C[10].1),
        f64::from_bits(C[11].1),
        f64::from_bits(C[12].1),
        f64::from_bits(C[13].1),
        f64::from_bits(C[14].1),
        f64::from_bits(C[15].1),
        f64::from_bits(C[16].1),
        f64::from_bits(C[17].1),
        f64::from_bits(C[18].1),
        f64::from_bits(C[19].1),
        f64::from_bits(C[20].1),
        f64::from_bits(C[21].1),
        f64::from_bits(C[22].1),
        f64::from_bits(C[23].1),
        f64::from_bits(C[24].1),
        f64::from_bits(C[25].1),
        f64::from_bits(C[26].1),
        f64::from_bits(C[27].1),
    );

    let mut p = DoubleDouble::mul_f64_add(r, p0, DoubleDouble::from_bit_pair(C[3]));
    p = DoubleDouble::mul_add(p, r, DoubleDouble::from_bit_pair(C[2]));
    p = DoubleDouble::mul_add(p, r, DoubleDouble::from_bit_pair(C[1]));
    p = DoubleDouble::mul_add(p, r, DoubleDouble::from_bit_pair(C[0]));

    let err = f_fmla(
        p.hi,
        f64::from_bits(0x3c50000000000000), // 2^-58
        f64::from_bits(0x3b90000000000000), // 2^-70
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);
    if ub == lb {
        return p.to_f64();
    }
    y1_transient_zone(x)
}

fn y1_transient_zone(x: f64) -> f64 {
    /*
    <<FunctionApproximations`
    ClearAll["Global`*"]
    f[x_]:= BesselY[1,x + 2.1971413260310170351490335626990]
    {approx,error} =MiniMaxApproximation[f[x],{x,{1.42 -  2.1971413260310170351490335626990, 2-  2.1971413260310170351490335626990 },27,0},WorkingPrecision->120]
    poly=error[[1]];
    coeffs=CoefficientList[poly,x];
    TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
     */
    const C: [(u64, u64); 28] = [
        (0x38d23216c8eac2ee, 0x3c631de77e55e9b0),
        (0x3c77f3a6718f8e03, 0x3fe0aa48442f0150),
        (0xbc52b77d09be8679, 0xbfbe56f82217b33a),
        (0x3c237af39fc9d759, 0xbfa0d2af4e922afe),
        (0xbc08241e95389bc3, 0xbf73a6dec2f9f739),
        (0x3c18eac6a35acd63, 0x3f7e671c82a1c956),
        (0x3be0a0f3c8908083, 0xbf65429d5dbc3bb0),
        (0x3be217fa58861600, 0x3f517abad71c26c0),
        (0x3bebc327d6c65514, 0xbf40b28b4ef33a56),
        (0xbbcb597d6bdd5992, 0x3f2ef0d150ec9934),
        (0x3bb1a480de696e07, 0xbf1c0758a86844be),
        (0xbb9c6ea6352ab84e, 0x3f0b7c88b58d9ef3),
        (0x3b7dbb78dfd868a9, 0xbee97ef9519bdcfd),
        (0xbb8d4519030f499d, 0x3f04815f08e7ad5e),
        (0x3bbbd70e1480e260, 0x3f10f348483b57bc),
        (0x3b7117cf4d2b6f3c, 0x3f231914389bb1bb),
        (0x3b8ca48beaf6a58d, 0x3f30b29e838345b4),
        (0x3bccfac65ce17cf9, 0x3f39a69e98f61897),
        (0xbbeae7e3065b09c9, 0x3f40b9511666fcf0),
        (0xbbe5cbddf691e7e6, 0x3f428cd8388e634b),
        (0xbbd91372412d1e1b, 0x3f414ba048d9e1d5),
        (0xbbb0781a70c6f715, 0x3f3acdcf66f1de95),
        (0xbba3ae83fd425494, 0x3f30f44ae6620bba),
        (0x3bc001d75da77b74, 0x3f21154a0a1f2161),
        (0xbb91c9afb1a1b874, 0x3f0a687b664cbac6),
        (0x3b8e0a06b9444963, 0x3eed7c3cbb4ba5d8),
        (0x3b4f0e9dfc915934, 0x3ec53ca23fdd0999),
        (0xbb0409258f8ffca8, 0x3e8de620acb51b2d),
    ];

    // this poly relative to first zero
    const ZERO: DoubleDouble =
        DoubleDouble::from_bit_pair((0xbc8bd1e50d219bfd, 0x400193bed4dff243));

    let r = DoubleDouble::full_add_f64(-ZERO, x);

    let p0 = f_polyeval13(
        r.to_f64(),
        f64::from_bits(C[15].1),
        f64::from_bits(C[16].1),
        f64::from_bits(C[17].1),
        f64::from_bits(C[18].1),
        f64::from_bits(C[19].1),
        f64::from_bits(C[20].1),
        f64::from_bits(C[21].1),
        f64::from_bits(C[22].1),
        f64::from_bits(C[23].1),
        f64::from_bits(C[24].1),
        f64::from_bits(C[25].1),
        f64::from_bits(C[26].1),
        f64::from_bits(C[27].1),
    );

    let mut p_e = DoubleDouble::mul_f64_add(r, p0, DoubleDouble::from_bit_pair(C[14]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[13]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[12]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[11]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[10]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[9]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[8]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[7]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[6]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[5]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[4]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[3]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[2]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[1]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(C[0]));

    let p = DoubleDouble::from_exact_add(p_e.hi, p_e.lo);
    let err = f_fmla(
        p.hi,
        f64::from_bits(0x3c20000000000000), // 2^-61
        f64::from_bits(0x3a90000000000000), // 2^-86
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);
    if ub == lb {
        return p.to_f64();
    }
    y1_transient_hard(x)
}

#[cold]
#[inline(never)]
fn y1_transient_hard(x: f64) -> f64 {
    /*
    <<FunctionApproximations`
    ClearAll["Global`*"]
    f[x_]:= BesselY[1,x + 2.1971413260310170351490335626990]
    {approx,error} =MiniMaxApproximation[f[x],{x,{1.42 -  2.1971413260310170351490335626990, 2-  2.1971413260310170351490335626990 },27,0},WorkingPrecision->120]
    poly=error[[1]];
    coeffs=CoefficientList[poly,x];
    TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
     */
    pub(crate) static C: [DyadicFloat128; 28] = [
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -184,
            mantissa: 0x98ef3bf2_af4d8048_c85b23ab_0bb72488_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -128,
            mantissa: 0x85524221_780a817f_3a6718f8_e03513ec_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -131,
            mantissa: 0xf2b7c110_bd99d256_efa137d0_cf1f7988_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -132,
            mantissa: 0x86957a74_9157ef64_286301b1_453792e2_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -135,
            mantissa: 0x9d36f617_cfb9c982_41e95389_bc32e8c2_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -135,
            mantissa: 0xf338e415_0e4ab31d_58d46b59_ac6792eb_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -136,
            mantissa: 0xaa14eaed_e1dd7f7a_f861bb7b_fbe6acb5_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -137,
            mantissa: 0x8bd5d6b8_e1360121_7fa58861_5ff9b614_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -138,
            mantissa: 0x85945a77_99d2ac87_9b052735_5d7f5f59_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -140,
            mantissa: 0xf7868a87_64c99c94_d0528454_cdccd44c_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -141,
            mantissa: 0xe03ac543_4225edcb_6fe432d2_3f11a8b0_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -142,
            mantissa: 0xdbe445ac_6cf79639_159cad54_7b1eda7c_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -144,
            mantissa: 0xcbf7ca8c_dee7e624_48720279_75757a17_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -142,
            mantissa: 0xa40af847_3d6aef15_d737e785_b3163322_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -141,
            mantissa: 0x879a4241_dabde37a_e1c2901c_4c06b1dc_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -140,
            mantissa: 0x98c8a1c4_dd8dd811_17cf4d2b_6f3c3910_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -139,
            mantissa: 0x8594f41c_1a2da01c_a48beaf6_a58cef9f_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -139,
            mantissa: 0xcd34f4c7_b0c4b9cf_ac65ce17_cf940c9c_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -138,
            mantissa: 0x85ca88b3_37e77ca3_039f349e_c6ddb06d_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -138,
            mantissa: 0x9466c1c4_731a5546_84412dc3_033a8ee7_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -138,
            mantissa: 0x8a5d0246_cf0ea66e_c8dbed2e_1e50b14f_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -139,
            mantissa: 0xd66e7b37_8ef4a77c_3f2c79c8_4757a40d_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -139,
            mantissa: 0x87a25733_105dcfb1_45f00af6_adb11b5f_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -140,
            mantissa: 0x88aa5050_f90b0a00_3aebb4ef_6e7e9ce1_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -142,
            mantissa: 0xd343db32_65d62ee3_6504e5e4_78c55965_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -144,
            mantissa: 0xebe1e5da_5d2ec3c1_40d72889_2c57e483_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -146,
            mantissa: 0xa9e511fe_e84cc8f8_74efe48a_c9a09ebc_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -150,
            mantissa: 0xef310565_a8d9675f_b6d38380_1abf92a6_u128,
        },
    ];
    const ZERO: DyadicFloat128 = DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -126,
        mantissa: 0x8c9df6a6_ff921721_70d796f3_2017e155_u128,
    };
    let r = DyadicFloat128::new_from_f64(x) - ZERO;

    let mut p = C[27];
    for i in (0..27).rev() {
        p = r * p + C[i];
    }
    p.fast_as_f64()
}

/// This method on small range searches for nearest zero or extremum.
/// Then picks stored series expansion at the point end evaluates the poly at the point.
#[inline]
pub(crate) fn y1_small_argument_fast(x: f64) -> f64 {
    // let avg_step = 51.03 / 33.0;
    // let inv_step = 1.0 / avg_step;
    //
    // println!("inv_step {}", inv_step);

    const INV_STEP: f64 = 0.6466784244562023;

    let fx = x * INV_STEP;
    const Y1_ZEROS_COUNT: f64 = (Y1_ZEROS.len() - 1) as f64;
    let idx0 = unsafe { fx.min(Y1_ZEROS_COUNT).to_int_unchecked::<usize>() };
    let idx1 = unsafe {
        fx.cpu_ceil()
            .min(Y1_ZEROS_COUNT)
            .to_int_unchecked::<usize>()
    };

    let found_zero0 = DoubleDouble::from_bit_pair(Y1_ZEROS[idx0]);
    let found_zero1 = DoubleDouble::from_bit_pair(Y1_ZEROS[idx1]);

    let dist0 = (found_zero0.hi - x).abs();
    let dist1 = (found_zero1.hi - x).abs();

    let (found_zero, idx, dist) = if dist0 < dist1 {
        (found_zero0, idx0, dist0)
    } else {
        (found_zero1, idx1, dist1)
    };

    if idx == 0 {
        return y1_near_zero_fast(x);
    }

    let close_to_zero = dist.abs() < 1e-3;

    let c = if close_to_zero {
        &Y1_COEFFS[idx - 1]
    } else {
        &Y1_COEFFS_REMEZ[idx - 1]
    };

    let r = DoubleDouble::full_add_f64(-found_zero, x);

    // We hit exact zero, value, better to return it directly
    if dist == 0. {
        return f64::from_bits(Y1_ZEROS_VALUES[idx]);
    }

    let p = f_polyeval22(
        r.hi,
        f64::from_bits(c[6].1),
        f64::from_bits(c[7].1),
        f64::from_bits(c[8].1),
        f64::from_bits(c[9].1),
        f64::from_bits(c[10].1),
        f64::from_bits(c[11].1),
        f64::from_bits(c[12].1),
        f64::from_bits(c[13].1),
        f64::from_bits(c[14].1),
        f64::from_bits(c[15].1),
        f64::from_bits(c[16].1),
        f64::from_bits(c[17].1),
        f64::from_bits(c[18].1),
        f64::from_bits(c[19].1),
        f64::from_bits(c[20].1),
        f64::from_bits(c[21].1),
        f64::from_bits(c[22].1),
        f64::from_bits(c[23].1),
        f64::from_bits(c[24].1),
        f64::from_bits(c[25].1),
        f64::from_bits(c[26].1),
        f64::from_bits(c[27].1),
    );

    let mut z = DoubleDouble::mul_f64_add(r, p, DoubleDouble::from_bit_pair(c[5]));
    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[4]));
    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[3]));
    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[2]));
    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[1]));
    z = DoubleDouble::mul_add(z, r, DoubleDouble::from_bit_pair(c[0]));
    let p = z;
    let err = f_fmla(
        p.hi,
        f64::from_bits(0x3c60000000000000), // 2^-57
        f64::from_bits(0x3c20000000000000), // 2^-61
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);
    if ub == lb {
        return p.to_f64();
    }
    y0_small_argument_moderate(r, c)
}

fn y0_small_argument_moderate(r: DoubleDouble, c0: &[(u64, u64); 28]) -> f64 {
    let c = &c0[15..];

    let p0 = f_polyeval13(
        r.to_f64(),
        f64::from_bits(c[0].1),
        f64::from_bits(c[1].1),
        f64::from_bits(c[2].1),
        f64::from_bits(c[3].1),
        f64::from_bits(c[4].1),
        f64::from_bits(c[5].1),
        f64::from_bits(c[6].1),
        f64::from_bits(c[7].1),
        f64::from_bits(c[8].1),
        f64::from_bits(c[9].1),
        f64::from_bits(c[10].1),
        f64::from_bits(c[11].1),
        f64::from_bits(c[12].1),
    );

    let c = c0;

    let mut p_e = DoubleDouble::mul_f64_add(r, p0, DoubleDouble::from_bit_pair(c[14]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[13]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[12]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[11]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[10]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[9]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[8]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[7]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[6]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[5]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[4]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[3]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[2]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[1]));
    p_e = DoubleDouble::mul_add(p_e, r, DoubleDouble::from_bit_pair(c[0]));

    let p = DoubleDouble::from_exact_add(p_e.hi, p_e.lo);
    let err = f_fmla(
        p.hi,
        f64::from_bits(0x3c30000000000000), // 2^-60
        f64::from_bits(0x3bf0000000000000), // 2^-64
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);
    if ub == lb {
        return p.to_f64();
    }
    y1_small_argument_hard(r, c)
}

#[cold]
#[inline(never)]
fn y1_small_argument_hard(r: DoubleDouble, c: &[(u64, u64); 28]) -> f64 {
    // if we're too close to zero taylor will converge faster and more accurate,
    // since remez, minimax and other almost cannot polynomial optimize near zero
    let mut p = DoubleDouble::from_bit_pair(c[27]);
    for i in (0..27).rev() {
        p = DoubleDouble::mul_add(r, p, DoubleDouble::from_bit_pair(c[i]));
        p = DoubleDouble::from_exact_add(p.hi, p.lo);
    }
    p.to_f64()
}

/*
   Evaluates:
   Y1 = sqrt(2/(PI*x)) * beta(x) * sin(x - 3*PI/4 - alpha(x))

   Discarding 1/2*PI gives:
   Y1 = sqrt(2/(PI*x)) * beta(x) * (-cos(x - PI/4 - alpha(x)))
*/
#[inline]
pub(crate) fn y1_asympt_fast(x: f64) -> f64 {
    const SQRT_2_OVER_PI: DoubleDouble = DoubleDouble::new(
        f64::from_bits(0xbc8cbc0d30ebfd15),
        f64::from_bits(0x3fe9884533d43651),
    );
    const MPI_OVER_4: DoubleDouble = DoubleDouble::new(
        f64::from_bits(0xbc81a62633145c07),
        f64::from_bits(0xbfe921fb54442d18),
    );

    let recip = if x.to_bits() > 0x7fd000000000000u64 {
        DoubleDouble::quick_mult_f64(DoubleDouble::from_exact_div_fma(4.0, x), 0.25)
    } else {
        DoubleDouble::from_recip(x)
    };

    let alpha = bessel_1_asympt_alpha_fast(recip);
    let beta = bessel_1_asympt_beta_fast(recip);

    let AngleReduced { angle } = rem2pi_any(x);

    // Without full subtraction cancellation happens sometimes
    let x0pi34 = DoubleDouble::full_dd_sub(MPI_OVER_4, alpha);
    let r0 = DoubleDouble::full_dd_add(angle, x0pi34);

    let m_cos = -cos_dd_small_fast(r0);
    let z0 = DoubleDouble::quick_mult(beta, m_cos);
    let r_sqrt = DoubleDouble::from_rsqrt_fast(x);
    let scale = DoubleDouble::quick_mult(SQRT_2_OVER_PI, r_sqrt);
    let r = DoubleDouble::quick_mult(scale, z0);
    let p = DoubleDouble::from_exact_add(r.hi, r.lo);
    let err = f_fmla(
        p.hi,
        f64::from_bits(0x3c40000000000000), // 2^-60
        f64::from_bits(0x3bf0000000000000), // 2^-87
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);
    if ub == lb {
        return p.to_f64();
    }
    y1_asympt(x, recip, r_sqrt, angle)
}

/*
   Evaluates:
   Y1 = sqrt(2/(PI*x)) * beta(x) * sin(x - 3*PI/4 - alpha(x))

   Discarding 1/2*PI gives:
   Y1 = sqrt(2/(PI*x)) * beta(x) * (-cos(x - PI/4 - alpha(x)))
*/
fn y1_asympt(x: f64, recip: DoubleDouble, r_sqrt: DoubleDouble, angle: DoubleDouble) -> f64 {
    const SQRT_2_OVER_PI: DoubleDouble = DoubleDouble::new(
        f64::from_bits(0xbc8cbc0d30ebfd15),
        f64::from_bits(0x3fe9884533d43651),
    );
    const MPI_OVER_4: DoubleDouble = DoubleDouble::new(
        f64::from_bits(0xbc81a62633145c07),
        f64::from_bits(0xbfe921fb54442d18),
    );

    let alpha = bessel_1_asympt_alpha(recip);
    let beta = bessel_1_asympt_beta(recip);

    // Without full subtraction cancellation happens sometimes
    let x0pi34 = DoubleDouble::full_dd_sub(MPI_OVER_4, alpha);
    let r0 = DoubleDouble::full_dd_add(angle, x0pi34);

    let m_cos = -cos_dd_small(r0);
    let z0 = DoubleDouble::quick_mult(beta, m_cos);
    let scale = DoubleDouble::quick_mult(SQRT_2_OVER_PI, r_sqrt);
    let r = DoubleDouble::quick_mult(scale, z0);
    let p = DoubleDouble::from_exact_add(r.hi, r.lo);
    let err = f_fmla(
        p.hi,
        f64::from_bits(0x3c30000000000000), // 2^-60
        f64::from_bits(0x3a80000000000000), // 2^-87
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);
    if ub == lb {
        return p.to_f64();
    }
    y1_asympt_hard(x)
}

/*
   Evaluates:
   Y1 = sqrt(2/(PI*x)) * beta(x) * sin(x - 3*PI/4 - alpha(x))

   Discarding 1/2*PI gives:
   Y1 = sqrt(2/(PI*x)) * beta(x) * (-cos(x - PI/4 - alpha(x)))
*/
#[cold]
#[inline(never)]
fn y1_asympt_hard(x: f64) -> f64 {
    const SQRT_2_OVER_PI: DyadicFloat128 = DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xcc42299e_a1b28468_7e59e280_5d5c7180_u128,
    };

    const MPI_OVER_4: DyadicFloat128 = DyadicFloat128 {
        sign: DyadicSign::Neg,
        exponent: -128,
        mantissa: 0xc90fdaa2_2168c234_c4c6628b_80dc1cd1_u128,
    };

    let x_dyadic = DyadicFloat128::new_from_f64(x);
    let recip = DyadicFloat128::accurate_reciprocal(x);

    let alpha = bessel_1_asympt_alpha_hard(recip);
    let beta = bessel_1_asympt_beta_hard(recip);

    let angle = rem2pi_f128(x_dyadic);

    let x0pi34 = MPI_OVER_4 - alpha;
    let r0 = angle + x0pi34;

    let m_cos = cos_f128_small(r0).negated();

    let z0 = beta * m_cos;
    let r_sqrt = bessel_rsqrt_hard(x, recip);
    let scale = SQRT_2_OVER_PI * r_sqrt;
    let p = scale * z0;
    p.fast_as_f64()
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_y1() {
        // 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
        assert_eq!(f_y1(0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000007291282546733975),
                   -873124540555277200000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.);
        assert_eq!(f_y1(2.1957931471395398), -0.0007023285780874727);
        assert_eq!(
            f_y1(f64::from_bits(0x571a31ffe2ff7e9fu64)),
            f64::from_bits(0x32e58532f95056ffu64)
        );
        assert_eq!(
            f_y1(f64::from_bits(0x400193bed4dff243)),
            0.00000000000000002513306678922122
        );
        assert_eq!(
            f_y1(f64::from_bits(0x3ffc513c569fe78e)),
            -0.24189760895998239
        );
        assert_eq!(
            f_y1(f64::from_bits(0x4192391e4c8faa60)),
            -0.000000000000000002572292246748134
        );
        assert_eq!(
            f_y1(f64::from_bits(0x403e9e480605283c)),
            -0.00000000000000001524456280251315
        );
        assert_eq!(
            f_y1(f64::from_bits(0x40277f9138d43206)),
            0.000000000000000006849807120770496
        );
        assert_eq!(f_y1(f64::INFINITY), 0.);
        assert!(f_y1(f64::NEG_INFINITY).is_nan());
        assert!(f_y1(f64::NAN).is_nan());
    }

    #[test]
    fn test_y1_edge_cases() {
        assert_eq!(f_y1(2.1904620854463985), -0.0034837351616785234);
        assert_eq!(f_y1(2.197142201034536), 4.5568985277260593e-7);
        assert_eq!(f_y1(1.4000000000000004), -0.4791469742327998);
        assert_eq!(f_y1(2.0002288794493848), -0.1069033735586767);
    }
}