pxfm 0.1.28

/*
 * // 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::alpha0::{
    bessel_0_asympt_alpha, bessel_0_asympt_alpha_fast, bessel_0_asympt_alpha_hard,
};
use crate::bessel::beta0::{
    bessel_0_asympt_beta, bessel_0_asympt_beta_fast, bessel_0_asympt_beta_hard,
};
use crate::bessel::i0::bessel_rsqrt_hard;
use crate::bessel::j0::j0_maclaurin_series;
use crate::bessel::y0_coeffs::Y0_COEFFS;
use crate::bessel::y0_coeffs_taylor::Y0_COEFFS_TAYLOR;
use crate::bessel::y0f_coeffs::{Y0_ZEROS, Y0_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::{sin_dd_small, sin_dd_small_fast, sin_f128_small};
use crate::sincos_reduce::{AngleReduced, rem2pi_any, rem2pi_f128};

/// Bessel of the second kind of order 0 (Y0)
pub fn f_y0(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 <= 0x4052d9999999999au64 {
        // x <= 75.4
        if xb <= 0x4000000000000000u64 {
            // x <= 2
            if xb <= 0x3ff599999999999au64 {
                // x < 1.35
                return y0_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 y0_transient_area_fast(x);
        }
        return y0_small_argument_fast(x);
    }

    y0_asympt_fast(x)
}

/**
Generated by SageMath:
Evaluates:
Y0(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))
expressed as:
Y0(x)=log(x)*W0(x) - Z0(x)
```python
from sage.all import *

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

# Define J0(x) Taylor expansion at x = 0
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)])

J0_series = j_series(0, x)

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

W0 = (d2/pi) * J0_series
Z0 = -gamma * (d2/pi) * J0_series + RealField(300)(2).log() * (d2/pi) * J0_series + (d2/pi) * z_series(x)

# see the series
print(W0)
print(Z0)
```
**/
#[inline]
fn y0_near_zero_fast(x: f64) -> f64 {
    const W: [(u64, u64); 15] = [
        (0xbc86b01ec5417056, 0x3fe45f306dc9c883),
        (0x3c66b01ec5417056, 0xbfc45f306dc9c883),
        (0xbc26b01ec5417056, 0x3f845f306dc9c883),
        (0xbbd67fe4a5feb897, 0xbf321bb945252402),
        (0x3b767fe4a5feb897, 0x3ed21bb945252402),
        (0xbaf5c2495706f745, 0xbe672db9f21b0f5f),
        (0x3a90c8209874dfad, 0x3df49a6c656d62ff),
        (0x3a12921e91b07dd0, 0xbd7ae90af76a4d0f),
        (0xb992921e91b07dd0, 0x3cfae90af76a4d0f),
        (0x39089b0d8a9228ca, 0xbc754331c053fdad),
        (0x3878d321ddfd3c6e, 0x3beb3749ebf0a0dd),
        (0x37e77548130d809b, 0xbb5cca5ae46eae67),
        (0xb73a848e7ca1c943, 0x3ac9976d3cd4293f),
        (0xb6c884706195a054, 0xba336206ff1ce731),
        (0xb6387a7d2389630d, 0x39995103e9f1818f),
    ];
    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]));

    const Z: [(u64, u64); 15] = [
        (0xbc5ddfd831a70821, 0x3fb2e4d699cbd01f),
        (0xbc6d93e63489aea6, 0xbfc6bbcb41034286),
        (0xbc1b88525c2e130b, 0x3f9075b1bbf41364),
        (0x3be097334e26e578, 0xbf41a6206b7b973d),
        (0x3b51c64a34c78cda, 0x3ee3e99794203bbd),
        (0xbb1c407b0f5b2805, 0xbe7bce4a600d3ea4),
        (0xbaa57d1e1e88c9ca, 0x3e0a6ee796b871b6),
        (0x3a3b6e7030a77899, 0xbd92393d82c6b2e4),
        (0x397fcfedacb03781, 0x3d131085da82054c),
        (0xb8e45f51f6118e46, 0xbc8f4ed4b492ebcc),
        (0xb89bd46046c3c8de, 0x3c04b7ac8a1b15d0),
        (0x37d1a206fb205e32, 0xbb769201941d0d49),
        (0x3782f38acbf23993, 0x3ae4987e587ab039),
        (0x36b691bdabf5672b, 0xba4ff1953e0a7c5b),
        (0x3636e1c8cd260e18, 0x39b55031dc5e1967),
    ];
    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]));
    let w_log = log_dd_fast(x);
    let p = DoubleDouble::mul_add(w, w_log, -z);
    let err = f_fmla(
        p.hi,
        f64::from_bits(0x3c50000000000000), // 2^-58
        f64::from_bits(0x3c30000000000000), // 2^-60
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);
    if ub == lb {
        return p.to_f64();
    }
    y0_near_zero(x, w_log)
}

