fix(ops): floor_int cancellation-free, compiled-precise saturating staircase

floor_int built its staircase via piecewise_linear, which sums slope-change
ReLUs Σ ±s·ReLU(x−b_i) in one linear projection. For a tall staircase at large
x this is catastrophic: each term is ~s·x and the fp32 partial sums overflow
2^24 (16.7M) so the +/- terms cancel into garbage. floor_int(40000.4, [0,65535],
s=1000) returned 0.0 (true 40000); even [0,255] staircases lost ~1-2 units near
the top of the range, unstable to 1-ULP input changes.

Rebuild as a sum of bounded saturating steps:
  floor(x) = min + Σ_k step_k,  step_k = clamp(s·(x−k)+1, 0, W) = relu(t)−relu(t−W)
  floor    = min + n − Σ_k relu(1 − step_k).
Each step_k is a 2-term per-output difference (not a sum over all breakpoints),
so the final Σ accumulates only bounded [0,1] terms — partial sums never exceed
n. The step is clamped to [0, W] (not left unbounded), keeping the stage-1
intermediate small so it does not amplify upstream compiled noise (an unbounded
relu(t) intermediate broke the linear_bin_index compile-fidelity probe). W is
sized max(2, 8·ulp(s·n)) — 2 for the common small-range case, larger only as s·n
approaches 2^24, where relu(t)−relu(t−W) would otherwise collapse to 0. Exact at
any magnitude/sharpness AND compiled-precise. Cost: one extra MLP sublayer (two
chained ReLUs). ceil_int = −floor_int(−x) inherits the fix.

thermometer_floor_div shares the slope-change-sum structure but is correct for
its documented integer-inputs-only contract (verified exact on large integers) —
not reworked.

Tests: test_floor_int_wide_range_large_magnitude (s=10/1000 over [0,65535]),
test_floor_int_byte_range_exact_high_sharpness, test_ceil_int_wide_range; full
arithmetic_ops + resampling (compiled probe) + value_type suites green. Noise
re-measured (the [-5,10] distributions never probed the bug — see
numerical_noise_findings.md).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Author: physicsrob

docs/numerical_noise.md

@@ -1,7 +1,7 @@
 # Numerical noise reference
 
 *Generated by `scripts/measure_op_noise.py`; do not edit by hand.*
-*Measured at commit `fe133f9` on `2026-06-02`.*
+*Measured at commit `2e04d93` on `2026-06-02`.*
 
 Every piecewise-linear op in `torchwright/ops/` is measured on one or more
 named input distributions characterising its expected input ranges.
@@ -296,7 +296,7 @@ Uniform continuous samples on [-5, 10] — most land in flat zones; measures bas
 | --- | --- |
 | Samples | 4096 |
 | Max abs error | 0.9982 |
-| Mean abs error | 0.04655 |
+| Mean abs error | 0.04654 |
 | p99 abs error | 0.8959 |
 | Max rel error | 0.953 (over 3828 samples with |ref| ≥ 1e-6) |
 | Mean rel error | 0.01849 |

docs/numerical_noise_findings.md

@@ -127,6 +127,71 @@ leakage of noise onto the validity bit gets multiplied by 1000 inside the
 attention logit. This op is not yet wired into the measurement pipeline
 and is a reasonable follow-up.
 
