diff --git a/.claude/blackboard.md b/.claude/blackboard.md index a27b987b..deff2874 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -1,3 +1,63 @@ +## 2026-10-27 — wasm `mul_add` certification audit (#378): relaxed-SIMD, proof obligations, X11 + +- **Relaxed SIMD, tested:** no source path selects on `relaxed-simd` (only the `simd_wasm.rs` module note mentions it). The parity wasm built with `+simd128` and with `+simd128,+relaxed-simd` has a byte-identical code section (621,383 bytes); only the `name` and `target_features` custom sections differ. Byte scan of the code section for relaxed-SIMD opcodes (`0xFD`, sub-opcode `0x100..0x113`): 0 hits; a control calling `f64x2_relaxed_madd` scans to 1 hit. The relaxed build passes under node 22. New arm `scripts/masking-parity.sh wasm-relaxed` and a CI step in `realization/wasm + scalar × wasm32`, as a guard against a relaxed path appearing later. What it establishes: the shipped arithmetic is the same under both flags for this compiler; it does not test an engine executing relaxed instructions, since none are emitted. +- **Proof obligations:** `.claude/knowledge/wasm-fma-proof-obligations.md` separates PROVEN (F32 path end to end; F64 no-overflow, Dekker condition `e_a + e_b >= -970` vs the guard's `-920`, TwoSum, the `tl != 0 ⇒ |th| >= 2^-922` lemma, round-to-odd through the midpoint argument incl. subnormal `v`, near-zero cancellation, the exact-product repair, bit-select, mixed pairs) from CITED (Veltkamp split, Dekker product, Knuth TwoSum), TESTED (compiler output per realization; the integer fallback `simd_fma::fma_f64` against the oracle) and OPEN (no mechanized proof). The earlier "~37 binades of margin per operand" is a measurement of where the unguarded algorithm breaks, not part of the proof; the proven margin is 50 binades of product exponent. Guard unchanged. +- **Bit-select audit:** masks are full 64-bit lanes from `f64x2` compares; both candidates are lane-wise; wasm has no FP exceptions, so discarded-lane NaNs cannot reach a kept lane. No code change. +- **X11 rewritten** in `coresearch-2026-10-11-simd-fleet.md`: current status first. The superseded wording, kept here as history: "`F32x16::mul_add` is not bit-identical across realizations: wasm without `relaxed-simd` is mul-then-add, with it the host chooses; x86 fuses. lance-graph D-WFL-DET must exclude it". +- **lance-graph (untouched):** D-WFL-DET (replay determinism) may rely on a verified single-rounding `mul_add` across realizations. D-OIF-4 defines its kernel as multiply then add and must not silently become a fused FMA. Different source-level semantics; neither changes. +- **Haiku:** `census-counter` ×1 (cross-lane ops in the FMA kernels and dispatchers of `simd_wasm.rs`). Re-ran 2 of its counts (F64 dispatcher `extract_lane` = 6; kernels, any cross-lane op = 0): 0 wrong. + +## 2026-10-27 — wasm `mul_add`: certified domain, exact-product repair, the open item closed + +- **Closes:** the 2026-10-11 open item "`mul_add` not bit-identical across realizations (wasm)" and coresearch row X11. Both were written before PR #374 (merged 2026-10-24, merge a78de7b). #374 made wasm `F32x16::mul_add` and `F64x8::mul_add` round once. No blackboard entry recorded it. +- **Fast domains:** + - `F32x16`: no guard. Every f32 input is in domain. The product of two f32 is exact in f64, and the f64 sum is rounded to odd before the single rounding to f32 (Boldo & Melquiond, IEEE TC 2008). The f64 sum cannot overflow for f32 operands. + - `F64x8` vector path: Dekker product + two-sum + round-to-odd, only when 2^-460 <= |a|,|b| <= 2^460 and |c| <= 2^960. Decided per v128 pair (2 lanes). +- **Change (path B, new local repair):** a lane whose a or b is a zero, an infinity or NaN has an exact product, so `a * b + c` rounds once. Such pairs now stay in the vector (per-lane select). This is the identity `simd_fma::fma_f64` starts with. Other pairs go lane by lane to the integer FMA `simd_fma::fma_f64` (path C, exact, no libm). +- **Guard margin, measured here:** scalar copy of the emulation vs hardware FMA on the x86 host, 200,000 samples per point. + - Failures start at products near 2^-1000 (operands near 2^-497). + - None for |c| up to f64::MAX, given the a,b bounds. + - All fail at a ~ 2^1000 (the split overflows). + - So the bounds have ~37 binades of margin per operand. The c bound is not needed for exactness; it excludes infinite/NaN c. +- **Independent reference:** the parity crate's 4608-bit integer oracle `oracle_fma`. It shares no arithmetic with any backend and calls no `mul_add`. It is itself checked against hardware FMA in `oracle_matches_hardware_fma`. +- **Gap found:** the 0x11xx group reached the F64 lane fallback only through random bit patterns. Widening the lower bound to 2^-600 stayed green. +- **Tests added:** round 5 (`edge_lanes_f64`) puts lanes on both sides of every bound, |c| == 2^960, mixed pairs, products below the subnormal range and above f64::MAX, and zeros/infinities/NaN in a or b against every kind of c. New anti-vacuity codes: + - 0x11F2: >= 256 lanes where the unguarded emulation (a scalar copy) is wrong. + - 0x11F3: >= 256 vector-path pairs. + - 0x11F4: >= 256 mixed pairs. + - 0x11F5: >= 128 repaired pairs. +- **Disable runs on the wasm arm (each restored), red:** + - guard removed (0x1159) + - lower bound 2^-520 (0x116c, the new edge round) + - a,b upper bound 2^530 (0x115a) + - per-pair `all_true` -> `any_true` (0x115c) + - f32 round-to-odd dropped (0x112e) + - round-to-odd step always up (0x112d) + - f64 round-to-odd dropped (0x1150) + - repair returns `c` (0x116e) + - repair drops `c` (0x116e) + - tiny a treated as exact product (0x1161) + - large or non-finite c treated as exact product (0x116e) +- **Disable runs, green, explained:** + - lower bound 2^-500: inside the measured margin, within ~3 binades of the onset. + - c bound 2^1023: no counterexample exists given the a,b bounds. + - Three narrowing disables of the repair predicate (zero not counted, both operands required, only a checked) stay green. The narrowed lanes fall back to the exact integer FMA. These are performance regressions, not correctness errors. +- **Performance on V8** (node 22, wasm32+simd128, best of 7, noise about +-30%): + - F64 vector path ~6.1-6.3 ns/element against ~1.1 ns unfused mul+add. Lane fallback ~50-58 ns/element. + - F32 ~3.1 ns against ~0.33 ns. + - Exact zeros before the repair: 10% zeros in a -> 19% of pairs on the lane path, ~9-13 ns/element. 50% zeros (ReLU-like) -> 75% of pairs, ~28-33 ns/element. + - After the repair: ~7.0 and ~9.3 ns/element. The in-domain path is unchanged. + - Random bit patterns: 95.8% of pairs out of domain. +- **Arms green after the change:** wasm, wasm-scalar, native (x86 AVX-512 host, hardware FMA), and `scripts/wasm-parity.sh`. Parity crate clippy `-D warnings` clean. +- **Not claimed:** + - NaN payload bits (any NaN matches any NaN). + - `reduce_sum`: a separate contract, not FMA. Summation order differs per backend; last-ulp differences are documented. + - Exhaustive coverage of the f32/f64 input space. Fixtures are random and adversarial only. +- **Open:** + - The wasm build of ndarray emits 3 dead-code warnings from `cpu_guard.rs` (`guard_before_main`, `guard_is_in_main_executable`, `contains`). Pre-existing and unrelated. No CI row runs clippy on wasm32. + - The narrowing direction of the repair predicate has no automated check. Correctness cannot see it; only the bench does. +- **Haiku:** `harvest-scout` ×1 (`F64x8`/`F32x16` `mul_add` call sites, 141 hits). Checked: `heel_f64x8.rs:73`, `signature_pde.rs:143` read in place; counts re-run for `heel_f64x8.rs` (8), `signature_pde.rs` (4), `reductions.rs` (5). 0 wrong, misses not measured. + ## 2026-10-14 — Haiku cards for supervised harvest and mapping - Five Haiku cards added (`harvest-scout`, `census-counter`, `quote-harvester`, `substrate-mapper`, `spec-editor`); rules in `.claude/rules/agent-model-tiers.md`, CLAUDE.md Agent Protocol item 7. No `guarded-executor`: only the orchestrator compiles here. diff --git a/.claude/knowledge/coresearch-2026-10-11-simd-fleet.md b/.claude/knowledge/coresearch-2026-10-11-simd-fleet.md index d6dea0e9..c0912469 100644 --- a/.claude/knowledge/coresearch-2026-10-11-simd-fleet.md +++ b/.claude/knowledge/coresearch-2026-10-11-simd-fleet.md @@ -88,7 +88,7 @@ non-FMA tiers stay LOOSE (correct, slow; no consumer on those CPUs named). | X3/X4 | ladybug's multi-slice image picks by `grep avx512f /proc/cpuinfo`; v4 also needs bw/cd/dq/vl. The exact selector already exists: the guard's own predicate | Q-HOST | | X6 | MedCare-rs and q2 ship `x86-64-v4` images to Railway; Railway publishes no CPU guarantee; GCP E2 mixes Broadwell (no AVX-512); MedCare's own handover says "v4 build never run"; MedCare measured the v4 gain as µs beside a 43 ms KDF | Q-HOST | | X7 | psABI: v4 = AVX512F/BW/CD/DQ/VL only; VNNI/BF16/AMX/VPOPCNTDQ are outside every level, so `config-avx512.toml` and holograph's `avx512vpopcntdq` are "beyond v4" | Q-TIER | -| X11 | `F32x16::mul_add` is not bit-identical across realizations: wasm without `relaxed-simd` is mul-then-add, with it the host chooses; x86 fuses. lance-graph D-WFL-DET must exclude it | Q-TIER | +| X11 | **CLOSED (#374, #378).** Wasm `F32x16::mul_add` and `F64x8::mul_add` round once, like every other realization. Each realization is checked bit for bit against the parity crate's integer oracle (any NaN matches any NaN), and `+relaxed-simd` emits the same code section as `+simd128`. A determinism gate (lance-graph D-WFL-DET) need not exclude a single-rounding `mul_add` on cross-realization grounds; D-OIF-4's separate mul-then-add contract is unaffected. `reduce_sum` is not covered: its summation order differs per backend. The superseded observation (wasm mul-then-add, host-chosen under `relaxed-simd`) is kept in the blackboard entries of 2026-10-27. | Q-TIER | | X16 | lance-graph `ndarray_bridge.rs` intrinsic blocks are all under `cfg(not(feature = "ndarray-hpc"))` and `ndarray-hpc` is default → dead in default builds; holograph `hamming.rs` is live (feature `simd` default, no ndarray dep); ladybug `collapse_gate.rs` is a per-call `is_x86_feature_detected!` + raw `_mm256_*` kernel in no sweep ledger | Q-SWEEP | | X17 | the MSRV issue still lists a `rev`-pin option forbidden by `E-CONSUMER-PINS-ON-INTERNAL-SIBLINGS-PROHIBITED-1`; `CARGO-COMPUTE-SUBSTRATE.md` §3 still shows `rev = ""`; lichess-rs says `cargo +1.95`; tesseract `rust.yml:59` says "1.95 floor" | Q-MSRV | | X18 | OGAR pulls `encryption` (+ argon2 `ndarray-simd`) from ndarray **git master** while lance-graph uses a path: two package ids in a binary linking both. surrealdb's `ndarray = "0.17.2"` is crates.io, a different package from the fork at the same version | Q-MSRV | diff --git a/.claude/knowledge/wasm-fma-proof-obligations.md b/.claude/knowledge/wasm-fma-proof-obligations.md new file mode 100644 index 00000000..24d8f158 --- /dev/null +++ b/.claude/knowledge/wasm-fma-proof-obligations.md @@ -0,0 +1,159 @@ +# wasm `mul_add`: proof obligations + +READ BY: anyone changing `F32x16::mul_add` / `F64x8::mul_add` in +`src/simd_wasm.rs`, the scalar emulation in `src/simd_fma.rs`, or the `0x11xx` +parity group in `crates/simd-masking-parity`. + +Contract: every lane returns `RN(a * b + c)`, binary32 or binary64, round to +nearest, ties to even, one rounding of the exact value, subnormals included. +Any NaN result matches any NaN: payload and NaN sign are not part of the +contract. + +The proofs below rest on one platform fact and three cited lemmas. Everything +else is derived here. Each item is marked **PROVEN** (derived here from the +lemmas), **CITED** (taken from the literature, not re-proven), **TESTED** +(executable witness only) or **OPEN**. + +## Platform fact + +WebAssembly scalar and SIMD float operations (`f64x2.add/sub/mul`, +`f32x4.demote_f64x2_zero`, `f64x2.promote_low_f32x4`) are IEEE 754 operations +with round to nearest, ties to even and gradual underflow. There is no flush to +zero, no trap and no status flag. Only NaN bit patterns are nondeterministic. +Rust emits no `contract` or reassociation flags, so `x * y + z` stays two +operations. + +`+relaxed-simd` does not change this code. No source path is selected on it, +and the code section of the parity binary is byte-identical with and without +it. A byte scan for relaxed-SIMD opcodes (`0xFD`, sub-opcode `0x100..0x113`) +finds 0 hits there, and 1 hit in a control that calls `f64x2_relaxed_madd` +(blackboard 2026-10-27, audit entry). **TESTED** for the current compiler; +`scripts/masking-parity.sh wasm-relaxed` re-checks it in CI. + +## Lemmas + +- **L1 (Veltkamp split)** — **CITED**. With `C = 2^27 + 1`, + `t = C*x; hi = t - (t - x); lo = x - hi` gives `hi + lo == x` exactly, `hi` and + `lo` each fit in 26 bits, provided `C*x` does not overflow. Underflow does not + break it. (Dekker 1971; Boldo, "Pitfalls of a full floating-point proof", + IJCAR 2006.) +- **L2 (Dekker product)** — **CITED**. With L1 splits of `a` and `b`, + `ph = RN(a*b)` and `pl = (((ah*bh - ph) + ah*bl) + al*bh) + al*bl` satisfy + `ph + pl == a*b` exactly, provided nothing overflows and + `e_a + e_b >= e_min + p - 1 = -970` (`e_x` the binary exponent of `x`, `x` in + `[2^e, 2^(e+1))`). (Same references.) +- **L3 (TwoSum)** — **CITED**. `s = x + y; yy = s - x; + e = (x - (s - yy)) + (y - yy)` gives `s + e == x + y` exactly for finite + `x, y` when nothing overflows, underflow included. (Knuth, TAOCP vol. 2.) +- **L4 (round to odd)** — **PROVEN** below. Let `w` be real and `v` the binary64 + value obtained by `add_round_to_odd`, at grid spacing `g` (the ulp at `v`, or + `2^-1074`). If `w` is representable, `v == w`. Otherwise `v` is the odd + multiple of `g` with `|w - v| < g`, so `w` and `v` lie strictly between the + same two consecutive even multiples of `g`. + + Proof. `add_round_to_odd(x, y)` takes `s = RN(x + y)` and the exact error `e` + (L3). If `e == 0`, `x + y == s` is representable. Otherwise, if `s` is odd it + is already the odd neighbour. If `s` is even, the code steps its bit pattern + one place toward `x + y` (down in magnitude when `e` and `s` differ in sign). + That gives the adjacent value on the side of the exact sum, which is odd; + across a power of two it is the largest value of the lower binade, which is + odd on the finer grid there. In every case the exact sum lies strictly + between `v - g` and `v + g`. +- **L5 (rounding through odd)** — **PROVEN**. Let `x = T + w`, where `T` lies + on a grid whose spacing is at least `8g`. Then + `RN_target(T + v) == RN_target(x)` for any target format whose rounding + boundaries (midpoints and the overflow threshold) near `x` are multiples of + `2g`. + + Proof. If `w` is representable, `T + v == x`. Otherwise every rounding + boundary `M` satisfies: `M - T` is an even multiple of `g`. By L4, `w` and + `v` lie strictly between the same two even multiples of `g`. So `x` and + `T + v` lie strictly on the same side of every boundary, and neither of them + is one. Hence both round to the same value. + +## F32x16 (no guard) + +- **Product exact** — **PROVEN**. A binary32 significand has 24 bits, so the + product has at most 48 bits (53 available). Its exponent is in `[-298, 256]`, + inside the binary64 normal range. +- **Sum** — **PROVEN**. `|a*b| <= 2^256` and `|c| < 2^128`, so the binary64 + sum cannot overflow and L3 applies. +- **Final rounding to binary32** — **PROVEN** by L5. Take `T = 0` and `w` the + exact sum. Every binary32 rounding boundary, including the subnormal + midpoints and the overflow threshold `2^128 - 2^103`, is a multiple of + `2^-150`. The binary64 grid at the value is at most `2^-52` times its + magnitude, which is at least 29 bits finer than binary32. That holds down to + the binary32 subnormal grid, since binary64 is still normal there. +- **Special values** — **PROVEN**. An infinite or NaN product, or an infinite + `c`, makes `s` non-finite. The step is then masked off, and `demote` returns + the IEEE infinity or NaN. Signed zeros: when `a*b + c` is exactly zero, + `e == 0`, so `s` is the IEEE-signed `RN` sum and is passed through unchanged. + +## F64x8 vector path + +Guard: `2^-460 <= |a|, |b| <= 2^460` and `|c| <= 2^960`, required on BOTH lanes +of a v128 pair. The guard excludes zero, infinite and NaN `a`, `b`, and +infinite or NaN `c`. + +- **No overflow anywhere** — **PROVEN**. `|C*x| < 2^488`; `|a*b| <= 2^920`; + `|c + ph| <= 2^960 + 2^920`. Every partial product in L2 is at most `|a*b|`. +- **L2 applies** — **PROVEN**. `e_a + e_b >= -920 >= -970`. This is 50 binades + of product exponent above what L2 needs (for example one operand at + `2^-460` and the other at `2^-510`). +- **TwoSum `(th, tl)` of `c + ph`** — **PROVEN** by L3 (finite, no overflow). +- **`tl != 0` implies `|th| >= 2^-922`** — **PROVEN**. If `c + ph` is inexact, + then `c` and `ph` are not opposite-signed within a factor of 2 (Sterbenz). + So `|c + ph| >= |ph| / 2 >= 2^-921`, and also `|ph| <= 2 |c + ph|`. +- **Final rounding** — **PROVEN**. Let `x = c + a*b = th + (tl + pl)`, and let + `v` be the round-to-odd of `w = tl + pl` (L4). + - Case `tl == 0`: `w = pl` is representable (L2), so `v == pl` and + `th + v == x` exactly. The final add then rounds once. This covers the + near-zero cancellations: an exact zero, and a subnormal or tiny result. + - Case `tl != 0`: we have `|tl| <= ulp(th)/2` and `|pl| <= ulp(ph)/2 <= 2 ulp(th)`. + So `|w| <= 2.5 ulp(th)`, and `g <= max(2^-50 ulp(th), 2^-1074) <= ulp(th)/8`, + because `ulp(th) >= 2^-974`. Binary64 midpoints near `x` are multiples of + `ulp(th)/4` (`ulp(th)/2`, or `ulp(th)/4` just below a power of two), so + they are multiples of `2g`. L5 with `T = th` gives `RN(th + v) == RN(x)`. + The result is at least `2^-922` in magnitude, so it is never subnormal. +- **Special-value repair** — **PROVEN**. If `a` or `b` is `±0`, `±inf` or NaN, + the binary64 product is exactly the IEEE product: `±0`, `±inf`, or NaN for + `0 * inf`. `a*b + c` is then one rounding of the exact sum, which is IEEE + `fusedMultiplyAdd`: invalid gives NaN, the zero-sum sign follows the + round-to-nearest rule, and a nonzero `c` is returned exactly. It is the same + identity `simd_fma::fma_f64` uses first. +- **Bit-select** — **PROVEN**. The masks come from `f64x2` comparisons, so each + 64-bit lane is all ones or all zeros. They stay so through `v128_and/or/not`. + `v128_bitselect` therefore returns, per lane, exactly one of the two computed + lane patterns. Both candidate computations are lane-wise: the kernels contain + no shuffle, swizzle or lane extract (Haiku census, verified). WebAssembly has + no float exceptions, so the NaN or garbage computed on a discarded lane has no + effect on the kept lane. +- **Mixed pairs** — **PROVEN**. The dispatcher picks one path per pair. The + vector path is entered only when every lane is in the guard domain or has an + exact product; all other pairs go lane by lane to `simd_fma::fma_f64`. + +## Margins: proven versus measured + +- **Proven sufficient conditions:** `e_a + e_b >= -970` for L2 (the guard gives + `-920`); `|C*x|` finite (`|x| < ~2^996`); `c + ph` not overflowing. Given the + guard on `a` and `b`, a finite `c` cannot make `c + ph` overflow, so the + `|c| <= 2^960` bound is not needed for exactness. It is kept because it + excludes infinite and NaN `c`, which TwoSum turns into NaN. +- **Measured** (scalar copy of the unguarded emulation vs hardware FMA, + 200,000 samples per point): failures start near products of `2^-1000` + (operands near `2^-497`), and every case fails at `|a| ~ 2^1000`. This is + evidence about where the algorithm breaks, not part of the proof. The guard + is not to be widened on its strength. + +## Still open + +- **OPEN:** no mechanized proof. L1–L3 are cited, not re-derived; L4, L5 and + the case analysis above are hand proofs. +- **OPEN:** `simd_fma::fma_f64` (the integer lane fallback) and + `simd_fma::fma_f32` are certified only by tests against the 4608-bit integer + oracle and against hardware FMA. `fma_f32` uses the same argument as the + F32x16 path above. +- **TESTED, not proven:** that each realization's compiler output matches the + source semantics. This rests on the parity arms (wasm, wasm-scalar, + wasm-relaxed, native, NEON, AVX, v2, nightly). Sampling covers no input space + exhaustively. diff --git a/.github/workflows/simd-matrix.yaml b/.github/workflows/simd-matrix.yaml index 7acea754..4d9a5ac8 100644 --- a/.github/workflows/simd-matrix.yaml +++ b/.github/workflows/simd-matrix.yaml @@ -381,6 +381,8 @@ jobs: run: bash scripts/masking-parity.sh wasm - name: masking parity (scalar realization = wasm32 without simd128) run: bash scripts/masking-parity.sh wasm-scalar + - name: masking parity (wasm, +simd128 +relaxed-simd) + run: bash scripts/masking-parity.sh wasm-relaxed nightly: # The `core::simd` realization behind the opt-in `nightly-simd` feature. diff --git a/crates/simd-masking-parity/src/lib.rs b/crates/simd-masking-parity/src/lib.rs index 30155a4e..b1d51969 100644 --- a/crates/simd-masking-parity/src/lib.rs +++ b/crates/simd-masking-parity/src/lib.rs @@ -2017,15 +2017,18 @@ fn check_mul_add_single_rounding() -> Result<(), u32> { // below the 53 bits after the product's error term pl, so `tl + pl` is // inexact and only its round-to-odd keeps the tie broken. Round 4 makes the // exact sum itself a midpoint, so the final rounding must pick the even - // neighbour (ties-away or ties-toward-zero fail it). + // neighbour (ties-away or ties-toward-zero fail it). Round 5 sits on the + // edges of the vector path's domain (see `edge_lanes_f64`). let mut unfused_differs = 0usize; + let (mut emulation_wrong, mut vector_pairs, mut mixed_pairs, mut exact_product_pairs) = + (0usize, 0usize, 0usize, 0usize); // (1 + 2^-26)(1 + 2^-27) = 1 + 2^-26 + 2^-27 + 2^-53: 2^-53 is half an ulp, a tie. let (ya, yb) = (1.0f64 + f64::EPSILON * 2f64.powi(26), 1.0f64 + f64::EPSILON * 2f64.powi(25)); - for round in 0..5120u32 { + for round in 0..6144u32 { let rin = |r: &mut SplitMix64| { f64::from_bits((r.next() & 0x800F_FFFF_FFFF_FFFF) | ((1023 - 400 + r.next() % 800) << 52)) }; - let (a, b, c): ([f64; 8], [f64; 8], [f64; 8]) = match round % 5 { + let (a, b, c): ([f64; 8], [f64; 8], [f64; 8]) = match round % 6 { 0 => ( core::array::from_fn(|_| rin(&mut rng)), core::array::from_fn(|_| rin(&mut rng)), @@ -2051,7 +2054,7 @@ fn check_mul_add_single_rounding() -> Result<(), u32> { core::array::from_fn(|_| f64::from_bits(rng.next())), core::array::from_fn(|_| f64::from_bits(rng.next())), ), - _ => { + 4 => { // a * b = (1 + 2^-26 + 2^-27 + 2^-53) * 2^k; each c cancels the middle // terms so the exact sum is a midpoint between two f64 neighbours: // lane 0: 1 + 2^-53 -> 1 (even below) @@ -2070,14 +2073,36 @@ fn check_mul_add_single_rounding() -> Result<(), u32> { core::array::from_fn(|i| neg(i) * cut[i % 4] * s), ) } + _ => edge_lanes_f64(&mut rng), }; + if round % 6 == 5 { + for p in 0..4 { + let ok = |i: usize| f64_vector_domain(a[i], b[i], c[i]); + match (ok(2 * p), ok(2 * p + 1)) { + (true, true) => vector_pairs += 1, + (false, false) => {} + _ => mixed_pairs += 1, + } + // In domain or with an exact product: the vector path's repair. + let exact = |i: usize| [a[i], b[i]].iter().any(|x| *x == 0.0 || !x.is_finite()); + let repaired = |i: usize| ok(i) || exact(i); + if repaired(2 * p) && repaired(2 * p + 1) && !(ok(2 * p) && ok(2 * p + 1)) { + exact_product_pairs += 1; + } + } + for i in 0..8 { + if !same64(emulated_fma_f64_unguarded(a[i], b[i], c[i]), oracle_fma_f64(a[i], b[i], c[i])) { + emulation_wrong += 1; + } + } + } let got = F64x8::from_array(a) .mul_add(F64x8::from_array(b), F64x8::from_array(c)) .to_array(); for i in 0..8 { let want = oracle_fma_f64(a[i], b[i], c[i]); if !same64(got[i], want) { - return Err(0x1140 + (round % 5) * 8 + i as u32); + return Err(0x1140 + (round % 6) * 8 + i as u32); } if !same64(a[i] * b[i] + c[i], want) { unfused_differs += 1; @@ -2087,9 +2112,191 @@ fn check_mul_add_single_rounding() -> Result<(), u32> { if unfused_differs < 256 { return Err(0x11F1); } + // Round 5 must reach both paths, split pairs between them, and contain + // lanes on which the emulation alone is wrong: those lanes fail if the + // domain guard of `F64x8::mul_add` is removed or widened. + if emulation_wrong < 256 { + return Err(0x11F2); + } + if vector_pairs < 256 { + return Err(0x11F3); + } + if mixed_pairs < 256 { + return Err(0x11F4); + } + if exact_product_pairs < 128 { + return Err(0x11F5); + } Ok(()) } +/// The domain of the wasm `F64x8::mul_add` vector path: +/// `2^-460 <= |a|, |b| <= 2^460` and `|c| <= 2^960` (NaN is outside). +fn f64_vector_domain(a: f64, b: f64, c: f64) -> bool { + let (lo, hi, chi) = (2f64.powi(-460), 2f64.powi(460), 2f64.powi(960)); + (lo..=hi).contains(&a.abs()) && (lo..=hi).contains(&b.abs()) && c.abs() <= chi +} + +/// The wasm vector path's emulated FMA without its domain guard, in scalar +/// f64: Dekker's product, two-sum with `c`, the error terms summed and rounded +/// to odd, one final rounding. Used only to show that the edge fixtures hold +/// lanes this algorithm gets wrong outside its domain. +fn emulated_fma_f64_unguarded(a: f64, b: f64, c: f64) -> f64 { + let split = |x: f64| { + let t = 134_217_729.0 * x; // 2^27 + 1 + let hi = t - (t - x); + (hi, x - hi) + }; + let two_sum = |x: f64, y: f64| { + let s = x + y; + let yy = s - x; + (s, (x - (s - yy)) + (y - yy)) + }; + let (ah, al) = split(a); + let (bh, bl) = split(b); + let ph = a * b; + let pl = (((ah * bh - ph) + ah * bl) + al * bh) + al * bl; + let (th, tl) = two_sum(c, ph); + let (s, e) = two_sum(tl, pl); + let bits = s.to_bits(); + let v = if e != 0.0 && s.is_finite() && bits & 1 == 0 { + f64::from_bits(if (bits ^ e.to_bits()) >> 63 == 1 { + bits - 1 + } else { + bits + 1 + }) + } else { + s + }; + th + v +} + +/// Lanes on and around the edges of [`f64_vector_domain`]. Pairs are the v128 +/// halves `F64x8::mul_add` dispatches on: +/// * lanes 0, 1: inside, at the lower and the upper bound of `|a|`, `|b|`; +/// * lanes 2, 3: one lane just below the lower bound (or `|c|` just above +/// `2^960`), one inside with `|c| == 2^960`, so the pair is mixed; +/// * lanes 4, 6: lanes the emulation gets wrong outside the domain (a product +/// whose error terms underflow; a product above `f64::MAX`), `|a|` just above +/// `2^460`, and a zero, infinity or NaN in `a` or `b` (the exact-product +/// repair, which keeps the pair in the vector) against every kind of `c`; +/// * lanes 5, 7: inside, so pairs (4, 5) and (6, 7) are mixed and a widened +/// guard sends the witness lanes to the vector path. +fn edge_lanes_f64(rng: &mut SplitMix64) -> ([f64; 8], [f64; 8], [f64; 8]) { + // A random sign and significand with biased exponent `e` (1..=2046). + let with_exp = |r: &mut SplitMix64, e: u64| f64::from_bits((r.next() & 0x800F_FFFF_FFFF_FFFF) | e << 52); + let sign = |r: &mut SplitMix64| if r.next() & 1 == 0 { 1.0 } else { -1.0 }; + let (lo, hi, chi) = (2f64.powi(-460), 2f64.powi(460), 2f64.powi(960)); + let lo_in = |r: &mut SplitMix64| { + if r.next() & 3 == 0 { + sign(r) * lo + } else { + with_exp(r, 1023 - 460) + } + }; + let hi_in = |r: &mut SplitMix64| { + if r.next() & 3 == 0 { + sign(r) * hi + } else { + with_exp(r, 1023 + 459) + } + }; + let mut a = [0f64; 8]; + let mut b = [0f64; 8]; + let mut c = [0f64; 8]; + // Lane 0: the smallest products of the vector path, a small c. + a[0] = lo_in(rng); + b[0] = lo_in(rng); + c[0] = { + let e = 1023 - 920 - 40 - rng.next() % 80; + with_exp(rng, e) + }; + // Lane 1: the largest products, nearly cancelled by c. + a[1] = hi_in(rng); + b[1] = hi_in(rng); + c[1] = -(a[1] * b[1]) * (1.0 + f64::EPSILON * (rng.next() % 64) as f64); + // Lane 2: |a| one binade below the bound, or |c| just above 2^960. + if rng.next() & 1 == 0 { + a[2] = with_exp(rng, 1023 - 461); + b[2] = lo_in(rng); + c[2] = { + let e = 1023 - 921 - 40 - rng.next() % 80; + with_exp(rng, e) + }; + } else { + a[2] = hi_in(rng); + b[2] = with_exp(rng, 1023 + 400); + c[2] = f64::from_bits(with_exp(rng, 1023 + 960).to_bits() | 1); + } + // Lane 3: inside, with |c| exactly 2^960. + a[3] = hi_in(rng); + b[3] = { + let e = 1023 + 400 + rng.next() % 60; + with_exp(rng, e) + }; + c[3] = sign(rng) * chi; + // Lane 4: operands of 2^-496 to 2^-545, so the product lies near or below + // the subnormal range, where Dekker's error terms underflow; a tiny c. + a[4] = { + let e = 1023 - 496 - rng.next() % 50; + with_exp(rng, e) + }; + b[4] = { + let e = 1023 - 496 - rng.next() % 50; + with_exp(rng, e) + }; + c[4] = f64::from_bits(rng.next() & 0x800F_FFFF_FFFF_FFFF); + // Lanes 5 and 7: inside, so each witness lane shares its v128 with an + // in-range lane and reaches the vector path if the guard is widened. + a[5] = lo_in(rng); + b[5] = hi_in(rng); + c[5] = { + let e = 1023 - 40 + rng.next() % 80; + with_exp(rng, e) + }; + a[7] = { + let e = 1023 - 200 + rng.next() % 400; + with_exp(rng, e) + }; + b[7] = { + let e = 1023 - 200 + rng.next() % 400; + with_exp(rng, e) + }; + c[7] = -(a[7] * b[7]); + // Lane 6, by turns: |a * b| in [2^1024, 2^1026) with c bringing the sum + // back into range; |a| just above 2^460; an infinity, NaN or zero. + match rng.next() % 3 { + 0 => { + a[6] = with_exp(rng, 1023 + 524); + b[6] = with_exp(rng, 1023 + 500); + c[6] = if (a[6] * b[6]).is_sign_negative() { + f64::MAX + } else { + -f64::MAX + }; + } + 1 => { + a[6] = f64::from_bits(with_exp(rng, 1023 + 460).to_bits() | 1); + b[6] = with_exp(rng, 1023 + 300); + c[6] = with_exp(rng, 1023 + 700); + } + _ => { + // A zero, infinity or NaN in a or b, against every kind of c. + let special = [f64::INFINITY, f64::NEG_INFINITY, f64::NAN, 0.0, -0.0]; + let s = special[(rng.next() % 5) as usize]; + let x = { + let e = 1023 + rng.next() % 64; + with_exp(rng, e) + }; + (a[6], b[6]) = if rng.next() & 1 == 0 { (s, x) } else { (x, s) }; + let cs = [0.0, -0.0, f64::INFINITY, f64::NEG_INFINITY, f64::NAN]; + let k = rng.next() % 7; + c[6] = if k < 5 { cs[k as usize] } else { with_exp(rng, 1023) }; + } + } + (a, b, c) +} + /// `I32x16::saturating_rounding_doubling_high_mul`, `rounding_divide_by_pot` and /// `multiply_by_quantized_multiplier` against an independent `i128` reference /// written from the definitions (no call into the scalar functions of diff --git a/scripts/masking-parity.sh b/scripts/masking-parity.sh index ccc9302c..05b7887e 100755 --- a/scripts/masking-parity.sh +++ b/scripts/masking-parity.sh @@ -7,6 +7,8 @@ # scripts/masking-parity.sh nightly `cargo +nightly --features nightly-simd` # scripts/masking-parity.sh wasm wasm32 +simd128 under node (simd_wasm arm) # scripts/masking-parity.sh wasm-scalar wasm32 WITHOUT simd128 under node = the scalar arm +# scripts/masking-parity.sh wasm-relaxed wasm32 +simd128 +relaxed-simd under node (simd_wasm arm; no source path +# selects on relaxed-simd, so this guards against one appearing) # scripts/masking-parity.sh neon-qemu aarch64 cross-build run under qemu-aarch64-static # scripts/masking-parity.sh avx-qemu AVX-without-AVX2 build (config-avx) under qemu -cpu SandyBridge # scripts/masking-parity.sh v2-qemu x86-64-v2 build (config-v2, scalar realization) under qemu -cpu Nehalem @@ -15,7 +17,7 @@ # realization on the native row) pass through CARGO_ARGS, see # scripts/codegen-witness.sh for why that form and not the env var. set -euo pipefail -ARM="${1:?native|nightly|wasm|wasm-scalar|neon-qemu|avx-qemu|v2-qemu}" +ARM="${1:?native|nightly|wasm|wasm-scalar|wasm-relaxed|neon-qemu|avx-qemu|v2-qemu}" ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" MANIFEST="$ROOT/crates/simd-masking-parity/Cargo.toml" TD="${CARGO_TARGET_DIR:-$ROOT/crates/simd-masking-parity/target}" @@ -50,8 +52,9 @@ case "$ARM" in cargo +nightly ${CARGO_ARGS:-} build --release --manifest-path "$MANIFEST" --bin simd-masking-parity --features nightly-simd "$TD/release/simd-masking-parity" ;; - wasm|wasm-scalar) + wasm|wasm-scalar|wasm-relaxed) FLAGS=""; [ "$ARM" = wasm ] && FLAGS="-C target-feature=+simd128" + [ "$ARM" = wasm-relaxed ] && FLAGS="-C target-feature=+simd128,+relaxed-simd" RUSTFLAGS="$FLAGS" cargo ${CARGO_ARGS:-} build --release --lib --manifest-path "$MANIFEST" --target wasm32-unknown-unknown echo "==> $ARM (RUSTFLAGS='$FLAGS') under node" node "$ROOT/crates/simd-masking-parity/run.mjs" "$TD/wasm32-unknown-unknown/release/simd_masking_parity.wasm" diff --git a/src/simd_wasm.rs b/src/simd_wasm.rs index 5b1274b7..46200966 100644 --- a/src/simd_wasm.rs +++ b/src/simd_wasm.rs @@ -45,10 +45,25 @@ //! IEEE TC 2008), including for subnormal results. //! - `F64x8`: Dekker's exact product, two-sum with `c`, a rounded-to-odd sum //! of the two error terms, then one rounding (the same paper's emulated -//! FMA). It needs no overflow or underflow in the split, so a v128 pair whose -//! lanes leave `2^-460 <= |a|, |b| <= 2^460`, `|c| <= 2^960` (or are NaN, -//! infinite or zero) is computed lane by lane with the integer FMA in -//! `simd_fma.rs` (no libm call, so no run-time CPU dispatch). +//! FMA). It needs no overflow or underflow in the split, so it runs only on +//! lanes with `2^-460 <= |a|, |b| <= 2^460` and `|c| <= 2^960`. A lane whose +//! `a` or `b` is a zero, an infinity or NaN has an exact product, so +//! `a * b + c` rounds once and is used instead, still in the vector. A v128 +//! pair with any other lane is computed lane by lane with the integer FMA +//! in `simd_fma.rs` (no libm call, so no run-time CPU dispatch). +//! +//! Inside the bounds the result is proven correctly rounded: +//! `.claude/knowledge/wasm-fma-proof-obligations.md`. The bounds are +//! sufficient, not tight. Dekker's product needs only `e_a + e_b >= -970` +//! (the guard gives `-920`). Given the bounds on `a` and `b`, a finite `c` +//! cannot make `c + a * b` overflow, so `|c| <= 2^960` serves only to +//! exclude infinite and NaN `c`, which the two-sum would turn into NaN. +//! Measured with the guard removed, against a hardware FMA: the emulation +//! first fails for products near `2^-1000` (operands near `2^-497`) and +//! for products above `f64::MAX`. That is evidence about where the +//! algorithm breaks, not a reason to widen the guard. +//! Measured on V8 (node 22): ~6-9 ns per element on the vector path, ~50-75 +//! ns on the lane path, ~1 ns for an unfused multiply-add. //! * **reduce_sum.** Uses a balanced tree reduction (pairwise v128 adds) //! rather than the scalar fallback's sequential left-fold, so the result //! can differ in the last ULP — the same reduction-order tolerance the @@ -562,6 +577,17 @@ pub mod wasm32_simd { if i64x2_all_true(ok) { return fma_f64_emulated(a, b, c); } + // A lane whose `a` or `b` is a zero, an infinity or NaN has an + // exact product (a signed zero, an infinity or NaN), so + // `a * b + c` rounds once: the IEEE result, and the same identity + // `simd_fma::fma_f64` starts with. Such pairs stay in the vector. + let ordinary = |x: v128| { + v128_and(f64x2_lt(f64x2_abs(x), f64x2_splat(f64::INFINITY)), f64x2_ne(x, f64x2_splat(0.0))) + }; + let exact_product = v128_not(v128_and(ordinary(a), ordinary(b))); + if i64x2_all_true(v128_or(ok, exact_product)) { + return v128_bitselect(fma_f64_emulated(a, b, c), f64x2_add(f64x2_mul(a, b), c), ok); + } let lane = |i: usize| match i { 0 => crate::simd_fma::fma_f64( f64x2_extract_lane::<0>(a),