alice_det_math/

lib.rs

//! Deterministic `f32` / `f64` transcendentals (cross-platform bit-exact).
//!
//! IEEE 754 guarantees that `+`, `-`, `*`, `/` and `sqrt` produce the same
//! bits on every platform Rust targets (SSE2+, aarch64, wasm32). It does
//! **not** cover `sin`, `exp`, `ln`, `powf`, `cbrt`, `hypot`, … — those come
//! from the platform `libm` and differ in the last ulp between macOS, glibc,
//! MSVC and wasm. A lockstep simulation that calls them diverges across
//! peers; a signed-distance field evaluated on two machines disagrees about
//! which side of a surface a point is on.
//!
//! Every function here is built from the guaranteed operations only —
//! integer range reduction (`trunc(t + (t < 0 ? -0.5 : 0.5))`, three
//! operations), fixed-degree polynomials evaluated in a fixed order, and
//! bit-level exponent construction — so the result is a pure function of
//! the input bits on every target. `mul_add` is not used
//! anywhere: it fuses into one rounding on FMA hardware and stays two
//! roundings elsewhere, which would break the guarantee.
//!
//! The crate is `no_std` (no `alloc` either). `round` and `sqrt`, which
//! `core` does not provide, are in [`ops`]; with the `std` feature `sqrt` is
//! the hardware instruction, without it a correctly rounded software one —
//! both return the same bits.
//!
//! With the `simd` feature, the `simd` module provides `wide::f32x8` versions
//! that are lane-for-lane bit-identical to the scalar functions.
//!
//! ```rust
//! use alice_det_math::{atan2, exp, sin_cos, sin_cos64};
//!
//! let (s, c) = sin_cos(1.0_f32);
//! assert_eq!(s.to_bits(), 0x3f57_6aa5); // the same bits on every platform
//! assert_eq!(c.to_bits(), 0x3f0a_5140);
//! assert_eq!(exp(1.0_f32).to_bits(), 0x402d_f854);
//! let a = atan2(0.0_f32, -1.0); // π, fdlibm special cases
//! assert_eq!(a, core::f32::consts::PI);
//!
//! // f64 over the whole range: 1e22 needs the Payne–Hanek reduction
//! let (s, c) = sin_cos64(1.0e22);
//! assert_eq!(s.to_bits(), 0xbfeb_453a_b76b_f397); // -0.8522008497671888
//! assert_eq!(c.to_bits(), 0x3fe0_be2c_ef01_c8f4); // 0.523214785395139
//! ```
//!
//! # Accuracy
//!
//! Measured over dense sweeps (`tests/accuracy.rs`). `f32` functions are
//! compared with a correctly rounded reference (`f64` libm rounded once to
//! `f32`), [`sin64`] / [`cos64`] with an independent correctly rounded
//! reference, the other `f64` functions with the platform `f64` libm — which itself
//! varies by ≤ 1 ulp between macOS / glibc / MSVC, so those bounds carry that
//! slack:
//!
//! | function | domain measured | max error |
//! |----------|-----------------|-----------|
//! | [`sin`] / [`cos`] / [`sin_cos`] | `\|x\| ≤ 100` | ≤ 2 ulp (Cephes single-precision coefficients), including every input within 64 ulp of a zero `k·π/2`; measured max 2 out to `\|x\| ≤ 8192` (`tests/accuracy_trig_zeros.rs`) |
//! | [`exp`] | `[-87, 88]` | ≤ 2 ulp |
//! | [`ln`] | `[1e-30, 1e30]` | ≤ 1 ulp (musl `logf` algorithm) |
//! | [`log2`] / [`log10`] | `[1e-30, 1e30]` | ≤ 1 ulp, measured max 0 ([`log2_64`] / [`log10_64`] rounded once; powers of two are exact) |
//! | [`cbrt`] | `[1e-30, 1e30]` | ≤ 1 ulp (bit-hack seed + 3 Newton steps) |
//! | [`hypot`] | finite | ≤ 1 ulp of `sqrt(x² + y²)` |
//! | [`powf`] | `x ∈ [1e-3, 1e3]`, `\|y\| ≤ 8` | ≤ 1 ulp (evaluated in `f64`, rounded once) |
//! | [`atan`] / [`atan2`] | `\|x\| ≤ 1e4`, all quadrants | ≤ 1 ulp (fdlibm `s_atanf` / `e_atan2f`, single precision) |
//! | [`asin`] / [`acos`] | `[-1, 1]` | ≤ 1 ulp (fdlibm in `f64`, rounded once) |
//! | [`tan`] | `\|x\| ≤ 100` | ≤ 1 ulp (fdlibm `k_sin`/`k_cos` in `f64`) |
//! | [`tanh`] | finite | ≤ 1 ulp (`exp64` based, `x` below 2^-14) |
//! | [`sin64`] / [`cos64`] / [`sin_cos64`] | every finite `f64` (Cody–Waite below `2^20·π/2`, Payne–Hanek above) | ≤ 1 ulp of correctly rounded, measured max 1 over 6234 points of an independent 2400-bit reference (musl `sin.c` / `cos.c`) |
//! | [`tan64`] | every finite `f64` (same reduction as [`sin64`]) | ≤ 1 ulp of correctly rounded, measured max 1 over 6310 points of an independent 2400-bit reference (fdlibm `k_tan.c`) |
//! | [`atan64`] / [`atan2_64`] | every finite `f64`, all quadrants | ≤ 1 ulp of correctly rounded, measured max 1 over 1758 points of an independent 2400-bit reference covering every branch of the reduction (fdlibm `s_atan.c` / `e_atan2.c`) |
//! | [`exp64`] / [`ln64`] | `[-745, 710]` / `[1e-300, 1e300]` | ≤ 2 ulp, measured max 1 against an independent 2400-bit reference (fdlibm algorithms) |
//! | [`log2_64`] / [`log10_64`] | `[1e-300, 1e300]` | ≤ 2 ulp of the platform libm, measured max 2 / 1 (`ln64` converted to the base; powers of two are exact for [`log2_64`]) |
//! | [`powf64`] | `x ∈ [1e-3, 1e3]`, `\|y\| ≤ 8` | ≤ 16 ulp measured 13 (`exp64(y·ln64 x)` with double-double argument) ⚠️ outside this domain the result can be *wrong*, not just less accurate — see [`powf64`] |
//! | [`round`] / [`round64`] / [`sqrt`] / [`sqrt64`] | all | exact (bit-identical to `f32::round` / `f32::sqrt`, except that every NaN is the canonical one) |
//!
//! The reduction subtracts `k·π/2` in four `f32` pieces, the third formed
//! exactly with Dekker's two-product (no fused multiply-add, so every target
//! computes the same bits). The first two products are exact only for
//! `|k| < 2¹³`, so large arguments to [`sin`] / [`cos`] (`|x| > 2¹³`) are
//! still deterministic but lose accuracy, exactly as `libm`'s own `sinf` does
//! without Payne–Hanek. Within the domain the remaining error is the Cephes
//! polynomial's (≤ 2 ulp). [`sin64`] / [`cos64`] reduce
//! every finite argument exactly (Payne–Hanek with 2/π to 1584 bits), so
//! their bound holds up to `f64::MAX`; the reference values are in
//! `tests/data/sin_cos64_reference.txt`, generated by
//! `scripts/gen_sin_cos64_reference.py` (`mpmath`, not a libm). Every other
//! `f64` function whose bound is asserted has its own table under
//! `tests/data/`, from the matching `scripts/gen_*_reference.py`, so no `f64`
//! bound depends on a platform `libm`. [`tan64`] is
//! deliberately not compared with a platform `libm` at all: one disagrees
//! with the correctly rounded value by ~1.0e5 ulp in the Payne–Hanek range,
//! so a bound on that difference would depend on the machine running it.
//!
//! # Guarantee scope
//!
//! Bit-exactness relies on the target honouring IEEE 754 for the basic
//! operations: `x86_64` (SSE2), `aarch64`, `wasm32` and every other Rust tier-1/2
//! target qualify. 32-bit x86 built for x87 (`i586`) without SSE2, or any
//! build with fast-math style flags, is outside the guarantee.
//! `tests/golden.rs` pins a SHA-256 of every function's output over a fixed
//! input grid; CI runs it on macOS (ARM / Intel), Linux (x86 / ARM), Windows
//! and wasm32, and with the SIMD fallback paths forced.
//!
//! # Policy for consumers
//!
//! List the `f32` / `f64` `libm` methods under `clippy.toml`
//! `disallowed-methods` and run clippy with `-D warnings`, so a stray
//! `x.sin()` on a float is a red gate; call these functions instead.