/**
Generated by SageMath:
Evaluates:
Y0(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))
expressed as:
Y0(x)=log(x)*W0(x) - Z0(x)
```python
from sage.all import *

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

# Define J0(x) Taylor expansion at x = 0
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)])

J0_series = j_series(0, x)

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

W0 = (d2/pi) * J0_series
Z0 = -gamma * (d2/pi) * J0_series + RealField(300)(2).log() * (d2/pi) * J0_series + (d2/pi) * z_series(x)

# see the series
print(W0)
print(Z0)
```
**/
#[cold]
#[inline(never)]
fn y0_near_zero(x: f64, w_log: DoubleDouble) -> f64 {
    const W: [(u64, u64); 15] = [
        (0xbc86b01ec5417056, 0x3fe45f306dc9c883),
        (0x3c66b01ec5417056, 0xbfc45f306dc9c883),
        (0xbc26b01ec5417056, 0x3f845f306dc9c883),
        (0xbbd67fe4a5feb897, 0xbf321bb945252402),
        (0x3b767fe4a5feb897, 0x3ed21bb945252402),
        (0xbaf5c2495706f745, 0xbe672db9f21b0f5f),
        (0x3a90c8209874dfad, 0x3df49a6c656d62ff),
        (0x3a12921e91b07dd0, 0xbd7ae90af76a4d0f),
        (0xb992921e91b07dd0, 0x3cfae90af76a4d0f),
        (0x39089b0d8a9228ca, 0xbc754331c053fdad),
        (0x3878d321ddfd3c6e, 0x3beb3749ebf0a0dd),
        (0x37e77548130d809b, 0xbb5cca5ae46eae67),
        (0xb73a848e7ca1c943, 0x3ac9976d3cd4293f),
        (0xb6c884706195a054, 0xba336206ff1ce731),
        (0xb6387a7d2389630d, 0x39995103e9f1818f),
    ];
    let x2 = DoubleDouble::from_exact_mult(x, x);
    let 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]),
    );

    const Z: [(u64, u64); 15] = [
        (0xbc5ddfd831a70821, 0x3fb2e4d699cbd01f),
        (0xbc6d93e63489aea6, 0xbfc6bbcb41034286),
        (0xbc1b88525c2e130b, 0x3f9075b1bbf41364),
        (0x3be097334e26e578, 0xbf41a6206b7b973d),
        (0x3b51c64a34c78cda, 0x3ee3e99794203bbd),
        (0xbb1c407b0f5b2805, 0xbe7bce4a600d3ea4),
        (0xbaa57d1e1e88c9ca, 0x3e0a6ee796b871b6),
        (0x3a3b6e7030a77899, 0xbd92393d82c6b2e4),
        (0x397fcfedacb03781, 0x3d131085da82054c),
        (0xb8e45f51f6118e46, 0xbc8f4ed4b492ebcc),
        (0xb89bd46046c3c8de, 0x3c04b7ac8a1b15d0),
        (0x37d1a206fb205e32, 0xbb769201941d0d49),
        (0x3782f38acbf23993, 0x3ae4987e587ab039),
        (0x36b691bdabf5672b, 0xba4ff1953e0a7c5b),
        (0x3636e1c8cd260e18, 0x39b55031dc5e1967),
    ];
    let 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]),
    );
    DoubleDouble::mul_add(w, w_log, -z).to_f64()
}

