open-standards/wai/wai-rs/src/float_repro.rs

//! Float functions whose results are fixed by this source, not by the platform.
//!
//! IEEE-754 requires addition, subtraction, multiplication, division, square root and
//! the conversions to be correctly rounded, so they give the same bits for the same
//! operands on every conforming target. A platform math library's transcendental
//! functions (`exp`, `ln`, `pow`, `tanh`, ...) are not held to that, and their last
//! bits differ between libraries. Rounding to an integer (`round`,
//! `round_ties_even`) is exactly specified, but whether it compiles to an instruction
//! or to a call depends on the target. With rustc 1.98, baseline x86_64 compiles
//! `f64::round`, `f32::round` and `f64::round_ties_even` to calls to the symbols
//! `round`, `roundf` and `rint`, and wasm32 compiles the first two to calls to `round`
//! and `roundf`; AArch64 uses instructions. What answers such a call is settled when
//! the executable is linked, not by the platform's math library as such: the x86_64
//! Linux test executable inspected defined `round`, `roundf` and `rint` itself, as weak
//! symbols from Rust's compiler builtins, and imported only `exp` and `expf` from the C
//! library. The compiler builtins for wasm32-wasip2 define the same three; which
//! definition a wasm32 link takes was not checked.
//!
//! The functions here use only those basic operations, integer conversions and bit
//! manipulation, evaluated in the order written. Rust neither fuses a multiply and an
//! add into one rounding nor reassociates unless the source asks for it, so each
//! function returns the same bits on every target that evaluates at the operand's own
//! precision (x86_64, AArch64, wasm32; not x87) in the IEEE-754 default
//! floating-point environment (round to nearest even, subnormals kept). The JPEG AI
//! exact decoder ([`crate::jpegai_decode`]) uses them so that no platform math library
//! is on its path.

// On the JPEG AI exact decode path: no platform math-library call (see clippy.toml).
#![deny(clippy::disallowed_methods)]

/// 2⁵²: every `f64` of at least this magnitude is an integer.
const TWO_52: f64 = 4_503_599_627_370_496.0;
/// 2²³: every `f32` of at least this magnitude is an integer.
const TWO_23: f32 = 8_388_608.0;

/// `x` rounded to the nearest integer, ties away from zero: the value `f64::round`
/// returns, including the sign of a zero result.
pub fn round_half_away(x: f64) -> f64 {
    if x.is_nan() || x.abs() >= TWO_52 {
        return x; // NaN, an infinity, or already an integer
    }
    let t = x as i64; // toward zero; exact, since |x| < 2^52
    let d = x - t as f64; // exact (Sterbenz): the fraction, in (−1, 1)
    let r = if d >= 0.5 { t + 1 } else if d <= -0.5 { t - 1 } else { t };
    (r as f64).copysign(x)
}

/// `x` rounded to the nearest integer, ties to even: the value `f64::round_ties_even`
/// returns, including the sign of a zero result.
pub fn round_half_even(x: f64) -> f64 {
    if x.is_nan() || x.abs() >= TWO_52 {
        return x;
    }
    let t = x as i64;
    let d = x - t as f64;
    let odd = t % 2 != 0;
    let r = if d > 0.5 || (d == 0.5 && odd) {
        t + 1
    } else if d < -0.5 || (d == -0.5 && odd) {
        t - 1
    } else {
        t
    };
    (r as f64).copysign(x)
}

/// [`round_half_away`] for `f32`: the value `f32::round` returns.
pub fn round_half_away_f32(x: f32) -> f32 {
    if x.is_nan() || x.abs() >= TWO_23 {
        return x;
    }
    let t = x as i32; // exact, since |x| < 2^23
    let d = x - t as f32; // exact: the fraction, in (−1, 1)
    let r = if d >= 0.5 { t + 1 } else if d <= -0.5 { t - 1 } else { t };
    (r as f32).copysign(x)
}

/// `2^e` for `−1022 ≤ e ≤ 1023`, exact, from its exponent bits.
fn pow2(e: i32) -> f64 {
    debug_assert!((-1022..=1023).contains(&e));
    f64::from_bits(((e + 1023) as u64) << 52)
}

/// ln 2 split so that `k · LN2_HI` is exact for |k| < 2²¹ (its low 21 significand bits
/// are zero); `LN2_HI + LN2_LO` is ln 2 to about 10⁻²⁶.
const LN2_HI: f64 = f64::from_bits(0x3FE6_2E42_FEE0_0000);
const LN2_LO: f64 = f64::from_bits(0x3DEA_39EF_3579_3C76);

