wai-quantum 0.3.23

//! Transcendental functions that give the same bits on every platform.
//!
//! The platform's libm may differ from another's in the last place. That is
//! enough to break a result promised bit-identical natively, in the browser and
//! at the edge, and it did: four layers of this crate computed different bits
//! natively and as a WASI component until they moved here. Everything below is
//! fdlibm's algorithm, as musl carries it, using only `+ − × ÷`, which Rust
//! never fuses.

/// `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)
}

// ---------------------------------------------------------------------------
// Natural logarithm
// ---------------------------------------------------------------------------

#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))]
const LN2_HI: f64 = 6.931_471_803_691_238_164_90e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))]
const LN2_LO: f64 = 1.908_214_929_270_587_700_02e-10;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))]
const LG1: f64 = 6.666_666_666_666_735_130e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))]
const LG2: f64 = 3.999_999_999_940_941_908e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))]
const LG3: f64 = 2.857_142_874_366_239_149e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))]
const LG4: f64 = 2.222_219_843_214_978_396e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))]
const LG5: f64 = 1.818_357_216_161_805_012e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))]
const LG6: f64 = 1.531_383_769_920_937_332e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))]
const LG7: f64 = 1.479_819_860_511_658_591e-1;

/// The natural logarithm, within an ulp of the true value: `−∞` at zero, NaN
/// below it.
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))]
pub fn ln(x: f64) -> f64 {
    let mut bits = x.to_bits();
    let mut hx = (bits >> 32) as u32;
    let mut k: i32 = 0;
    let mut x = x;
    if hx < 0x0010_0000 || hx >> 31 != 0 {
        if bits << 1 == 0 {
            return f64::NEG_INFINITY;
        }
        if hx >> 31 != 0 {
            return f64::NAN;
        }
        // subnormal: scale up
        k -= 54;
        x *= 18_014_398_509_481_984.0; // 2⁵⁴
        bits = x.to_bits();
        hx = (bits >> 32) as u32;
    } else if hx >= 0x7ff0_0000 {
        return x;
    } else if hx == 0x3ff0_0000 && bits << 32 == 0 {
        return 0.0;
    }
    // reduce x into [√2/2, √2]
    hx = hx.wrapping_add(0x3ff0_0000 - 0x3fe6_a09e);
    k += (hx >> 20) as i32 - 0x3ff;
    hx = (hx & 0x000f_ffff) + 0x3fe6_a09e;
    let x = f64::from_bits((u64::from(hx) << 32) | (bits & 0xffff_ffff));
    let f = x - 1.0;
    let hfsq = 0.5 * f * f;
    let s = f / (2.0 + f);
    let z = s * s;
    let w = z * z;
    let t1 = w * (LG2 + w * (LG4 + w * LG6));
    let t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7)));
    let r = t2 + t1;
    let dk = f64::from(k);
    s * (hfsq + r) + dk * LN2_LO - hfsq + f + dk * LN2_HI
}

// ---------------------------------------------------------------------------
// Arctangent
// ---------------------------------------------------------------------------

#[cfg(feature = "quantum_qfhrr")]
const ATANHI: [f64; 4] = [
    4.636_476_090_008_060_935_15e-1,
    7.853_981_633_974_482_789_99e-1,
    9.827_937_232_473_290_540_82e-1,
    1.570_796_326_794_896_558_00,
];
#[cfg(feature = "quantum_qfhrr")]
const ATANLO: [f64; 4] = [
    2.269_877_745_296_168_709_24e-17,
    3.061_616_997_868_383_017_93e-17,
    1.390_331_103_123_099_845_16e-17,
    6.123_233_995_736_766_035_87e-17,
];
#[cfg(feature = "quantum_qfhrr")]
const AT: [f64; 11] = [
    3.333_333_333_333_293_180_27e-1,
    -1.999_999_999_987_648_324_76e-1,
    1.428_571_427_250_346_637_11e-1,
    -1.111_111_040_546_235_578_80e-1,
    9.090_887_133_436_506_561_96e-2,
    -7.691_876_205_044_829_994_95e-2,
    6.661_073_137_387_531_206_69e-2,
    -5.833_570_133_790_573_486_45e-2,
    4.976_877_994_615_932_360_17e-2,
    -3.653_157_274_421_691_552_70e-2,
    1.628_582_011_536_578_236_23e-2,
];

/// The arctangent, within an ulp of the true value.
#[cfg(feature = "quantum_qfhrr")]
pub fn atan(x: f64) -> f64 {
    let ix = (x.to_bits() >> 32) as u32;
    let negative = ix >> 31 != 0;
    let ix = ix & 0x7fff_ffff;
    if ix >= 0x4410_0000 {
        // |x| ≥ 2⁶⁶
        if x.is_nan() {
            return x;
        }
        let z = ATANHI[3];
        return if negative { -z } else { z };
    }
    let mut x = x;
    let id: i32;
    if ix < 0x3fdc_0000 {
        // |x| < 0.4375
        if ix < 0x3e40_0000 {
            return x;
        }
        id = -1;
    } else {
        x = x.abs();
        if ix < 0x3ff3_0000 {
            if ix < 0x3fe6_0000 {
                id = 0;
                x = (2.0 * x - 1.0) / (2.0 + x);
            } else {
                id = 1;
                x = (x - 1.0) / (x + 1.0);
            }
        } else if ix < 0x4003_8000 {
            id = 2;
            x = (x - 1.5) / (1.0 + 1.5 * x);
        } else {
            id = 3;
            x = -1.0 / x;
        }
    }
    let z = x * x;
    let w = z * z;
    let s1 = z * (AT[0] + w * (AT[2] + w * (AT[4] + w * (AT[6] + w * (AT[8] + w * AT[10])))));
    let s2 = w * (AT[1] + w * (AT[3] + w * (AT[5] + w * (AT[7] + w * AT[9]))));
    if id < 0 {
        return x - x * (s1 + s2);
    }
    let i = id as usize;
    let z = ATANHI[i] - (x * (s1 + s2) - ATANLO[i] - x);
    if negative { -z } else { z }
}

