1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
//! 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.
/// 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: = ;
pub use ;
/// 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 use ;
pub use ;
pub use ;