cinque.rs

const PIO2_HI: f64 = core::f64::consts::FRAC_PI_2;

// pi/2 - nearest-f64(pi/2). Important for cos_q near pi/2.
const PIO2_LO: f64 = 6.123_233_995_736_766e-17;

// SLEEF full-quadrant polynomial coefficients.
const A0: f64 =  8.333_333_333_333_329_748_238_15e-3;
const A1: f64 = -1.984_126_984_126_961_628_068_09e-4;
const A2: f64 =  2.755_731_922_391_987_476_304_16e-6;
const A3: f64 = -2.505_210_837_635_020_458_107_55e-8;
const A4: f64 =  1.605_904_306_056_645_016_290_54e-10;
const A5: f64 = -7.647_122_191_181_588_332_884_84e-13;
const A6: f64 =  2.810_099_727_108_632_000_912_51e-15;
const A7: f64 = -7.972_559_550_090_378_688_919_52e-18;

const S1: f64 = -1.666_666_666_666_666_574_148_08e-1;

#[inline(always)]
fn sin_quadrant_kernel(x: f64) -> f64 {
    let z  = x * x;
    let z2 = z * z;
    let z4 = z2 * z2;

    // Four independent FMAs.
    let p0 = A1.mul_add(z, A0);
    let p1 = A3.mul_add(z, A2);
    let p2 = A5.mul_add(z, A4);
    let p3 = A7.mul_add(z, A6);

    // Estrin reduction.
    let q0 = p1.mul_add(z2, p0);
    let q1 = p3.mul_add(z2, p2);
    let p  = q1.mul_add(z4, q0);
    let p  = p.mul_add(z, S1);

    (x * z).mul_add(p, x)
}

#[inline(always)]
pub fn sin_q(x: f64) -> f64 {
    debug_assert!(x >= 0.0 && x <= PIO2_HI);
    sin_quadrant_kernel(x)
}

#[inline(always)]
pub fn cos_q(x: f64) -> f64 {
    debug_assert!(x >= 0.0 && x <= PIO2_HI);
    sin_quadrant_kernel((PIO2_HI - x) + PIO2_LO)
}