#![cfg_attr(not(feature = "std"), no_std)]
#![cfg_attr(docsrs, feature(doc_cfg))]

/// Identifies the numeric behaviour of this crate: it changes whenever any
/// function here would return different bits for the same input, and
/// otherwise only when the set of functions it covers grows.
///
/// ```text
/// 135c407d4249c1c200e4d1be79a12e848ce5c355fc3a918da965849e697ef000
/// ```
///
/// It is the SHA-256 of the per-function bit pins in `tests/golden.rs` —
/// each of which is itself a SHA-256 over one function's output on a fixed
/// input grid — folded in ascending order of function name, each as the
/// name's length in four big-endian bytes, then the name, then its 32-byte
/// pin. Because the inputs are the recorded pins and nothing else, the value
/// is a property of this crate alone: no platform `libm`, no timing, no
/// environment enters it, and it is the same on every target the pins
/// reproduce on.
///
/// A caller that derives an identifier from a formula plus its parameters can
/// mix this in, so that the identifier also covers *how the arithmetic was
/// evaluated*. Two runs that agree on formula, parameters and this value
/// computed the same bits; if this value differs, they did not, even when
/// everything else matches.
///
/// `tests/golden.rs` recomputes it from the pins and fails if the two
/// disagree; it also reads every `pub fn` this crate defines out of its own
/// sources and fails if one has no pin, so the constant cannot drift from the
/// behaviour it names. The SIMD kernels are covered by parity with the scalar
/// functions rather than by pins of their own: they are feature-gated, and an
/// identifier that changed with the feature set would not identify the
/// arithmetic.
///
/// # Re-recording it, in this order
///
/// After an intentional change to any function's output, update the
/// per-function hashes in `tests/golden.rs` **first**, and only then recompute
/// this value. The fold reads the table as it is written in the source, so
/// running the re-pin before the table is updated reproduces the old value and
/// looks as though the identifier failed to notice the change.
pub const SEMANTICS_ID: [u8; 32] = [
    0x13, 0x5c, 0x40, 0x7d, 0x42, 0x49, 0xc1, 0xc2, 0x00, 0xe4, 0xd1, 0xbe, 0x79, 0xa1, 0x2e, 0x84,
    0x8c, 0xe5, 0xc3, 0x55, 0xfc, 0x3a, 0x91, 0x8d, 0xa9, 0x65, 0x84, 0x9e, 0x69, 0x7e, 0xf0, 0x00,
];

