wai-quantum 0.3.36

//! `sin` and `cos`, for the layers that rotate: fdlibm's algorithm as musl
//! carries it, using only `+ − × ÷`.

/// `sin` and `cos` of `x`, bit-identical on every IEEE-754 platform. The
/// algorithm is fdlibm's (as musl carries it): Cody–Waite reduction by π/2 in
/// three parts, then minimax polynomials on [−π/4, π/4]. Only `+ − × ÷` are
/// used, and Rust never fuses them, so the bits do not depend on the target.
/// Accurate to within an ulp of the true value.
///
/// For `|x| > 2²⁰` (and for non-finite `x`) both are NaN: the reduction used
/// here is exact only below that, and a NaN is better than a value that
/// depends on the platform.
pub fn sin_cos(x: f64) -> (f64, f64) {
    if !(x.abs() <= 1_048_576.0) {
        return (f64::NAN, f64::NAN);
    }
    let ix = (x.to_bits() >> 32) as u32 & 0x7fff_ffff;
    if ix <= 0x3fe9_21fb {
        // |x| ≤ π/4
        if ix < 0x3e40_0000 {
            // |x| < 2⁻²⁷: sin x = x and cos x = 1 to the last bit.
            return (x, 1.0);
        }
        return (k_sin(x, 0.0, false), k_cos(x, 0.0));
    }
    let (n, y0, y1) = rem_pio2(x, ix);
    let (s, c) = (k_sin(y0, y1, true), k_cos(y0, y1));
    match n & 3 {
        0 => (s, c),
        1 => (c, -s),
        2 => (-s, -c),
        _ => (-c, s),
    }
}

const S1: f64 = -1.666_666_666_666_663_243_48e-1;
const S2: f64 = 8.333_333_333_322_489_461_24e-3;
const S3: f64 = -1.984_126_982_985_794_931_34e-4;
const S4: f64 = 2.755_731_370_707_006_767_89e-6;
const S5: f64 = -2.505_076_025_340_686_341_95e-8;
const S6: f64 = 1.589_690_995_211_550_102_21e-10;

fn k_sin(x: f64, y: f64, tail: bool) -> f64 {
    let z = x * x;
    let w = z * z;
    let r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6);
    let v = z * x;
    if !tail {
        x + v * (S1 + z * r)
    } else {
        x - ((z * (0.5 * y - v * r) - y) - v * S1)
    }
}

const C1: f64 = 4.166_666_666_666_660_190_37e-2;
const C2: f64 = -1.388_888_888_887_410_957_49e-3;
const C3: f64 = 2.480_158_728_947_672_941_78e-5;
const C4: f64 = -2.755_731_435_139_066_330_35e-7;
const C5: f64 = 2.087_572_321_298_174_827_90e-9;
const C6: f64 = -1.135_964_755_778_819_482_65e-11;

fn k_cos(x: f64, y: f64) -> f64 {
    let z = x * x;
    let w = z * z;
    let r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6));
    let hz = 0.5 * z;
    let w = 1.0 - hz;
    w + (((1.0 - w) - hz) + (z * r - x * y))
}

const TOINT: f64 = 1.5 / f64::EPSILON;
const PIO4: f64 = 7.853_981_633_974_482_789_99e-1;
const INVPIO2: f64 = 6.366_197_723_675_813_824_33e-1;
const PIO2_1: f64 = 1.570_796_326_734_125_614_17;
const PIO2_1T: f64 = 6.077_100_506_506_192_249_32e-11;
const PIO2_2: f64 = 6.077_100_506_303_965_976_60e-11;
const PIO2_2T: f64 = 2.022_266_248_795_950_631_54e-21;
const PIO2_3: f64 = 2.022_266_248_711_166_455_80e-21;
const PIO2_3T: f64 = 8.478_427_660_368_899_569_97e-32;

/// `x = n·π/2 + (y0 + y1)` with `|y0 + y1| ≤ π/4`, for `|x| ≤ 2²⁰`.
fn rem_pio2(x: f64, ix: u32) -> (i64, f64, f64) {
    let mut f = (x * INVPIO2 + TOINT) - TOINT;
    let mut r = x - f * PIO2_1;
    let mut w = f * PIO2_1T;
    // Matters only under directed rounding; kept so the port is exact.
    if r - w < -PIO4 {
        f -= 1.0;
        r = x - f * PIO2_1;
        w = f * PIO2_1T;
    } else if r - w > PIO4 {
        f += 1.0;
        r = x - f * PIO2_1;
        w = f * PIO2_1T;
    }
    let mut y0 = r - w;
    let ex = (ix >> 20) as i32;
    let ey = ((y0.to_bits() >> 52) & 0x7ff) as i32;
    if ex - ey > 16 {
        let t = r;
        w = f * PIO2_2;
        r = t - w;
        w = f * PIO2_2T - ((t - r) - w);
        y0 = r - w;
        let ey = ((y0.to_bits() >> 52) & 0x7ff) as i32;
        if ex - ey > 49 {
            let t = r;
            w = f * PIO2_3;
            r = t - w;
            w = f * PIO2_3T - ((t - r) - w);
            y0 = r - w;
        }
    }
    let y1 = (r - y0) - w;
    (f as i64, y0, y1)
}