/// `1/n!` for `n = 1..=13`, each correctly rounded by the compiler.
const INV_FACT: [f64; 13] = [
    1.0,
    0.5,
    0.16666666666666666,
    0.041666666666666664,
    0.008333333333333333,
    0.001388888888888889,
    0.0001984126984126984,
    2.48015873015873e-05,
    2.7557319223985893e-06,
    2.755731922398589e-07,
    2.505210838544172e-08,
    2.08767569878681e-09,
    1.6059043836821613e-10,
];

/// `e^x`, from basic operations only.
///
/// Range reduction `x = k·ln2 + r` with `|r| ≤ ln2/2` (the first subtraction is exact,
/// by the split of ln 2), then `e^r = 1 + r·p(r)` with `p` the degree-12 Taylor
/// polynomial of `(e^r − 1)/r` in Horner form, then an exact scaling by `2^k`. On the
/// reduced range, truncating the series leaves the polynomial within 1.4·10⁻¹⁷
/// relative of `(e^r − 1)/r`, and `1 + r·p(r)` within 5.8·10⁻¹⁸ relative of `e^r` (both
/// at `r = −ln2/2`), far below one rounding (1.1·10⁻¹⁶). So the result's error is a
/// few rounding errors; a test bounds it against the host math library. It is not
/// correctly rounded, and it need not be: what it guarantees is the same bits on every
/// target.
pub fn exp(x: f64) -> f64 {
    if x.is_nan() {
        return x;
    }
    if x > 709.782_712_893_384 {
        return f64::INFINITY; // above ln(f64::MAX)
    }
    if x < -745.133_219_101_941_2 {
        return 0.0; // below ln of half the least subnormal
    }
    let k = round_half_away(x * std::f64::consts::LOG2_E); // x/ln2; |k| ≤ 1075
    let r = (x - k * LN2_HI) - k * LN2_LO;
    let mut p = INV_FACT[12];
    for &c in INV_FACT[..12].iter().rev() {
        p = p * r + c;
    }
    let er = 1.0 + r * p; // in [0.70, 1.42]
    let k = k as i32;
    if k > 1023 {
        er * 2.0 * pow2(k - 1) // k ≤ 1024: one rounding, which may overflow to ∞
    } else if k < -1022 {
        er * pow2(k + 1000) * pow2(-1000) // exact, then one rounding into the subnormals
    } else {
        er * pow2(k)
    }
}

#[cfg(test)]
#[allow(clippy::disallowed_methods, reason = "the tests compare with the standard library's functions")]
mod tests {
    use super::*;

    /// Deterministic sample of f64 bit patterns (an LCG), for sweeps.
    fn lcg(seed: u64) -> impl Iterator<Item = u64> {
        let mut s = seed;
        std::iter::from_fn(move || {
            s = s.wrapping_mul(6_364_136_223_846_793_005).wrapping_add(1_442_695_040_888_963_407);
            Some(s)
        })
    }

    /// The distance in units in the last place between two finite f64 of one sign.
    fn ulps(a: f64, b: f64) -> u64 {
        (a.to_bits() as i64 - b.to_bits() as i64).unsigned_abs()
    }

    fn same(a: f64, b: f64) -> bool {
        a.to_bits() == b.to_bits() || (a.is_nan() && b.is_nan())
    }

    fn edge_values() -> Vec<f64> {
        let mut v = vec![0.0, -0.0, 0.5, -0.5, 1.5, -1.5, 2.5, -2.5, 0.49999999999999994,
            -0.49999999999999994, 4503599627370495.5, -4503599627370495.5, 2251799813685247.5,
            TWO_52, -TWO_52, f64::MAX, f64::MIN, f64::MIN_POSITIVE, f64::from_bits(1),
            -f64::from_bits(1), f64::INFINITY, f64::NEG_INFINITY, f64::NAN];
        for i in -4000..=4000 {
            v.push(i as f64 / 8.0);
        }
        v
    }

    #[test]
    fn rounding_matches_the_exactly_specified_std_functions() {
        let mut xs = edge_values();
        xs.extend(lcg(7).take(200_000).map(f64::from_bits));
        // and values with a small fraction part, where the ties are
        xs.extend(lcg(11).take(200_000).map(|b| (b >> 11) as f64 / 64.0 - 70_368_744_177_664.0));
        for &x in &xs {
            assert!(same(round_half_away(x), x.round()), "round_half_away({x:e})");
            assert!(same(round_half_even(x), x.round_ties_even()), "round_half_even({x:e})");
            let xf = x as f32;
            let (a, b) = (round_half_away_f32(xf), xf.round());
            assert!(a.to_bits() == b.to_bits() || (a.is_nan() && b.is_nan()), "round_half_away_f32({xf:e})");
        }
    }