pub mod semantics;
pub use semantics::{check_semantics, SemanticsCheck, PREVIOUS_SEMANTICS_IDS};

pub mod double;
pub mod metric;
pub mod ops;
#[cfg(feature = "simd")]
#[cfg_attr(docsrs, doc(cfg(feature = "simd")))]
pub mod simd;
pub mod single;

/// The `f32` kernel constants (Cephes / musl), for code generators that must
/// emit the *same* law elsewhere (a JIT, a shader) and be checked against
/// these functions bit for bit.
pub mod consts {
    pub use crate::single::{
        COS_P, EXP_HI, EXP_LO, EXP_P, FRAC_2_PI, LG1, LG2, LG3, LG4, LN2_HI, LN2_LO, LOG2E, PIO2_1,
        PIO2_2, PIO2_3, PIO2_4, SIN_P, SPLIT_MASK, SQRT2_BITS,
    };
}

pub use double::{
    acos64, asin64, atan2_64, atan64, cos64, exp64, ln64, log10_64, log2_64, powf64, sin64,
    sin_cos64, tan64,
};
pub use ops::{round, round64, sqrt, sqrt64};
pub use single::{
    acos, asin, atan, atan2, cbrt, cos, exp, hypot, ln, log10, log2, powf, powi, sin, sin_cos, tan,
    tanh,
};