open-standards/wai/wai-quantum/src/repro.rs

//! 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.

#[cfg(any(
    feature = "quantum_spd",
    feature = "quantum_pauli",
    feature = "quantum_mps",
    feature = "quantum_vml",
    feature = "quantum_phasor",
    feature = "quantum_kernel",
    feature = "quantum_tn",
    feature = "quantum_sv",
    feature = "quantum_tdvp",
    feature = "quantum_nqs"
))]
mod trig;
#[cfg(any(
    feature = "quantum_spd",
    feature = "quantum_pauli",
    feature = "quantum_mps",
    feature = "quantum_vml",
    feature = "quantum_phasor",
    feature = "quantum_kernel",
    feature = "quantum_tn",
    feature = "quantum_sv",
    feature = "quantum_tdvp",
    feature = "quantum_nqs"
))]
pub use trig::sin_cos;

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

#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg"))]
const LN2_HI: f64 = 6.931_471_803_691_238_164_90e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg"))]
const LN2_LO: f64 = 1.908_214_929_270_587_700_02e-10;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg"))]
const LG1: f64 = 6.666_666_666_666_735_130e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg"))]
const LG2: f64 = 3.999_999_999_940_941_908e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg"))]
const LG3: f64 = 2.857_142_874_366_239_149e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg"))]
const LG4: f64 = 2.222_219_843_214_978_396e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg"))]
const LG5: f64 = 1.818_357_216_161_805_012e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg"))]
const LG6: f64 = 1.531_383_769_920_937_332e-1;
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg"))]
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", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg"))]
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
}

// ---------------------------------------------------------------------------
// Exponential
// ---------------------------------------------------------------------------

#[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_nqs"))]
const EXP_P: [f64; 5] = [
    1.666_666_666_666_660_190_37e-1,
    -2.777_777_777_701_559_338_42e-3,
    6.613_756_321_437_934_361_17e-5,
    -1.653_390_220_546_525_153_90e-6,
    4.138_136_797_057_238_460_39e-8,
];

/// `eˣ`, bit-identical on every IEEE-754 platform: fdlibm's algorithm
/// (reduction by `ln 2` in two parts, a rational approximation on
/// `[−ln2/2, ln2/2]`), using only `+ − × ÷`. Within an ulp of the true value.
#[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_nqs"))]
pub fn exp(x: f64) -> f64 {
    const LN2_HI_E: f64 = 6.931_471_803_691_238_164_90e-1;
    const LN2_LO_E: f64 = 1.908_214_929_270_587_700_02e-10;
    const INV_LN2: f64 = 1.442_695_040_888_963_387_00;
    if x.is_nan() {
        return x;
    }
    if x > 7.097_827_128_933_839_730_96e2 {
        return f64::INFINITY;
    }
    if x < -7.451_332_191_019_411_084_20e2 {
        return 0.0;
    }
    let ax = x.abs();
    let (hi, lo, k): (f64, f64, i32);
    if ax > 0.5 * core::f64::consts::LN_2 {
        let kk = if ax < 1.5 * core::f64::consts::LN_2 {
            if x < 0.0 { -1 } else { 1 }
        } else {
            (INV_LN2 * x + if x < 0.0 { -0.5 } else { 0.5 }) as i32
        };
        let t = f64::from(kk);
        hi = x - t * LN2_HI_E;
        lo = t * LN2_LO_E;
        k = kk;
    } else if ax < 3.725_290_298_461_914e-9 {
        // |x| < 2⁻²⁸: 1 + x to the last bit.
        return 1.0 + x;
    } else {
        hi = x;
        lo = 0.0;
        k = 0;
    }
    let r = hi - lo;
    let t = r * r;
    let c = r - t * (EXP_P[0] + t * (EXP_P[1] + t * (EXP_P[2] + t * (EXP_P[3] + t * EXP_P[4]))));
    if k == 0 {
        return 1.0 - ((r * c) / (c - 2.0) - r);
    }
    let y = 1.0 - ((lo - (r * c) / (2.0 - c)) - hi);
    // y · 2ᵏ, in two steps where 2ᵏ alone would leave the normal range.
    let pow2 = |e: i32| f64::from_bits(((e + 1023) as u64) << 52);
    if k > 1023 {
        y * pow2(1023) * pow2(k - 1023)
    } else if k < -1021 {
        y * pow2(k + 1000) * pow2(-1000)
    } else {
        y * pow2(k)
    }
}

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

#[cfg(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))]
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(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))]
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(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))]
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(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))]
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(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))]
const PI: f64 = 3.141_592_653_589_793_116;
#[cfg(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))]
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(any(feature = "quantum_qfhrr", feature = "quantum_nqs"))]
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::*;

    #[cfg(any(feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_nqs"))]
    #[test]
    fn exp_is_within_an_ulp_of_the_platform() {
        let mut s = 0x9e37_79b9_7f4a_7c15u64;
        for i in 0..400_000 {
            s ^= s << 13;
            s ^= s >> 7;
            s ^= s << 17;
            let u = (s >> 11) as f64 / (1u64 << 53) as f64;
            let span = [1e-10, 0.3, 0.7, 2.0, 20.0, 700.0, 740.0][i % 7];
            let x = (2.0 * u - 1.0) * span;
            assert!(ulps(exp(x), x.exp()) <= 1, "exp({x:e}) = {:e} vs {:e}", exp(x), x.exp());
        }
        assert_eq!(exp(0.0), 1.0);
        assert_eq!(exp(f64::INFINITY), f64::INFINITY);
        assert_eq!(exp(f64::NEG_INFINITY), 0.0);
        assert!(exp(f64::NAN).is_nan());
        assert_eq!(exp(800.0), f64::INFINITY);
        assert_eq!(exp(-800.0), 0.0);
    }

    #[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", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg"))]
    #[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]);
    }
}