    /// Accuracy against the host math library, over every input the JPEG AI gain
    /// scaler can see (`n · GAIN_LOG_K / 128` for |n| ≤ 8192, which covers every
    /// gain-table entry with every 12-bit β) and a sweep of the whole finite range. This
    /// bounds the error; it is not what makes `exp` portable (that is its construction,
    /// pinned below).
    ///
    /// The assertion is about the host's `exp` as much as about this one. The gate below
    /// selects by target: the test runs on macOS arm64, on Linux arm64 and x86_64 with
    /// glibc, and on WASI, and is ignored, with the reason, everywhere else, Windows
    /// included; the pinned test below still runs there. A target is not a host. The
    /// gate cannot see the C library's version or the processor, so the test also runs
    /// against glibc versions and x86_64 processors that were not measured. Which hosts
    /// were measured, and the worst case on each, is in `extensions/determinism-tiers.md`
    /// §6. A failure on a host not listed there says that its `exp` and this one are more
    /// than 2 ULPs apart, not that this module's bits changed; the pinned test checks
    /// those.
    #[test]
    #[cfg_attr(not(any(
        all(target_os = "macos", target_arch = "aarch64"),
        all(target_os = "linux", target_env = "gnu", any(target_arch = "x86_64", target_arch = "aarch64")),
        all(target_arch = "wasm32", target_os = "wasi"),
    )), ignore = "the host math library's exp has not been measured on this platform; \
                  float_repro's pinned-output test covers portability here")]
    fn exp_is_within_two_ulps_of_the_host_math_library() {
        let mut worst = (0u64, 0.0f64);
        let mut check = |x: f64| {
            let (a, b) = (exp(x), x.exp());
            if b == 0.0 || b.is_infinite() || b < f64::MIN_POSITIVE {
                assert!(a == b || (b < f64::MIN_POSITIVE && ulps(a, b) <= 2), "exp({x:e}) = {a:e}, host {b:e}");
                return;
            }
            let u = ulps(a, b);
            if u > worst.0 {
                worst = (u, x);
            }
        };
        for n in -8192i64..=8192 {
            check(n as f64 * crate::jpegai_pipeline::GAIN_LOG_K / 128.0);
        }
        for i in -48_000i64..=45_440 {
            check(i as f64 / 64.0);
        }
        for b in lcg(3).take(300_000) {
            check((b >> 11) as f64 / (1u64 << 53) as f64 * 1454.9 - 745.1);
        }
        eprintln!("float_repro::exp vs host exp: worst {} ULP(s), at x = {:e}", worst.0, worst.1);
        assert!(worst.0 <= 2, "exp differs from the host by {} ULPs at x = {:e}", worst.0, worst.1);
        assert!(same(exp(0.0), 1.0) && same(exp(f64::NAN), f64::NAN));
        assert_eq!(exp(710.0), f64::INFINITY);
        assert_eq!(exp(-746.0), 0.0);
    }

    /// The portability claim itself: these functions' bits over a fixed grid, hashed
    /// and pinned. It needs no model or input file, so it runs wherever the library
    /// tests run: the conformance workflow's library-test step runs it on each leg
    /// (Linux x86_64, Linux arm64, macOS arm64, Windows x86_64), and a target that
    /// computes any of these functions differently fails here.
    #[test]
    fn outputs_are_pinned_bit_for_bit() {
        let mut h = blake3::Hasher::new();
        // A fixed grid, not the gain scaler's inputs: the factor is the gain unit's
        // `log_k` literal from before its correction, kept so the pinned hash still
        // names the same computation.
        for n in -8192i64..=8192 {
            h.update(&exp(n as f64 * 0.200_365_493_161_539_4 / 128.0).to_bits().to_le_bytes());
        }
        for i in -48_000i64..=45_440 {
            h.update(&exp(i as f64 / 64.0).to_bits().to_le_bytes());
        }
        for b in lcg(3).take(100_000) {
            h.update(&exp((b >> 11) as f64 / (1u64 << 53) as f64 * 1454.9 - 745.1).to_bits().to_le_bytes());
        }
        // NaN is left out: its payload through a conversion is not fixed on every target.
        for &x in edge_values().iter().filter(|x| !x.is_nan()) {
            h.update(&round_half_away(x).to_bits().to_le_bytes());
            h.update(&round_half_even(x).to_bits().to_le_bytes());
            h.update(&round_half_away_f32(x as f32).to_bits().to_le_bytes());
        }
        assert_eq!(h.finalize().to_hex().as_str(),
            "80174a4b00930982e9e2ea4ffdc47805455d4cc59797b41eb562a1cca55b0cfa",
            "float_repro output changed");
    }
}