perf(bz): word-sized modulus arithmetic for the lattice tier's lift and CLD stages

[kim-em]

Current state. #8543 deliverable 1 landed (#8670): the lattice-tier loop's .degenerate arm now reads the all-ones certificate off the same CLD lattice build the classifier already ran, byte-identical, 1.27×/1.62×/1.82× at SD4/SD5/SD6. This issue is deliverable 3 of #8543, scoped as its own work item: word-sized modulus arithmetic in the lattice tier's lift and CLD stages.

Where the time goes. Recovery-pipeline sub-stage profile at the first k ≥ floor (carica-class machine, hex_lattice_spike-style):

input	lift	CLD build	LLL+cut	RREF
SD4 deg16	1.9ms	0.39ms	3.04ms (57%)	0.05ms
SD5 deg32	15.4ms	9.3ms	75.8ms (75%)	0.05ms
SD6 deg64	91ms	98.9ms	2161ms (92%)	0.06ms

The LLL reduction dominates and is integer Gram-Schmidt (not modular), so ZMod64/Montgomery does not apply to it — that is out of scope here (a separate word-sized/float LLL path would be its own issue). The word-sizable modular hot spots are the lift and CLD build, and they matter most at the SD4 / r=16 crossover, where the mod-pᵃ modulus fits 64 bits. At SD5/SD6 the modulus exceeds 2⁶⁴ and the big-integer path must stay; the win there is negligible (and those cases are already "won" versus the Isabelle comparator).

Deliverables.

1. Word-sized poly-mod-pᵃ layer. No polynomial type over a composite word-sized modulus exists yet (FpPoly is prime-only). Scalar Montgomery over an odd word-sized modulus does exist — HexArith.MontCtx needs only p % 2 = 1, and an odd prime power pᵃ < 2⁶⁴ qualifies. Build a small WordPolyMod-style layer (UInt64 residue arrays mod m, Montgomery multiply, monic division, reduce) with proven arithmetic. This is the landable first unit.
2. Route cldQuotientMod through it under a guard. cldQuotientMod f g p a = reduceModPow (f · g') p a → DensePoly.divMod _ g → reduce, all poly arithmetic mod pᵃ. Dispatch on Odd p ∧ pᵃ < 2⁶⁴; big-integer fallback otherwise. Prove cldQuotientMod = cldQuotientModWord under the guard so the CLD columns, indicators, and conformance stay byte-identical.
3. Same pattern for the lift. Route multifactorLiftQuadratic / quadraticHenselStep's mod-p^(2^i) poly arithmetic through the word-sized layer while p^(2^i) < 2⁶⁴, big-integer continuation past that. (Note the coupling: this touches HexHensel, which is more heavily proof-coupled than the CLD stage — land deliverable 2's data separately, and prefer the CLD stage first.)

Verification. Conformance byte-identical (hexbz_emit_fixtures, bz_flint.py, bz_trace_gate.py); hex_lattice_spike core ratios recorded before/after per deliverable, with special attention to SD4 (r=16); runFactorLatticeAdvSwinnertonDyerSD{3,4}Checksum stay inside the verify budget; no change to factorFast semantics; no sorry/axiom; the CI-gated HexBerlekampZassenhausMathlib layer stays green.

Context. Descends from #8543 deliverable 3 (that umbrella issue is closed; deliverable 1 landed as #8670, deliverable 2 was deprioritized after measurement showed SD5/SD6 do a single lift+recover step, and deliverable 4 — re-measuring the classical/lattice crossover and revisiting the dispatch threshold — should follow once this lands, since it is what actually moves the r=16 seam). The scalar Montgomery infrastructure is in HexArith/Montgomery/Context.lean and HexModArith/HotLoop.lean.

🤖 Prepared with Claude Code

[kim-em]

Before building the WordPolyMod layer I measured the actual operating point, and the picture differs materially from the issue's framing. Posting the evidence and a re-scope recommendation rather than landing a large proof-only layer against the stated premise.

The SD4 modulus is ~2^45, not near 2^64

The issue says the win "matters most at the SD4 / r=16 crossover, where the mod-p^a modulus fits 64 bits." At the actual recovery point that #8395 early-termination stops at, the modulus is much smaller than the precision cap suggests. Instrumenting factorLattice on the hex_lattice_spike inputs:

input	prime p	exponent a at recovery	modulus p^a	bits
SD4 deg16	11	13	11^13 = 34522712143931	~45
SD5 deg32	19	23	19^23	~95

