diff --git a/tests/reference/upernet/helpers.py b/tests/reference/upernet/helpers.py
index bdd1bfb..d34667e 100644
--- a/tests/reference/upernet/helpers.py
+++ b/tests/reference/upernet/helpers.py
@@ -1,11 +1,14 @@
-"""Small NumPy and PyTorch helpers the ported decoder tests use, copied unchanged (apart from
-their imports) from research modules archived at the tag ``research-archive-2026-10-12``:
+"""Small NumPy and PyTorch helpers copied unchanged (apart from imports and the historical
+fixture's name) from research modules archived at ``research-archive-2026-10-12``,
+plus a new portable rounding probe:
 
 - ``pixel_ids_numpy``: ``upernet/earth_engine/feature_interface.py``;
 - ``fold``: ``upernet/earth_engine/normalization.py`` (BatchNorm folded into a conv weight;
   the decoder runs explicit BatchNorm instead, and this shows why);
-- ``grouped_probe_inputs``: ``upernet/earth_engine/validate_fusion.py`` (the Phase 16 fixture
-  that separates the fusion formulations);
+- ``grouped_probe_inputs_historical``: ``upernet/earth_engine/validate_fusion.py``
+  (the original Phase 16 fixture, unchanged);
+- ``grouped_probe_inputs``: a new portable fixture whose double sums are independent
+  of accumulation order;
 - ``torch_resize_matrix``: ``upernet/reference/m12_architecture_contract.py`` (PyTorch's own
   bilinear resize matrix, read off one-hot inputs).
 """
@@ -34,6 +37,27 @@ def fold(weight: np.ndarray, gamma, beta, mean, var, eps: float = BN_EPS) -> tup
 
 
 def grouped_probe_inputs():
+    """A one-window fixture (7 x 7 tokens, four groups of two channels) whose fusion conv
+    separates A's float32 partials from B and the monolithic convolution. Channel 0
+    is 2^40 + 1 (group 1) - 2^40 (group 2); channel 1 is 2^40 + 1 (group 3) + 1 -
+    2^40 (group 4). All products and every subtotal are integers of magnitude at most
+    2^41 + 2, exactly representable in double in any accumulation order. Monolithic
+    and B give (1, 2); A rounds away the units in each partial and gives (0, 0).
+    This is an offline probe, not the historical live fixture's 2^53 example."""
+    x = np.zeros((7, 7, 8), np.float32)
+    x[3, 3, [0, 1, 2]] = [2.0 ** 20, 1.0, -(2.0 ** 20)]
+    x[3, 3, [4, 5, 6, 7]] = [2.0 ** 20, 1.0, 1.0, -(2.0 ** 20)]
+    w = np.zeros((2, 8, 3, 3), np.float32)
+    w[0, [0, 1, 2], 1, 1] = [2.0 ** 20, 1.0, 2.0 ** 20]
+    w[1, [4, 5, 6, 7], 1, 1] = [2.0 ** 20, 1.0, 1.0, 2.0 ** 20]
+    params = {R.WEIGHT: w, f"{R.PREFIX}.1.weight": np.ones(2, np.float32),
+              f"{R.PREFIX}.1.bias": np.zeros(2, np.float32),
+              f"{R.PREFIX}.1.running_mean": np.zeros(2, np.float32),
+              f"{R.PREFIX}.1.running_var": np.ones(2, np.float32)}
+    return x, w, params
+
+
+def grouped_probe_inputs_historical():
     """A one-window fixture (7 x 7 tokens, four groups of two channels) whose fusion conv
     reveals each formulation's rounding.  Output channel 0 at the centre token is
     2^40 (group 1) + 1 (group 1) - 2^40 (group 2); channel 1 is 2^53 + 1 (group 3) + 1 -
diff --git a/tests/test_upernet_fusion_optimized.py b/tests/test_upernet_fusion_optimized.py
index 6b51c1b..509cc80 100644
--- a/tests/test_upernet_fusion_optimized.py
+++ b/tests/test_upernet_fusion_optimized.py
@@ -66,16 +66,53 @@ def test_grouped_double_differs_only_at_rounding_boundaries(case):
 
 
 def test_rounding_probe_fixture_separates_the_formulations():
-    """Phase 16 probe (run live by the archived validate_fusion): channel 0 = 2^40 + 1 - 2^40 across groups,
-    channel 1 = 2^53 + 1 | + 1 - 2^53: monolithic (1, 0), A (0, 0), B (1, 1)."""
+    """Exact double sums in any order: monolithic and B (1, 2), rounded A (0, 0).
+
+    Each product is an integer and the sum of absolute products is at most 2^41 + 2,
+    below double's exact-integer limit. Float32 partials still lose the unit terms.
+    """
     from reference.upernet.helpers import grouped_probe_inputs
     x, w, params = grouped_probe_inputs()
     resized = {R.RESIZED[n]: x[..., 2 * g:2 * g + 2].astype(np.float64) for g, n in enumerate(R.PYRAMID)}
-    assert list(MF.conv3x3(x, w)[3, 3]) == [1.0, 0.0]
+    assert list(MF.conv3x3(x, w)[3, 3]) == [1.0, 2.0]
+    assert list(MU.grouped_rounded(resized, params)["grouped_sum"][3, 3]) == [0.0, 0.0]
+    assert list(MU.grouped_double(resized, params)["double_sum"][3, 3]) == [1.0, 2.0]
+
+
+def test_historical_rounding_probe_grouped_formulations():
+    """Original 2^53 probe: real grouped A (0, 0), B (1, 1).
+
+    Each group dot has at most two nonzero terms, so its result is independent of
+    the BLAS accumulation order. The four partials are added in MU's fixed order.
+    """
+    from reference.upernet.helpers import grouped_probe_inputs_historical
+    x, _, params = grouped_probe_inputs_historical()
+    resized = {R.RESIZED[n]: x[..., 2 * g:2 * g + 2].astype(np.float64) for g, n in enumerate(R.PYRAMID)}
     assert list(MU.grouped_rounded(resized, params)["grouped_sum"][3, 3]) == [0.0, 0.0]
     assert list(MU.grouped_double(resized, params)["double_sum"][3, 3]) == [1.0, 1.0]
 
 
+def test_historical_rounding_probe_with_explicit_sequential_accumulation():
+    """Historical Phase 16 monolithic dot in sequential order gives (1, 0).
+
+    This result depends on accumulation order: sequential 2^53 + 1 + 1 - 2^53
+    loses both units, while regrouping can retain them. Use fixed scalar additions,
+    never a BLAS/MKL reduction whose SIMD path may reassociate the double sum.
+    A and B are checked through the real grouped code in the test above.
+    """
+    from reference.upernet.helpers import grouped_probe_inputs_historical
+    x, w, _ = grouped_probe_inputs_historical()
+    # Only the centre tap is nonzero; preserve its channel order explicitly.
+    products = x[3, 3].astype(np.float64) * w[:, :, 1, 1].astype(np.float64)
+    monolithic = []
+    for channel in products:
+        total = 0.0
+        for term in channel:
+            total = total + float(term)
+        monolithic.append(float(np.float32(total)))
+    assert monolithic == [1.0, 0.0]
+
+
 def test_order_study(case):
     """Phase 17: B is insensitive to the group order except at rounding boundaries; A's result
     depends on the order in a fraction of the outputs; all within 1e-6 of float64."""