+### `floor_int` staircase is cancellation-free (saturating steps), not a slope-change ReLU sum
+
+`floor_int` (and `ceil_int = -floor_int(-x)`) previously realised its
+staircase through `piecewise_linear`, which sums slope-change ReLUs
+`Σ ±s·ReLU(x − b_i)` in **one** linear projection. For a tall staircase
+evaluated at large `x` this is catastrophic: each term has magnitude
+`~s·x`, and the partial sums during the fp32 dot-product accumulation
+overflow float32's `2^24 ≈ 16.7M` exact-integer limit, so the alternating
+`+`/`−` terms cancel into garbage. `floor_int(40000.4, [0,65535], s=1000)`
+returned `0.0` (true `40000`); even `[0,255]` staircases lost ~1–2 units
+near the top of the range. **The measured noise distributions never caught
+this** — `floor_uniform_neg5_10` / `floor_near_boundary_10` sample only
+`[-5,10]`, far inside the safe regime — which is the cautionary part: a
+clean per-op number does not certify an op outside its measured input
+range.
+
+The op now sums **bounded saturating** steps, `floor(x) = min + Σ_k step_k`
+with `step_k = clamp(s(x−k)+1, 0, W) = relu(t) − relu(t−W)`, counted as
+`floor = min + n − Σ_k relu(1 − step_k)`. Each `step_k` is a 2-term
+per-output difference (not a sum over all breakpoints), so the final `Σ`
+accumulates only bounded `[0,1]` terms — partial sums never exceed `n`. Two
+subtleties matter:
+
+- **Bound the step, don't leave it unbounded.** An earlier version output the
+  raw `relu(t)` (∈ `[0, ~s·x]`) and relied on `1 − relu(1 − relu(t))` to
+  saturate. That is reference-exact but the *intermediate* `relu(t)` carries
+  upstream compiled noise amplified by `s` (≈×10), which broke the
+  `linear_bin_index` compile-fidelity probe (`test_resampling_primitives`) on
+  an internal node even though the output was correct. Clamping the step to
+  `[0, W]` keeps the intermediate small (compiled-precise) with `≥ W−1` slack
+  below the stage-2 threshold of 1.
+- **Size `W` to the magnitude.** `relu(t) − relu(t−W)` collapses to 0 once `t`
+  and `t−W` round to the same float32 (`ulp(t) ≥ W`), so `W = max(2,
+  8·ulp(s·n))` — 2 for the common small-range case, larger only as `s·n`
+  approaches `2^24`.
+
+Cost: one extra MLP sublayer (two chained ReLUs) and a width-`n` intermediate.
+Exact at any magnitude/sharpness AND compiled-precise; pinned by
+`test_floor_int_wide_range_large_magnitude`,
+`test_floor_int_byte_range_exact_high_sharpness`, and the existing
+`linear_bin_index` compile probe. Note `thermometer_floor_div` shares the
+slope-change-sum structure but is correct for its documented **integer-only**
+contract (integer inputs land in flat zones); it is not reworked.
+
+### Known gap: the `map_select` table-lookup staircases have the same 2^24 cancellation
+
+An audit for the `floor_int` cancellation class found one untriggered sibling.
+`_table_lookup_index_staircase` / `_table_lookup_row_vector` /
+`_table_lookup_column_mask` (`torchwright/ops/map_select.py`) build a
+`piecewise_linear` with `2·(n−1)+1` breakpoints and `input_scale=sharpness`
+(default 100), feeding `table_lookup_2d`/`table_lookup_3d`. Same single
+alternating-sign ReLU sum, so the same overflow:
+`_table_lookup_index_staircase(n=256, sharpness=100)` at `x=100000` returns
+`0.0` (should clamp to 255 — the index is never clamped before the staircase),
+and `n=20000, sharpness=1000` collapses an in-range `x=19999` to `0.0`
+(`s·max_breakpoint ≈ 2e7 > 2^24`). **Untriggered today**: current table
+consumers use small tables with in-range, pre-rounded integer indices, and the
+DOOM texture/palette lookup track is still a Phase-F stub. Fix before that
+track lands large tables — clamp the index into `[0, n−1]` *before* the
+staircase, and/or apply the saturating-step reformulation. The audit also
+flagged `compare` as a milder 2-term variant that fails silently for `|x| ≳
+2^24/sharpness` (default-sharpness boundary ~1.68e6; live callers stay far
+under it), and confirmed `square`/`exp`/`multiply_2d`/`low_rank_2d` safe
+(monotone slopes or `input_scale=1`).
+
 ### Rel-error reporting: ramp-zone artefacts
 
 For `compare`, `floor_int`, `linear_bin_index`, and the boolean ops

docs/op_noise_data.json