#[cfg(feature = "quantum_qfhrr")]
const PI: f64 = 3.141_592_653_589_793_116;
#[cfg(feature = "quantum_qfhrr")]
const PI_LO: f64 = 1.224_646_799_147_353_177_2e-16;

/// The angle of the point `(x, y)`, in `(−π, π]`, within an ulp of the true
/// value.
#[cfg(feature = "quantum_qfhrr")]
pub fn atan2(y: f64, x: f64) -> f64 {
    if x.is_nan() || y.is_nan() {
        return x + y;
    }
    let (xb, yb) = (x.to_bits(), y.to_bits());
    let (ix, lx) = ((xb >> 32) as u32, xb as u32);
    let (iy, ly) = ((yb >> 32) as u32, yb as u32);
    if (ix.wrapping_sub(0x3ff0_0000) | lx) == 0 {
        return atan(y); // x = 1
    }
    let m = ((iy >> 31) & 1) | ((ix >> 30) & 2);
    let (ix, iy) = (ix & 0x7fff_ffff, iy & 0x7fff_ffff);
    if (iy | ly) == 0 {
        return match m {
            0 | 1 => y,
            2 => PI,
            _ => -PI,
        };
    }
    if (ix | lx) == 0 {
        return if m & 1 != 0 { -PI / 2.0 } else { PI / 2.0 };
    }
    if ix == 0x7ff0_0000 {
        return if iy == 0x7ff0_0000 {
            match m {
                0 => PI / 4.0,
                1 => -PI / 4.0,
                2 => 3.0 * PI / 4.0,
                _ => -3.0 * PI / 4.0,
            }
        } else {
            match m {
                0 => 0.0,
                1 => -0.0,
                2 => PI,
                _ => -PI,
            }
        };
    }
    if ix.wrapping_add(64 << 20) < iy || iy == 0x7ff0_0000 {
        return if m & 1 != 0 { -PI / 2.0 } else { PI / 2.0 };
    }
    let z = if m & 2 != 0 && iy.wrapping_add(64 << 20) < ix { 0.0 } else { atan((y / x).abs()) };
    match m {
        0 => z,
        1 => -z,
        2 => PI - (z - PI_LO),
        _ => (z - PI_LO) - PI,
    }
}

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

    #[allow(dead_code)]
    fn ulps(a: f64, b: f64) -> u64 {
        if a == b {
            return 0;
        }
        (a.to_bits() as i64).abs_diff(b.to_bits() as i64)
    }

    /// A deterministic spread of inputs over many magnitudes.
    #[allow(dead_code)]
    fn inputs() -> impl Iterator<Item = f64> {
        let mut s = 0x2545_f491_4f6c_dd1du64;
        (0..200_000).map(move |i| {
            s ^= s << 13;
            s ^= s >> 7;
            s ^= s << 17;
            let u = (s >> 11) as f64 / (1u64 << 53) as f64;
            let scale = [1e-300, 1e-8, 1e-3, 0.5, 1.0, 3.0, 1e3, 1e12, 1e300][i % 9];
            u * scale
        })
    }

    #[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))]
    #[test]
    fn ln_is_within_an_ulp_of_the_platform() {
        for x in inputs().filter(|&x| x > 0.0) {
            assert!(ulps(ln(x), x.ln()) <= 1, "ln({x:e}) = {:e} vs {:e}", ln(x), x.ln());
        }
        assert_eq!(ln(1.0), 0.0);
        assert_eq!(ln(0.0), f64::NEG_INFINITY);
        assert!(ln(-1.0).is_nan());
        assert_eq!(ln(f64::INFINITY), f64::INFINITY);
        assert!(ulps(ln(f64::MIN_POSITIVE / 8.0), (f64::MIN_POSITIVE / 8.0).ln()) <= 1, "subnormal");
    }

    #[cfg(feature = "quantum_qfhrr")]
    #[test]
    fn atan2_is_within_an_ulp_of_the_platform() {
        let xs: Vec<f64> = inputs().take(4000).collect();
        for (i, &a) in xs.iter().enumerate() {
            let b = xs[(i * 7 + 3) % xs.len()];
            for (y, x) in [(a, b), (-a, b), (a, -b), (-a, -b)] {
                assert!(ulps(atan2(y, x), y.atan2(x)) <= 1, "atan2({y:e}, {x:e}) = {:e} vs {:e}", atan2(y, x), y.atan2(x));
            }
        }
        assert_eq!(atan2(0.0, -1.0), PI);
        assert_eq!(atan2(1.0, 0.0), PI / 2.0);
        assert_eq!(atan(f64::INFINITY), ATANHI[3]);
    }
}