simd(wasm): certify the F64x8::mul_add domain, keep exact-product lanes in the vector

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

.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 = "<pinned>"`; 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 |

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

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