(The full precision cap B = latticePrecisionCap f is a coefficient magnitude, not the exponent; the exponent is precisionForCoeffBound B p = ceilLogP p (2B+1), and the doubling schedule stops at the first level clearing bhksRecoveryFloor core, well short of the cap. For SD4 the cap-precision exponent would be 398, i.e. p^a ~1375 bits, but recovery lands at exponent 13.)

So the word-sized path applies to SD4 (11^13 < 2^64) but the story is "a targeted SD4 CLD/lift optimization at ~45-bit moduli," not "needed for 64-bit precision pressure."

The win is real: 1.8x on CLD arithmetic, byte-identical

I wrote an unproven native-word Montgomery prototype of cldQuotientMod (residues in Montgomery form over Array UInt64, built on HexArith.MontCtx.mulMont, full-range modular add/sub, monic long division) and compared it against the current Int/DensePoly path at m = 11^13, SD4-shaped inputs (f deg 16, 8 divisors of deg 2), 20000 iterations x8 factors:

modulus m = 11^13 = 34522712143931
correctness: OK (byte-identical)
  ref   (Int/DensePoly): 3022.6 ms
  proto (Montgomery)   : 1679.4 ms
  speedup: 1.80x

The output is byte-identical and the prototype is 1.8x faster. The reason it wins here where the packed-UInt64 experiment in HexPolyFp/Packed.lean (#8642) measured neutral: at p < 2^31 both paths do a single-word % p, so representation does not matter. At ~45-bit residues the products are ~90 bits, which exceed Lean's unboxed-Int boundary (2^62) and heap-allocate as GMP temporaries on every coefficient multiply. Native Montgomery keeps them in words (mulFull + REDC, no allocation). That is the actual lever, and it is genuinely different from the #8642 regime.

But the end-to-end payout is bounded, and deliverable 2 is the real cost

• End-to-end at SD4, the issue's own sub-stage profile has LLL at 57% (3.04ms of ~5.4ms) and it is out of scope. Halving lift+CLD (2.3ms) saves ~1.15ms, i.e. roughly 18% at SD4.
• SD5/SD6 moduli exceed 2^64, so the word path does not apply there, and those cases are already won versus the comparator. So the realistic scope is "SD4-shaped inputs only," not a general lattice-tier unlock.
• Deliverable 1 (the proven WordPolyMod layer) delivers zero runtime benefit until deliverable 2 wires it in. As a standalone PR it is dead code plus maintenance surface.
• Deliverable 2's byte-identical equality proof cldQuotientMod = cldQuotientModWord is the dominant cost and risk, not a wiring step. There is no DensePoly.map/ring-hom transport in the repo today. It needs: a word ring WordMod m (odd m < 2^64) as a CommRing with Montgomery-form representation, a ring hom to ZMod m, and proofs that DensePoly mul/sub/reduceModPow and monic divMod all commute with the coefficient map. Under the no-sorry/no-axiom rule that is a large from-scratch formalization, and HexBerlekampZassenhausMathlib/CLDColumnBound.lean depends on the current definitional shape of cldQuotientMod (a rfl at line ~422), so the rewrite has to preserve or replace that.

Recommendation

I did not land a deliverable-1-only PR (dead code) or bash out the full transport proof in one session. Suggest re-scoping before the next code unit, so the first PR retires proof risk rather than adding unused machinery. Concretely, pick one of:

1. Re-confirm that a proof-only foundational WordPolyMod PR is wanted given ~18%-at-SD4-only economics, with a named owner for the deliverable-2 transport proof. Then I build it.
2. Feasibility-first: make the first PR the single hardest bridge, a proven DensePoly.map-commutes-with-monic-divMod (and with reduceModPow) transport, in the exact shape cldQuotientMod needs. If that is clean, deliverable 1 is justified; if it is ugly, we have saved building a layer that cannot be connected affordably.
3. Fold this into deliverable 4 (re-measuring the classical/lattice crossover): the r=16 seam is what actually moves the win, and a ~18%-of-SD4 CLD/lift speedup may be better spent after the crossover threshold is revisited.

The prototype (unproven, throwaway) and the instrumentation are reproducible; happy to attach the full source. Leaving this claimable for a re-scoped version per the "directives are hypotheses" doctrine.

🤖 Prepared with Claude Code

[kim-em]

Reproduction artifacts for the measurement above (throwaway, unproven, not committed).

To reproduce: drop CldProto.lean at bench/HexBench/CldProto.lean, add a lean_exe hex_cld_proto (srcDir bench, root HexBench.CldProto, moreLinkArgs := #["-lgmp"]), lake build hex_cld_proto, run it. It checks byte-identity against Hex.cldQuotientMod and times both paths.

CldProto.lean (native-Montgomery prototype of cldQuotientMod)

/-
Copyright (c) 2026 Lean FRO, LLC. All rights reserved.
Released under Apache 2.0 license as described in the file LICENSE.
Authors: Kim Morrison
-/

import HexBerlekampZassenhaus
import HexArith.Montgomery.Context

/-!
Throwaway premise-check microbenchmark for issue #8691: does a native-word
Montgomery `cldQuotientMod` beat the current `Int`/`DensePoly` path at the
*actual* SD4 recovery modulus (p^a = 11^13 ~ 2^45)?  Not proven, not wired;
exists only to measure the win before committing to a large proven layer.
-/

open Hex

/-- The reference CLD quotient (exactly `Hex.cldQuotientMod`). -/
def refCld (f g : ZPoly) (p a : Nat) : ZPoly := cldQuotientMod f g p a

namespace Proto

/-- Modular add for residues below `m`, `m < 2^63` so `a + b` never overflows. -/
@[inline] def addMod (m a b : UInt64) : UInt64 :=
  let s := a + b
  if s ≥ m then s - m else s

/-- Modular subtract for residues below `m`. -/
@[inline] def subMod (m a b : UInt64) : UInt64 :=
  if a ≥ b then a - b else (m - b) + a

/-- Reduce an `Int` coefficient to `[0, m)`, then into Montgomery form. -/
@[inline] def ofInt (mUInt : UInt64) (ctx : MontCtx mUInt) (m : Nat) (c : Int) : UInt64 :=
  let r : Nat := (c % (m : Int)).toNat % m
  ctx.toMont (UInt64.ofNat r)

/-- Montgomery-form coefficient array (padded to `len`) for an `Int` poly mod `m`. -/
def toMontArr (mUInt : UInt64) (ctx : MontCtx mUInt) (m : Nat) (f : ZPoly) (len : Nat) :
    Array UInt64 :=
  (Array.range len).map (fun i => ofInt mUInt ctx m (f.coeff i))

/-- Derivative in Montgomery form: `g'[i] = (i+1) * g[i+1]`. -/
def derivMont (mUInt : UInt64) (ctx : MontCtx mUInt) (g : Array UInt64) : Array UInt64 :=
  if g.size ≤ 1 then #[] else
  (Array.range (g.size - 1)).map (fun i =>
    ctx.mulMont (g[i+1]!) (ctx.toMont (UInt64.ofNat ((i+1) % mUInt.toNat))))

/-- Schoolbook convolution in Montgomery form, result length `a.size + b.size - 1`. -/
def mulMontArr (mUInt : UInt64) (ctx : MontCtx mUInt) (a b : Array UInt64) : Array UInt64 :=
  if a.size == 0 || b.size == 0 then #[] else
  let n := a.size + b.size - 1
  Id.run do
    let mut out := (Array.range n).map (fun _ => (0 : UInt64))
    for i in [0:a.size] do
      for j in [0:b.size] do
        let prod := ctx.mulMont (a[i]!) (b[j]!)
        out := out.set! (i+j) (addMod mUInt (out[i+j]!) prod)
    return out

/-- Trim trailing zeros. -/
def trim (a : Array UInt64) : Array UInt64 :=
  let rec go (n : Nat) : Nat :=
    match n with
    | 0 => 0
    | n+1 => if a[n]! == 0 then go n else n+1
  a.extract 0 (go a.size)

/-- Monic long division in Montgomery form (leading true-coeff of `g` is 1, so
the pivot needs no inverse): return the quotient (mont form). -/
def divMonicMont (mUInt : UInt64) (ctx : MontCtx mUInt) (num g : Array UInt64) : Array UInt64 :=
  let num := trim num
  let g := trim g
  let dg := g.size - 1
  if num.size < g.size then #[] else
  let qlen := num.size - dg
  Id.run do
    let mut rem := num
    let mut quot := (Array.range qlen).map (fun _ => (0 : UInt64))
    for step in [0:qlen] do
      let k := num.size - 1 - step
      let qi := k - dg
      let c := rem[k]!
      quot := quot.set! qi c
      for t in [0:g.size] do
        let sub := ctx.mulMont c (g[t]!)
        rem := rem.set! (qi + t) (subMod mUInt (rem[qi+t]!) sub)
    return quot

/-- Native prototype of `cldQuotientMod`: reduce(f·g') / g, all mod `m`. -/
def cldMont (mUInt : UInt64) (ctx : MontCtx mUInt) (m : Nat) (f g : ZPoly) : Array Nat :=
  let degF := f.degree?.getD 0
  let degG := g.degree?.getD 0
  let fA := toMontArr mUInt ctx m f (degF + 1)
  let gA := toMontArr mUInt ctx m g (degG + 1)
  let gD := derivMont mUInt ctx gA
  let num := mulMontArr mUInt ctx fA gD
  let quot := divMonicMont mUInt ctx num gA
  quot.map (fun c => (ctx.fromMont c).toNat)

end Proto

/-- Deterministic pseudo-random monic ZPoly of degree `d`, coeffs mod 2^40. -/
def randMonic (seed d : Nat) : ZPoly :=
  let coeffs := (Array.range d).map (fun i =>
    Int.ofNat (((seed * 6364136223846793005 + i * 1442695040888963407 + 12345)
      % 1099511627776)))
  DensePoly.ofCoeffs (coeffs.push 1)

def main (_args : List String) : IO Unit := do
  let out ← IO.getStdout
  let p := 11
  let a := 13
  let m := p ^ a
  let mUInt := UInt64.ofNat m
  IO.println s!"modulus m = {p}^{a} = {m}  (mUInt={mUInt})"
  let f := randMonic 7 16
  let gs := (Array.range 8).map (fun i => randMonic (100 + i) 2)
  have hodd : mUInt % 2 = 1 := by decide
  let ctx := MontCtx.mk mUInt hodd
  let mut ok := true
  for g in gs do
    let refQ := (refCld f g p a).toArray.map (fun c => c.toNat)
    let protoQ := Proto.cldMont mUInt ctx m f g
    let refQt := (refQ.toList.reverse.dropWhile (· == 0)).reverse
    let protoQt := (protoQ.toList.reverse.dropWhile (· == 0)).reverse
    if refQt != protoQt then
      ok := false
      IO.println s!"MISMATCH ref={refQt} proto={protoQt}"
  let status := if ok then "OK (byte-identical)" else "FAIL"
  IO.println s!"correctness: {status}"
  out.flush
  let iters := 20000
  let t0 ← IO.monoNanosNow
  let (sinkRef, _) ← IO.lazyPure (fun _ => Id.run do
    let mut s : UInt64 := 0
    for _ in [0:iters] do
      for g in gs do
        s := s + UInt64.ofNat ((refCld f g p a).coeff 0).natAbs
    (s, ()))
  let t1 ← IO.monoNanosNow
  let (sinkProto, _) ← IO.lazyPure (fun _ => Id.run do
    let mut s : UInt64 := 0
    for _ in [0:iters] do
      for g in gs do
        s := s + UInt64.ofNat ((Proto.cldMont mUInt ctx m f g)[0]!)
    (s, ()))
  let t2 ← IO.monoNanosNow
  let refMs := (t1 - t0).toFloat / 1.0e6
  let protoMs := (t2 - t1).toFloat / 1.0e6
  IO.println s!"iters={iters} (x8 factors each call)"
  IO.println s!"  ref   (Int/DensePoly): {refMs} ms  (sink={sinkRef})"
  IO.println s!"  proto (Montgomery)   : {protoMs} ms  (sink={sinkProto})"
  IO.println s!"  speedup: {refMs / protoMs}x"
  out.flush

[kim-em]

Re-taking this: building it via @[csimp] (add cldQuotientModWord + guarded cldQuotientModFast, prove cldQuotientMod = cldQuotientModFast) so existing proofs and the CLDColumnBound.lean rfl shape stay untouched, with the equality proved in the *Mathlib layer over ZMod (p^a) via monic-division uniqueness. Starting with the word kernel + arithmetic spec.

[kim-em]

Phases 1 and 1b landed in #8697 "feat(bz): word-sized Montgomery residue ring and CLD kernel": the proven WordMod odd-m < 2^64 Montgomery ring, and cldQuotientModWord? with conformance #guards pinning it byte-identical to cldQuotientMod (moduli 7^5, 11^13, 5^9) and declining on 11^400. No sorry/axiom; whole graph green.

Phase 2 (the byte-identical proof, its own PR next): WordMod.toZMod : WordMod ctx → ZMod (p^a) bridge + CommRing via Function.Injective.commRing + ≃+* ZMod; transport both quotients into Polynomial (ZMod (p^a)) through DensePoly.toPolynomial and prove dense monic divMod corresponds to Polynomial.divByMonic; conclude cldQuotientMod = cldQuotientModWord? under guard. Phase 3 then wires it with @[csimp] and re-measures SD4.