@@ -1,6 +1,6 @@
 {
   "schema_version": 1,
-  "commit": "fe133f9",
+  "commit": "2e04d93",
   "measured_at": "2026-06-02",
   "ops": [
     {
@@ -254,12 +254,12 @@
           "name": "floor_near_boundary_10",
           "description": "Samples within ±1/step_sharpness of integer k=5 — stresses the ramp zone where `floor_int` interpolates.",
           "n_samples": 4096,
-          "max_abs_error": 0.999268,
+          "max_abs_error": 0.999279,
           "mean_abs_error": 0.128753,
-          "p99_abs_error": 0.970386,
-          "max_rel_error": 0.249817,
+          "p99_abs_error": 0.970371,
+          "max_rel_error": 0.24982,
           "mean_rel_error": 0.0321882,
-          "p99_rel_error": 0.242596,
+          "p99_rel_error": 0.242593,
           "rel_valid_samples": 4096,
           "worst_input": {
             "x": 4.99993
@@ -269,12 +269,12 @@
           "name": "floor_uniform_neg5_10",
           "description": "Uniform continuous samples on [-5, 10] — most land in flat zones; measures baseline `floor_int` accuracy.",
           "n_samples": 4096,
-          "max_abs_error": 0.998169,
-          "mean_abs_error": 0.0465452,
+          "max_abs_error": 0.998161,
+          "mean_abs_error": 0.0465444,
           "p99_abs_error": 0.895853,
           "max_rel_error": 0.953022,
-          "mean_rel_error": 0.0184898,
-          "p99_rel_error": 0.46109,
+          "mean_rel_error": 0.0184896,
+          "p99_rel_error": 0.461094,
           "rel_valid_samples": 3828,
           "worst_input": {
             "x": 5.99982

tests/ops/test_arithmetic_ops.py

@@ -895,6 +895,52 @@ def test_ceil_int():
         ), f"ceil({v}) = {expected}, got {result.item()}"
 
 
+def test_floor_int_wide_range_large_magnitude():
+    """Regression: the old slope-change ReLU staircase summed ~n terms of
+    magnitude ~s*x in one projection, overflowing fp32's 2^24 exact-integer
+    limit and collapsing to ~0 above x ~ 2^24/s (e.g. floor_int(40000.4,
+    [0,65535], s=1000) returned 0.0). The saturating-step staircase keeps every
+    partial sum in [0, n] so it stays exact at any magnitude/sharpness."""
+    import math
+
+    # Mid-bin values (0.4 above an integer, clear of any ramp for s >= 10) so the
+    # only thing under test is the staircase summation, not near-boundary ramps.
+    x = create_input("x", 1)
+    for s in (10.0, 1000.0):
+        f = floor_int(x, min_value=0, max_value=65535, sharpness=s)
+        for v in (10.4, 1000.4, 16800.4, 40000.4, 62195.4, 65534.4):
+            r = f.compute(n_pos=1, input_values={"x": torch.tensor([[v]])}).item()
+            assert abs(r - math.floor(v)) < 0.01, (
+                f"floor_int({v}, [0,65535], s={s}) should be {math.floor(v)}, "
+                f"got {r}"
+            )
+
+
+def test_floor_int_byte_range_exact_high_sharpness():
+    """Exactness guard over the [0,255] byte range at high sharpness — the DOOM
+    digit-quad high byte floor(q/256). The saturating-step form must stay exact
+    across the whole range (the old form's partial sums grew with the step count
+    and position and could lose ~1-2 units near the top)."""
+    import math
+
+    x = create_input("x", 1)
+    f = floor_int(x, min_value=0, max_value=255, sharpness=1000.0)
+    for v in (0.4, 1.4, 100.6, 127.95, 200.3, 242.953, 254.6, 255.0):
+        r = f.compute(n_pos=1, input_values={"x": torch.tensor([[v]])}).item()
+        assert abs(r - math.floor(v)) < 0.01, f"floor_int({v}, [0,255]) got {r}"
+
+
+def test_ceil_int_wide_range_large_magnitude():
+    """ceil_int = -floor_int(-x); inherits the saturating-step fix."""
+    import math
+
+    x = create_input("x", 1)
+    c = ceil_int(x, min_value=0, max_value=65535)
+    for v in (10.4, 1000.4, 40000.4, 62195.98):
+        r = c.compute(n_pos=1, input_values={"x": torch.tensor([[v]])}).item()
+        assert abs(r - math.ceil(v)) < 0.01, f"ceil_int({v}) should be {math.ceil(v)}, got {r}"
+
+
 def test_signed_multiply():
     """Test all sign combinations."""
     a = create_input("a", 1)