Two papers submitted the same day, September 22, 2026, decompose the same failure from two directions: greedy decoding — argmax at temperature zero, seed fixed, the configuration everyone labels "deterministic" — does not reproduce. The first, Greedy Decoding Is Not Precision-Invariant (Du, Khan, Zhou, Liu, Chakrabarti, Suya, Li; arXiv:2609.26621, accepted at TMLR), shows the intra-precision axis: the same checkpoint, prompt, and algorithm produce different token sequences in BF16 versus FP16 on the same GPU1. The second, Accelerating the Mitigation of LLM Inference Nondeterminism Across GPU Architectures (Cooper, Jeong, Jeon, Young, Kim; arXiv:2609.25624), shows the cross-device axis: the same model and stack on an A100, L40S, and H100 diverge because frameworks select different GEMM kernels — hence different floating-point reduction orders — per architecture2. The shared root is the oldest fact in numerical computing: addition is not associative, so the order in which a dot product is summed changes its result at the last bits, and greedy argmax has zero tolerance for a last-bit rank flip.
The first paper's numbers: across six models (1.1B–7B, four families — Llama, Qwen, Mistral, OLMoE-MoE; divergence additionally characterized at 12B) and three benchmarks (GSM8K, HumanEval, MBPP), 49–100% of greedy generations diverge between BF16 and FP16. On TinyLlama-1.1B only 30–41% come out identical; one flipped token has a mean downstream length difference of 34 tokens; on Qwen2.5-3B GSM8K, 19% of prompts flip final-answer correctness while aggregate accuracy moves 1pp — the drift hides inside benchmark scores1. Standard reproducibility tooling does not help: torch.use_deterministic_algorithms and a fixed CUBLAS workspace make repeated runs bit-identical within each format (BF16 self-consistency is 100%) but cannot reconcile two formats that genuinely compute different logits.
1. The mechanism: a margin event, not an error-magnitude event
The paper's central finding is a negative result about the transformer body. Hidden-state error between the two precisions accumulates — roughly linearly across all 22 TinyLlama layers, from L2 of 0.002 at layer 0 to about 0.92 at layer 21 — but it is present in every prompt, flipping or not. At the divergence step, Diverged and Agreed prompts show indistinguishable body-layer L2 (1.022 vs 1.135) and logit L2 (4.928 vs 5.559); what separates them by two-plus orders of magnitude is the top-two logit margin at the lm_head: 0.039 versus 5.7331. Twenty-two layers of accumulated body error do not predict which steps flip. The margin does.
The paper formalizes the flip with an exact algebraic condition (Proposition 1): given logit vectors z and z-hat from two configurations with top token v(1), the argmax flips if and only if some competitor v satisfies
∆z(v) − ∆z(v(1)) > z(v(1)) − z(v),
that is, if the directional perturbation between the two candidates exceeds the original gap between them1. For the runner-up v(2), this is precisely the condition for v(2) to overtake v(1). The proposition is elementary — it is a rearrangement — but it cleanly separates the two quantities that matter: the margin (a property of the model's decision) and the directional perturbation (a property of the arithmetic). Empirically, the per-logit perturbation scale on TinyLlama is σz ≈ 0.026 (logit L2 of 4.63 across a 32,000-token vocabulary), giving directional differences a scale of about √2 · σz ≈ 0.037. Flips concentrate where the margin sits below that scale; a decision rule built on the directional statistic classifies flip versus non-flip steps correctly on 88–89% of steps across two model families1.
The second quantitative piece explains where the perturbation comes from: the lm_head itself. Under a residual-connection model where per-layer rounding stays confined to a small non-residual branch (residual ratio α ≈ 0.05), the body-to-head perturbation ratio is roughly L·α/√d — for TinyLlama, 22 · 0.05/√2048 ≈ 0.024, i.e. the head matmul contributes about 97.6% of the per-logit perturbation and the 22-layer body about 2.4%. This is a heuristic estimate with stated assumptions, not a theorem, but it matches the measured decomposition and, as Section 3 shows, correctly predicts what happens when you push FP32 compute beyond the head1.
2. Runnable: non-associativity, the flipping argmax, and the margin criterion
Three cells reproduce the paper's logic in miniature. All are plain Python, seeded, and their outputs below are from real runs.
Cell 1 — the root cause. Floating-point addition rounds after every operation, so summation order is part of the arithmetic: the same sum in different association orders gives different answers, and a rounding error's absolute size grows with magnitude, so there is no global epsilon.
import struct, math
def to_f32(x):
return struct.unpack('f', struct.pack('f', x))[0]
a, b, c = 1.0e8, -1.0e8, math.pi
print("three-term sum, a=1e8, b=-1e8, c=pi, computed in IEEE double:")
print(" (a + b) + c =", repr((a + b) + c))
print(" a + (b + c) =", repr(a + (b + c)))
print(" equal? ", (a + b) + c == a + (b + c))
terms = [16777216.0, 1.0, 1.0, -16777216.0] # exact sum = 2.0
sA = to_f32(to_f32(to_f32(terms[0] + terms[1]) + terms[2]) + terms[3])
sB = to_f32(terms[0] + to_f32(terms[1] + to_f32(terms[2] + terms[3])))
print("\nfour-term sum in float32, two reduction orders (exact value 2.0):")
print(" values :", terms)
print(" serial order :", sA)
print(" pairwise tree:", sB)
print(" identical? ", sA == sB)
one = to_f32(1.0)
ulp1 = to_f32(one + 1.1920929e-7) - one
big = to_f32(1.0e8)
ulpb = to_f32(big + 8.0) - big
print("\nULP of float32 at 1.0 :", ulp1)
print("ULP of float32 at 1e8 :", ulpb)Output from a real run:
three-term sum, a=1e8, b=-1e8, c=pi, computed in IEEE double:
(a + b) + c = 3.141592653589793
a + (b + c) = 3.141592651605606
equal? False
four-term sum in float32, two reduction orders (exact value 2.0):
values : [16777216.0, 1.0, 1.0, -16777216.0]
serial order : 0.0
pairwise tree: 2.0
identical? False
ULP of float32 at 1.0 : 1.1920928955078125e-07
ULP of float32 at 1e8 : 8.0The four-term float32 case is the paper's whole story in miniature: 2^24 + 1 is not representable in float32, so the serial order swallows both +1 terms into the large accumulator and returns 0.0, while the pairwise tree groups the small terms first and returns the exact 2.0. Different reduction order, different bits — and a GEMM kernel's split-K factor, tile size, and accumulator count are exactly a choice of association order.
Cell 2 — the flip and the cascade. The paper's flips are not representable-tie artifacts (that is Intervention A's null result, Section 3); they come from low-precision arithmetic: a d-term dot product accumulated in BF16 versus FP16 lands differently, and if the top-two margin is smaller than that difference, the argmax reverses. This cell emulates the lm_head dot product in both formats — every multiply and every accumulate rounded to the format — over 300 random near-tie trials, then plays one flipping trial out as a greedy decode.
import struct, math, random
def to_f16(x):
return struct.unpack('e', struct.pack('e', x))[0]
def to_bf16(x):
w = struct.unpack('I', struct.pack('f', x))[0]
w += 0x7FFF + ((w >> 16) & 1) # round-to-nearest-even on the low 16 bits
w &= 0xFFFF0000
return struct.unpack('f', struct.pack('I', w))[0]
def dot_q(wrow, h, fmt):
r = to_bf16 if fmt == 'bf16' else to_f16
acc = r(0.0)
for wk, hk in zip(wrow, h):
acc = r(acc + r(r(wk) * r(hk)))
return acc
random.seed(20260925)
d = 64
ntrials = 300
nflip = 0
flip_example = None
for trial in range(ntrials):
h = [random.gauss(0, 1) for _ in range(d)]
w1 = [random.gauss(0, 1) for _ in range(d)]
w2 = [random.gauss(0, 1) for _ in range(d)]
t1 = sum(w*hk for w, hk in zip(w1, h)) # true (float64) logits
t2 = sum(w*hk for w, hk in zip(w2, h))
# shift w2[0] so the true top-two margin lands within +/-0.005
target = random.uniform(1e-4, 5e-3) * random.choice([1, -1])
w2[0] += (t1 - target - t2) / h[0]
t2 = sum(w*hk for w, hk in zip(w2, h))
b1, b2 = dot_q(w1, h, 'bf16'), dot_q(w2, h, 'bf16')
f1, f2 = dot_q(w1, h, 'fp16'), dot_q(w2, h, 'fp16')
flipped = (b1 > b2) != (f1 > f2)
if flipped:
nflip += 1
if flip_example is None:
flip_example = (t1 - t2, b1 - b2, f1 - f2, h, w1, w2)
tm, bm, fm, h, w1, w2 = flip_example
print(f"{ntrials} toy lm_head dot products (d={d}), true top-two margin within +/-0.005:")
print(f" BF16-vs-FP16 argmax flips: {nflip} of {ntrials} ({100*nflip/ntrials:.1f}%)")
print("\nfirst flipping trial (used as step 3 below):")
print(f" true margin {tm:+.5f} BF16 margin {bm:+.5f} FP16 margin {fm:+.5f}")
winner_bf = "D" if dot_q(w1, h, 'bf16') > dot_q(w2, h, 'bf16') else "E"
winner_fp = "D" if dot_q(w1, h, 'fp16') > dot_q(w2, h, 'fp16') else "E"
# toy continuation: the model wants different tokens after D vs after E
cont = {"D": ["A", "B"], "E": ["E", "C"]}
seq_bf = ["A", "C", winner_bf] + cont[winner_bf]
seq_fp = ["A", "C", winner_fp] + cont[winner_fp]
print("\ngreedy decode, same model, same prompt, identical hardware:")
print(" BF16 sequence:", seq_bf)
print(" FP16 sequence:", seq_fp)
first_div = next((t + 1 for t, (x, y) in enumerate(zip(seq_bf, seq_fp)) if x != y), None)
print(" first divergence at step:", first_div)
print(" single flipped token at step 3 cascades; tails never re-converge:", seq_bf != seq_fp)Output from a real run:
300 toy lm_head dot products (d=64), true top-two margin within +/-0.005:
BF16-vs-FP16 argmax flips: 155 of 300 (51.7%)
first flipping trial (used as step 3 below):
true margin +0.00182 BF16 margin -0.03711 FP16 margin +0.00415
greedy decode, same model, same prompt, identical hardware:
BF16 sequence: ['A', 'C', 'E', 'E', 'C']
FP16 sequence: ['A', 'C', 'D', 'A', 'B']
first divergence at step: 3
single flipped token at step 3 cascades; tails never re-converge: TrueTwo things to read out. First, the flip rate: with the true margin held under 0.005 — the regime the paper's Figure 1 shows for diverged prompts (median margin at the measurement step about 10−5, with 30 of 59 divergent prompts at an exact BF16 tie) — about half of head computations flip between formats. Outside that margin band, flips are essentially impossible. Second, the cascade: after step 3 the two "arms" condition on different tokens. Every later step can have an enormous, perfectly safe margin — the toy's continuation steps are decided by gaps of 8–9 logits — and the sequences still never re-converge, because greedy decoding has no mechanism to merge trajectories. This is the paper's 34-token mean length divergence, in five tokens.
Cell 3 — the paper's margin criterion, made concrete. Proposition 1 compares the margin against a directional perturbation. Treat the perturbation as a zero-mean random variable with scale √2·σz and the margin as a draw from a heavy-left-tailed distribution, and flip probability becomes a function of the margin alone.
import math, random, bisect
sigma = 0.026 # per-logit RMS perturbation (paper, TinyLlama)
sigma_dir = math.sqrt(2) * sigma # directional top-two scale
random.seed(3)
n = 200000
abs_deltas = sorted(abs(random.gauss(0, sigma_dir)) for _ in range(n))
def flip_prob(m): # P(|directional perturbation| > margin)
return 1.0 - bisect.bisect_right(abs_deltas, m) / n
random.seed(4)
margins = sorted(random.lognormvariate(-1.0, 1.5) for _ in range(20000))
print(f"directional perturbation scale sqrt(2)*sigma_z = {sigma_dir:.4f}")
print("top-two margin band P(flip at this step)")
bands = [(0.0,0.001),(0.001,0.01),(0.01,0.037),(0.037,0.1),(0.1,0.37),(0.37,1.0),(1.0,3.7)]
for lo, hi in bands:
print(f" {lo:>6g} - {hi:<6g} {flip_prob((lo+hi)/2):.4f}")
tau = 1e-3
print(f"fraction of sampled margins below tau=1e-3 : "
f"{sum(1 for m in margins if m < tau)/len(margins):.4f}")
print(f"fraction below directional scale {sigma_dir:.4f} : "
f"{sum(1 for m in margins if m < sigma_dir)/len(margins):.4f}")
p = sum(flip_prob(m) for m in margins) / len(margins)
print(f"unconditional per-step flip probability : {p:.4f}")
L = 170
print(f"with ~{L} decoding steps, P(no flip anywhere) = {(1.0 - p) ** L:.4f}")Output from a real run:
directional perturbation scale sqrt(2)*sigma_z = 0.0368
top-two margin band P(flip at this step)
0 - 0.001 0.9893
0.001 - 0.01 0.8814
0.01 - 0.037 0.5256
0.037 - 0.1 0.0625
0.1 - 0.37 0.0000
0.37 - 1 0.0000
1 - 3.7 0.0000
fraction of sampled margins below tau=1e-3 : 0.0001
fraction below directional scale 0.0368 : 0.0635
unconditional per-step flip probability : 0.0493
with ~170 decoding steps, P(no flip anywhere) = 0.0002The paper's measurement says the same with real distributions: margins are bimodal across five-plus orders of magnitude — diverged prompts at median 0 (30/59 exact BF16 ties) versus agreed at about 5.9. The per-step flip probability is tiny at typical steps and enormous in the left tail, which is why the paper's mitigation can gate on an observable (the current margin) and trigger on only 0.7–1.4% of steps while catching most flippable ones1. The second paper measures the same structure across GPUs and derives the sobering corollary: below 0.1 nats the margin CDF is close to linear, so the exposed fraction of steps is proportional to logit noise — "every order of magnitude of extra precision buys exactly one order of magnitude fewer flips and no threshold below which flips stop"2.
3. Five predictions, including one that sounds wrong
The mechanism is worth little without testable consequences, and the paper frames its intervention ablation as five predictions made before the experiments:
- Integer logit quantization (bins of width τ before argmax): zero effect — divergence is value-level, not a tie artifact.
- Temperature sharpening (logits times α above 1): zero effect at any α — margin and perturbation scale identically, their ratio invariant.
- Top-K FP32 recomputation works for all K at or above 2 — a σz-scale perturbation cannot bridge the rank-2-to-rank-3 gap (median about 0.5).
- Full-vocabulary FP32 lm_head recomputation matches top-K, at 4000× lower per-trigger cost when K = 8.
- Extending FP32 scope beyond the lm_head makes agreement worse — re-arithmetic at safe steps changes the logit landscape and introduces new divergence downstream (RMSNorm+lm_head drops TinyLlama EAR 63% → 44%; ungated τ=0 recomputation scores 54–60% at about 2.5× latency, below the gated 55–63% at under 4%)1.
All five were confirmed. Prediction 5 is the anti-hype lesson in miniature: more FP32 compute is not monotone improvement. The repair works because it is narrow and gated — it changes only the steps the margin flag identifies as flippable. Blanket precision upgrades perturb exactly the safe steps that were agreeing, and agreement drops1.
The gated intervention (Intervention C): compute native-precision logits; if the top-two margin is below τ = 10−3, recompute the full lm_head in FP32 (weights and hidden state upcast on the fly, no persistent FP32 copies). On TinyLlama this lifts exact agreement from 41% to 63% on GSM8K (36→61 HumanEval, 30→55 MBPP) at under 4% latency (+1.4% latency, +11% peak memory; triggers on 0.7–1.4% of steps). Gated C fires on exactly 30 of 100 GSM8K prompts and converts 22 — the stratification predicts the headline lift without tuning, and correctly declines to fire on the 29 prompts whose flip is driven by accumulated upstream error a single-step recompute cannot reconcile1.
4. The applicability map, and where the fix dies
Cross-model at batch size 1 on A10G, the lift lands in three tiers: +22pp TinyLlama-1.1B, +36pp Llama-3.2-3B (Tier 1); +3pp Qwen2.5-3B, +8pp Mistral-7B, +10pp OLMoE-1B-7B (Tier 2); 0pp on DS-R1-Distill-Qwen-7B, whose 100/100 prompts diverge with median first-divergence step 6 — the tested BF16-saturated Qwen variants produce FP16-body NaNs and first-token divergence a head-scoped repair cannot reach (Tier 3). The paper hypothesizes — explicitly as a hypothesis — that effectiveness tracks training-time precision stability rather than scale1. Hardware moves the numbers: on L4 and A100 the tested lifts are +12 to +21pp, with Qwen2.5-3B at +3pp on A10G but +20pp on L4 — the A10G is the outlier, and the headline range is one point in a hardware-dependent range.
Two boundaries are hard. First, batch: the lift degrades from +17pp at batch 1 to 0pp at batch 8, because batching changes reduction orders and effectively creates a third trajectory rather than a larger gap between two — measured hidden-state gap for batch-1 vs batch-8 (0.0127) is statistically indistinguishable from the working BF16-vs-FP16 configuration, so the gap itself is not what kills it. The method is a low-batch (up to 4) single-stream tool in the paper's own scope statement. Second, end-to-end FP8: quantizing only the lm_head to FP8-E4M3 (body shared) leaves a rescaled-gate repair effective — exact agreement on TinyLlama, +56pp on Qwen2.5-3B at τ = 0.25 — but quantizing all 155 linear modules makes divergence body-dominated and the same repair recovers only +1pp. Widening the gate does not help. There is also a hard floor no head-side fix can cross: even both arms computing in FP32 reach only 0.88 agreement when one stores BF16 weights and the other FP16, because the stored mantissas differ — about 12pp of irreducible weight-truncation gap1. The paper's own summary is the right one: a partial mitigation, not a determinism guarantee.
5. The second axis: kernel choice across architectures
The first paper held hardware fixed and varied precision. The second paper varies the hardware and finds a deeper channel2. Frameworks select GEMM kernels per architecture — cuBLAS heuristics keyed to SM count, tiles, split-K factors — and each kernel implies a different parallel reduction order, accumulated with tensor-core arithmetic whose internal behavior (truncated significand alignment, carry-out, intermediate width) is officially undocumented and demonstrably different across Volta, Ampere, Ada, and Hopper. A trap sits inside the obvious fix: "FP32" GEMMs on modern NVIDIA GPUs typically execute on tensor cores in TF32, which rounds inputs to a 10-bit significand, so a nominally FP32 pipeline silently inherits the unspecified arithmetic unless TF32 is disabled — and framework-level switches do not propagate into Triton-generated kernels2.
The state-of-the-art mitigation (LayerCast, Yuan et al., NeurIPS 2025 oral — the GEMM paper's baseline) stores weights in BF16 and upcasts each weight matrix to a transient FP32 copy for computation3. It is effective but doubly flawed: the GEMM reads the FP32 copy, so weight-memory traffic — the quantity that sets decode latency — stays at FP32 width; and the vendor GEMM still picks different reduction orders per device, so cross-GPU reproducibility stays statistical: on the paper's probe shapes, the cast-then-cuBLAS path was bitwise-identical on no shape across A100, L40S, and H100 — A100 and L40S disagreed on all seven, with 88–99% of output elements differing2.
The fix pins every device-dependent degree of freedom with four design rules: (R1) every dot product runs as IEEE-754 FMA on CUDA cores — bit-specified for given operands on every architecture — never tensor cores; (R2) no autotuning: kernel configurations are compile-time constants selected by problem shape alone (rerunning the same search on three GPUs picked the same configuration for only 3 of 9 representative shapes, and rerunning on one device changed its own answer on 1 of 9); (R3) deterministic split-K: partial sums combined in fixed ascending segment order, no atomics; (R4) batch invariance by construction — the K-segmentation depends only on K, never the batch dimension, so a request's outputs are bitwise independent of co-scheduled load for batch sizes 1 through 642. Weights load at 16-bit width and upcast to FP32 in registers, halving weight traffic.
The result is reproducibility by construction rather than by margin: the linear layers' reduction order is a pure function of problem shape, so cross-architecture agreement reduces to one premise — that each GPU implements IEEE-754 FP32 correctly — and the probe confirms bitwise-identical linear outputs across A100, L40S, and H100 on all seven shapes. End-to-end in vLLM, unmitigated BF16 diverges on 30.81–100% of problems (median about 85%) across GPU pairs; the FP32-compute baseline shrinks that under 1% but leaves residues up to 0.51% on five of twelve cells; the fixed-order kernels are at or below the baseline everywhere and zero in eleven of twelve — the single nonzero cell (0.08%, Qwen3-4B GSM8K) is one token at an exact FP32 tie, traceable to an unpinned layer-0 attention kernel on H100, not the linear layers. And the performance claim, verified at the source: faster end-to-end than the cast-then-GEMM state of the art in all 36 model-benchmark-GPU configurations, by 1.17–1.43× on A100 and H100 and 1.6–3.1× on the bandwidth-lean L40S, plus 1–1.4 GiB weight savings converted into KV-cache capacity on untied-embedding models2. The scope matters: that comparison is against FP32-compute mitigation at batch 32, models 3B–14B, single GPU — not against unmitigated BF16 throughput, and not against tensor-core GEMMs, which this design deliberately never uses (its IEEE-FMA Triton GEMM reaches about 0.85× cuBLAS at prefill).
6. One root, two axes, two different fixes
Re-reading both papers together makes the structure clean. Floating-point non-associativity is the invariant root; the papers differ in which degree of freedom orders the additions. In the precision paper, hardware and kernels are fixed and the formats differ — BF16's 7-bit mantissa (machine epsilon about 3.9 · 10−3) versus FP16's 10 bits (about 4.9 · 10−4), an 8× gap — so the same weights produce different logits at every step. In the GEMM paper the formats are fixed and the kernel selection differs per device, so the same weights, same format, same software produce different reduction orders on different GPUs. Both channels terminate in the same event: a directional perturbation exceeding a small top-two margin at the decision, amplified by the lm_head projection, cascading through autoregressive conditioning.
The fixes differ accordingly, and neither is "more precision everywhere." The precision paper repairs where (only the lm_head, only low-margin steps) — and shows repairing more broadly is worse. The GEMM paper repairs how the order is fixed (make the reduction order a pure function of problem shape) — and shows that merely shrinking the noise to FP32 scale leaves a statistical, not constructed, guarantee. Each is the honest response to its axis: gating on the margin exploits the bimodality the margin distribution gives you; pinning the reduction order exploits the fact that IEEE-754 FMA is bit-specified for given operands. Between them they also bracket the limits: the head-scoped repair fails at batch 8 and under end-to-end FP8; the fixed-order GEMM leaves attention, normalization, and sampler kernels unpinned (H100's attention path still shifted margins at a median 7.9 · 10−6 nats, with 9.6% of positions bitwise equal versus 98.4% between A100 and L40S)2.
7. Anti-hype: what neither paper claims
Both papers are careful about scope, and the care is the useful part.
Scale. The precision paper's evidence runs 1.1B to 7B across six models, plus a divergence-only probe at 12B (Mistral-Nemo). The full intervention comparison is not evaluated at 12B or 70B+ — the paper says so verbatim in its limitations. Whether 400B-class frontier models show the same margin-event structure, or the same +22pp-style recoverability, is an extrapolation the paper does not make. What does transfer by construction is the mechanism: the exact flip condition is format-pair-agnostic, and the paper demonstrates the same low-margin pattern in head-isolated FP8 and FP16-vs-FP32 settings.
"More precision makes it worse" needs the prediction framing. The counterintuitive result is specifically extending FP32 scope beyond the lm_head at safe steps: prediction 5 of the paper's five pre-registered predictions, confirmed as 63% → 44% on TinyLlama (RMSNorm+lm_head) and replicated on DeepSeek-R1-7B at τ = 0. It is not a claim that FP32 inference is worse than BF16 — the FP32 oracle is 100% agreement, and FP16 is the better FP32 proxy (FP16-alone agrees with FP32 on 90% of sequences versus 42% for BF16 on TinyLlama). The lesson is about selectivity: untargeted re-arithmetic at agreeing steps creates divergence that was not there.
The mitigation is partial, and its boundaries are the finding. +22–36pp on A10G means 55–67% agreement, not 100%. The residual has a measured decomposition: roughly 12pp is the BF16-vs-FP16 weight-truncation floor (different stored mantissas), irreducible without a shared storage format; batch-8 and end-to-end-FP8 regimes are body-dominated and need body-scope composition; accuracy is preserved within tested margins (TOST non-inferiority on Qwen2.5-3B at n = 600, p = 0.009) but Mistral-7B shows an unmatched −3 to −5pp point estimate that n = 100 cannot resolve. Applications needing bit-exact replay should pin the full numerical configuration or use an end-to-end higher-precision reference — the paper's own recommendation1.
Bitwise kernel determinism is scoped to the linear layers. The GEMM paper's guarantee covers the GEMMs (the dominant FLOPs and the identified divergence channel) — attention, rotary embeddings, normalization, and the sampler remain vendor kernels with architecture-dependent paths, and the one surviving end-to-end divergence cell is exactly there. Cross-vendor portability, where software differs and not just hardware, is stated as an open problem, and the 1.17–3.1× claim is specifically versus the FP32-compute state of the art, not versus unmitigated BF162.
The deployment question both papers converge on: stop treating "greedy, temperature 0, seed fixed" as a reproducibility contract. Record the full numerical configuration — format pair, batch size, kernel choices, architecture — or engineer the order of additions to be invariant. Determinism is not a property greedy decoding gives you; it is a property you build, either by gating on margins or by pinning reduction orders, and each approach comes with a measured map of where it stops working.
Footnotes
Footnotes
-
Du, Gaoyuan; Khan, Anam Nawaz; Zhou, Rex; Liu, Xiaoyang; Chakrabarti, Deepayan; Suya, Fnu; Li, Xueping — Greedy Decoding Is Not Precision-Invariant: Cross-Precision Output Divergence in LLM Inference, arXiv:2609.26621, cs.LG, submitted September 22, 2026, accepted at TMLR (September 2026). Full PDF v1 verified at primary: six models 1.1B-7B plus 12.2B divergence-only probe; GSM8K/HumanEval/MBPP; TinyLlama EAR 41/36/30%, 49-100% divergence range, mean cascade 34 tokens, Qwen2.5-3B 19% correctness flips at 1pp aggregate gap; Prop. 1 exact directional flip condition; sigma_z about 0.026, directional scale about 0.037, directional rule 88-89% step accuracy; body L2 1.022 vs 1.135 indistinguishable, margin 0.039 vs 5.733; Diverged median margin 0 with 30/59 exact BF16 ties vs Agreed 2.69; L*alpha/sqrt(d) about 0.024, 97.6% head contribution; five predictions all confirmed (A zero, sharpening zero, K = 2 suffices, full matches top-K, scope extension 63% to 44%, ungated tau = 0 at 54-60% and 2.5x latency); gated C: tau = 1e-3, 0.7-1.4% trigger, +22pp GSM8K 41 to 63%, +1.4% latency +11% memory, fires 30/100 converts 22; tiers +36/+8/+10/+3/0; L4/A100 +12 to +21pp, Qwen2.5-3B +3 A10G vs +20 L4; batch 17/7/7/0pp at bs 1/2/4/8, +5pp composed with global FP32 at bs 8; head-isolated FP8 rescaled gate 100% TinyLlama / +56pp Qwen at tau = 0.25; end-to-end FP8 (155 modules) +1pp; weight-truncation floor about 12pp (FP32-compute pair at 0.88); TOST p = 0.009 on Qwen2.5-3B n = 600; Mistral-7B -3 to -5pp point estimates at McNemar p = 0.125/0.250; 100% within-format self-consistency under deterministic flags: https://arxiv.org/abs/2609.26621 ↩ ↩2 ↩3 ↩4 ↩5 ↩6 ↩7 ↩8 ↩9 ↩10 ↩11 ↩12 ↩13
-
Cooper, Liam; Jeong, Shinnung; Jeon, Hyeran; Young, Jeffrey; Kim, Hyesoon — Accelerating the Mitigation of LLM Inference Nondeterminism Across GPU Architectures, arXiv:2609.25624, cs.AR/cs.LG, submitted September 22, 2026 (Georgia Tech; UC Merced; source at github.com/lpc97667/rf; the system's stylized name does not survive PDF text extraction and is not reproduced here). Full PDF v1 verified at primary: R1-R4 rules as stated in text; reduction order a pure function of (M,N,K), decode bucket M up to 64; bitwise-identical linear outputs across A100/L40S/H100 on 7/7 probe shapes vs 0/7 for cast-then-cuBLAS (A100-L40S disagree on all seven, 88-99% of elements differing); unmitigated BF16 30.81-100% cross-GPU divergence; FP32-compute residual up to 0.51% on 5 of 12 cells; proposed kernels zero on 11 of 12 cells, the twelfth 0.08% from an unpinned H100 layer-0 attention kernel at an exact FP32 tie; margin tail near-linear below 0.1 nats, rho 0.06-0.17 per nat, no safe precision short of exactness; 1.17-3.1x end-to-end vs FP32-compute SOTA in all 36 configurations (1.17-1.43x A100/H100, 1.6-3.1x L40S); weight traffic halved; 1-1.4 GiB savings on untied-embedding models; Llama-3.2-3B, Qwen3-4B, DeepSeek-R1-Distill-Llama-8B (Qwen3-14B memory profiling), GSM8K/MATH500/AIME24/GPQA-Diamond, vLLM 0.8, batch 32; prefill IEEE-FMA GEMM about 0.85x cuBLAS; cross-vendor portability open: https://arxiv.org/abs/2609.25624 ↩ ↩2 ↩3 ↩4 ↩5 ↩6 ↩7 ↩8 ↩9
-
Yuan, Jiayi; Li, Hao; Ding, Xinheng; Xie, Wenya; Li, Yu-Jhe; Zhao, Wentian; Wan, Kun; Shi, Jing; Hu, Xia; Liu, Zirui — Understanding and Mitigating Numerical Sources of Nondeterminism in LLM Inference, NeurIPS 2025 (oral): BF16 storage with per-layer just-in-time FP32 upcast, divergence below 3.4% of problems as characterized at arXiv:2609.25624 (the primary for this guide). Not independently re-verified at the NeurIPS source. ↩