/**
Path for transient area between 1.35 to 2.
**/
#[inline]
pub(crate) fn y0_transient_area_fast(x: f64) -> f64 {
    /**
    Polynomial generated by Wolfram:
    ```text
    <<FunctionApproximations`
    ClearAll["Global`*"]
    f[x_]:= BesselY[0,x + 2.1971413260310170351490335626990]
    {approx,error} = MiniMaxApproximation[f[x],{x,{1.35 - 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] = [
        (0xbc689e111675434b, 0x3fe0aa48442f014b),
        (0x396ffb11562e8c70, 0x3cc1806f07aceb3c),
        (0xbc6156edff56513d, 0xbfd0aa48442f0030),
        (0xbc278dbff5ee7db4, 0x3fa439fac165db2d),
        (0x3c1f9463c2023663, 0x3f80d2af4ebc8fc4),
        (0xbbee09e5733a1236, 0x3f4f716488aebd9c),
        (0xbbf261a1f255ddf4, 0xbf5444bd2a7e6254),
        (0x3bd4f543544f1fe7, 0x3f384c300e8674d8),
        (0xbb96ef8d95fe049b, 0xbf217a0fc83af41e),
        (0x3ba2be82573ae98d, 0x3f0dbc664048e495),
        (0xbb942b15646c85f2, 0xbef8522f83e4a3e3),
        (0x3b7a127725ba4606, 0x3ee775c010ce4146),
        (0x3ae4a02f0b2a18e2, 0x3e7d1d7cf40f9697),
        (0x3b8fcf1a3d27236b, 0x3eea9c226c21712d),
        (0xbb70b3aa0a1e9ffb, 0x3ef8f237831ec74b),
        (0x3baa6c24261245f3, 0x3f08ebfea3ea469e),
        (0x3bb5fa1b8c4587c4, 0x3f1474022b90cbda),
        (0x3b69545db8d098d1, 0x3f1d153dc04c51c0),
        (0x3bc68eab6520d21b, 0x3f2198a4578cb6ca),
        (0xbbc255734bc49c8b, 0x3f2212febf7cecdd),
        (0x3bb8dd02722339f5, 0x3f1f314deec17049),
        (0x3bbb6ef8f04b26a2, 0x3f1657d051699088),
        (0x3b878233fbf501dc, 0x3f0a1a422dafcef6),
        (0xbb73730138d1dbc2, 0x3ef8423cd021f1dd),
        (0x3b7e2a0d7009d709, 0x3ee145cae37afe1b),
        (0x3b6af21aeaba4e57, 0x3ec1bc74380f6d7b),
        (0xbb3607fb9242657f, 0x3e977341fc10cdc8),
        (0xbac747923880f651, 0x3e5e30218bc1fee3),
    ];

    const ZERO: DoubleDouble =
        DoubleDouble::from_bit_pair((0xbc8bd1e50d219bfd, 0x400193bed4dff243));

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

    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(0x3b10000000000000), // 2^-78
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);
    if ub != lb {
        return y0_transient_area_moderate(x);
    }
    p.to_f64()
}

/**
Path for transient area between 1.35 to 2.
**/
fn y0_transient_area_moderate(x: f64) -> f64 {
    /**
    Polynomial generated by Wolfram:
    ```text
    <<FunctionApproximations`
    ClearAll["Global`*"]
    f[x_]:= BesselY[0,x + 2.1971413260310170351490335626990]
    {approx,error} = MiniMaxApproximation[f[x],{x,{1.35 - 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] = [
        (0xbc689e111675434b, 0x3fe0aa48442f014b),
        (0x396ffb11562e8c70, 0x3cc1806f07aceb3c),
        (0xbc6156edff56513d, 0xbfd0aa48442f0030),
        (0xbc278dbff5ee7db4, 0x3fa439fac165db2d),
        (0x3c1f9463c2023663, 0x3f80d2af4ebc8fc4),
        (0xbbee09e5733a1236, 0x3f4f716488aebd9c),
        (0xbbf261a1f255ddf4, 0xbf5444bd2a7e6254),
        (0x3bd4f543544f1fe7, 0x3f384c300e8674d8),
        (0xbb96ef8d95fe049b, 0xbf217a0fc83af41e),
        (0x3ba2be82573ae98d, 0x3f0dbc664048e495),
        (0xbb942b15646c85f2, 0xbef8522f83e4a3e3),
        (0x3b7a127725ba4606, 0x3ee775c010ce4146),
        (0x3ae4a02f0b2a18e2, 0x3e7d1d7cf40f9697),
        (0x3b8fcf1a3d27236b, 0x3eea9c226c21712d),
        (0xbb70b3aa0a1e9ffb, 0x3ef8f237831ec74b),
        (0x3baa6c24261245f3, 0x3f08ebfea3ea469e),
        (0x3bb5fa1b8c4587c4, 0x3f1474022b90cbda),
        (0x3b69545db8d098d1, 0x3f1d153dc04c51c0),
        (0x3bc68eab6520d21b, 0x3f2198a4578cb6ca),
        (0xbbc255734bc49c8b, 0x3f2212febf7cecdd),
        (0x3bb8dd02722339f5, 0x3f1f314deec17049),
        (0x3bbb6ef8f04b26a2, 0x3f1657d051699088),
        (0x3b878233fbf501dc, 0x3f0a1a422dafcef6),
        (0xbb73730138d1dbc2, 0x3ef8423cd021f1dd),
        (0x3b7e2a0d7009d709, 0x3ee145cae37afe1b),
        (0x3b6af21aeaba4e57, 0x3ec1bc74380f6d7b),
        (0xbb3607fb9242657f, 0x3e977341fc10cdc8),
        (0xbac747923880f651, 0x3e5e30218bc1fee3),
    ];

    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_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(0x3c10000000000000), // 2^-62
        f64::from_bits(0x3a30000000000000), // 2^-91
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);
    if ub != lb {
        return y0_transient_area_hard(x);
    }
    p.to_f64()
}

#[cold]
#[inline(never)]
fn y0_transient_area_hard(x: f64) -> f64 {
    const ZERO: DyadicFloat128 = DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -126,
        mantissa: 0x8c9df6a6_ff921721_70d796f3_2017e155_u128,
    };
    let r = DyadicFloat128::new_from_f64(x) - ZERO;
    /*
    Polynomial generated by Wolfram:
    ```text
    <<FunctionApproximations`
    ClearAll["Global`*"]
    f[x_]:= BesselY[0,x + 2.1971413260310170351490335626990]
    {approx,error} = MiniMaxApproximation[f[x],{x,{1.35 - 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}]]
    ```
    */
    static C: [DyadicFloat128; 28] = [
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -128,
            mantissa: 0x85524221_780a573b_0f774c55_e5a946dc_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -178,
            mantissa: 0x8c03783d_6759e3ff_622ac5d1_8df27811_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -129,
            mantissa: 0x85524221_78018115_6edff565_13cf55ab_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -132,
            mantissa: 0xa1cfd60b_2ed96743_9200508c_125cf3d8_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -134,
            mantissa: 0x86957a75_e47e21f9_463c2023_663178df_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -138,
            mantissa: 0xfb8b2445_75ecdc3e_c35198bd_b9330775_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -137,
            mantissa: 0xa225e953_f312a24c_343e4abb_be73338f_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -139,
            mantissa: 0xc2618074_33a6c29e_a86a89e3_fcde0075_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -140,
            mantissa: 0x8bd07e41_d7a0f05b_be3657f8_126b9715_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -142,
            mantissa: 0xede33202_4724aa57_d04ae75d_31ad69cd_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -143,
            mantissa: 0xc2917c1f_251f1a85_62ac8d90_be4550ff_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -144,
            mantissa: 0xbbae0086_720a31a1_27725ba4_605e0479_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -151,
            mantissa: 0xe8ebe7a0_7cb4b852_80bc2ca8_63864ee0_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -144,
            mantissa: 0xd4e11361_0b896bf9_e347a4e4_6d642f8e_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -143,
            mantissa: 0xc791bc18_f63a577a_62afaf0b_002753e2_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -142,
            mantissa: 0xc75ff51f_5234f34d_8484c248_be6beef2_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -141,
            mantissa: 0xa3a0115c_865ed2bf_437188b0_f87883dc_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -141,
            mantissa: 0xe8a9ee02_628e0019_545db8d0_98d12eff_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -140,
            mantissa: 0x8cc522bc_65b652d1_d56ca41a_43537f6a_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -140,
            mantissa: 0x9097f5fb_e766e5b5_5196876c_6ea798ce_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -141,
            mantissa: 0xf98a6f76_0b824b1b_a04e4467_3ea487e1_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -141,
            mantissa: 0xb2be828b_4c84436d_df1e0964_d44f420f_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -142,
            mantissa: 0xd0d2116d_7e77b0bc_119fdfa8_0ee2dfee_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -143,
            mantissa: 0xc211e681_0f8ee764_67f63971_21f1a199_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -144,
            mantissa: 0x8a2e571b_d7f0d9e2_a0d7009d_70969fa2_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -146,
            mantissa: 0x8de3a1c0_7b6bdb5e_435d5749_cadf3edd_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -149,
            mantissa: 0xbb9a0fe0_866e3d3f_008db7b3_5029fe59_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -153,
            mantissa: 0xf1810c5e_0ff717a2_e1b71dfc_26babb9f_u128,
        },
    ];
    let mut z = C[27];
    for i in (0..27).rev() {
        z = r * z + C[i];
    }
    z.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 y0_small_argument_fast(x: f64) -> f64 {
    let x_abs = x;

    // let avg_step = 74.607799 / 47.0;
    // let inv_step = 1.0 / avg_step;

    const INV_STEP: f64 = 0.6299609508652038;

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

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

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

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

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

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

    let c = if is_too_close_to_zero {
        &Y0_COEFFS_TAYLOR[idx - 1]
    } else {
        &Y0_COEFFS[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(Y0_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 y0_small_argument_moderate(r, c);
    }
    z.to_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.
fn y0_small_argument_moderate(r: DoubleDouble, c: &[(u64, u64); 28]) -> f64 {
    let c0 = &c[15..];

    let p0 = f_polyeval13(
        r.to_f64(),
        f64::from_bits(c0[0].1),
        f64::from_bits(c0[1].1),
        f64::from_bits(c0[2].1),
        f64::from_bits(c0[3].1),
        f64::from_bits(c0[4].1),
        f64::from_bits(c0[5].1),
        f64::from_bits(c0[6].1),
        f64::from_bits(c0[7].1),
        f64::from_bits(c0[8].1),
        f64::from_bits(c0[9].1),
        f64::from_bits(c0[10].1),
        f64::from_bits(c0[11].1),
        f64::from_bits(c0[12].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(0x3c30000000000000), // 2^-60
        f64::from_bits(0x3a30000000000000), // 2^-91
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);
    if ub != lb {
        return y0_small_argument_hard(r, c);
    }
    p.to_f64()
}

#[cold]
#[inline(never)]
fn y0_small_argument_hard(r: DoubleDouble, c: &[(u64, u64); 28]) -> f64 {
    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:
   Y0 = sqrt(2/(PI*x)) * beta(x) * sin(x - PI/4 - alpha(x))
*/
#[inline]
pub(crate) fn y0_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_0_asympt_alpha_fast(recip);
    let beta = bessel_0_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 = sin_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^-59
        f64::from_bits(0x3c10000000000000), // 2^-62
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);

    if ub == lb {
        return p.to_f64();
    }
    y0_asympt(x, recip, r_sqrt, angle)
}

/*
   Evaluates:
   Y0 = sqrt(2/(PI*x)) * beta(x) * sin(x - PI/4 - alpha(x))
*/
fn y0_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_0_asympt_alpha(recip);
    let beta = bessel_0_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 = sin_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(0x3a20000000000000), // 2^-93
    );
    let ub = p.hi + (p.lo + err);
    let lb = p.hi + (p.lo - err);

    if ub == lb {
        return p.to_f64();
    }
    y0_asympt_hard(x)
}

/*
   Evaluates:
   Y0 = sqrt(2/(PI*x)) * beta(x) * sin(x - PI/4 - alpha(x))
*/
#[cold]
#[inline(never)]
fn y0_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_0_asympt_alpha_hard(recip);
    let beta = bessel_0_asympt_beta_hard(recip);

    let angle = rem2pi_f128(x_dyadic);

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

    let m_sin = sin_f128_small(r0);

    let z0 = beta * m_sin;
    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_y0() {
        //ULP should be less than 0.5, but it was 1017.1969361449036, on 3 result 0.37685001001284685, using f_y0 and MPFR 0.3768500100127904
        assert_eq!(f_y0(3.), 0.3768500100127904);
        assert_eq!(f_y0(0.906009703874588), 0.01085796448629276);
        assert_eq!(f_y0(80.), -0.05562033908977);
        assert_eq!(f_y0(5.), -0.30851762524903376);
        assert_eq!(
            f_y0(f64::from_bits(0x3fec982eb8d417ea)),
            -0.000000000000000023389279284062102
        );
        assert!(f_y0(f64::NAN).is_nan());
        assert_eq!(f_y0(f64::INFINITY), 0.);
        assert!(f_y0(f64::NEG_INFINITY).is_nan());
    }

    #[test]
    fn test_y0_edge_values() {
        assert_eq!(f_y0(0.8900000000138676), -0.0031519646708080126);
        assert_eq!(f_y0(0.8900000000409116), -0.0031519646469294936);
        assert_eq!(f_y0(98.1760435789366), 0.0000000000000056889416242533015);
        assert_eq!(
            f_y0(91.8929453121571802176),
            -0.00000000000000007281665706677893
        );
        assert_eq!(
            f_y0(f64::from_bits(0x6e7c1d741dc52512u64)),
            f64::from_bits(0x2696f860815bc669)
        );
        assert_eq!(f_y0(f64::from_bits(0x3e04cdee58a47edd)), -13.58605001628649);
        assert_eq!(
            f_y0(0.89357696627916749),
            -0.000000000000000023389279284062102
        );
    }
}