diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 7fcfdffed..6dab3995c 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -162,8 +162,8 @@ jobs:
       # fresh against the restored libs, which is cheap.
       - name: Define the hex-dev cache + build target sets
         run: |
-          echo "HEX_LIB_TARGETS=HexBasic HexTruncatedSeries HexTruncatedSeriesMathlib HexArith HexPoly HexMvPoly HexModArith HexGF2 HexPolyZ HexRoots HexResultant HexInterval HexIntervalExperiment HexIntervalMathlib HexIntervalMathlibExperiment HexIntervalReplayProbe HexIntervalMathlibReplayProbe HexRealRootsMathlibReplayProbe HexRCFProofProbe HexPolyFp HexGFqRing HexGFqField HexBerlekamp HexHensel HexConway HexGFq HexPrimality HexIntFactor HexPrimalityKernelProbe HexIntFactorKernelProbe HexMvGcdKernelProbe HexBerlekampZassenhaus HexRealRoots HexMatrix HexRowReduce HexDeterminant HexBareiss HexHermite HexSmith HexCharPoly HexMinPoly HexPolySmith HexGramSchmidt HexLLL HexMatrixMathlib HexHermiteMathlib HexSmithMathlib HexSmithTests HexCharPolyMathlib HexMinPolyMathlib HexPolySmithMathlib HexGramSchmidtMathlib HexLLLMathlib HexBerlekampZassenhausMathlib HexBerlekampZassenhausMathlibProofProbe HexBerlekampMathlibProofProbe HexPrimalityMathlib HexIntFactorMathlib HexPolyFpMathlib HexFactorizationModules HexRealRootsMathlib HexGF2Mathlib HexGFqMathlib HexMvPolyMathlib HexMvPolyMathlibProofProbe HexRCF HexRCFTests HexRootsMathlib HexResultantMathlib HexNumberField HexNumberFieldMathlib HexNumberFieldTower HexNumberFieldTowerMathlib HexReleaseTests HexReleaseExamples HexAggregateCheck" >> "$GITHUB_ENV"
-          echo "HEX_EXE_TARGETS=hextruncatedseries_bench hexarith_bench hexpoly_bench hexpolysmith_bench hexmvpoly_bench hexmvgcd_bench hexsparsepoly_bench hexpolyz_bench hexpolyfp_bench hexmodarith_bench hexmodular_bench hexgf2_bench hexgfqring_bench hexgfqfield_bench hexgfq_bench hexhensel_bench hexprimality_bench hexintfactor_bench hexberlekamp_bench hexbz_bench hexconway_bench hexmatrix_bench hexstrassen_compare hexdeterminant_bench hexbareiss_bench hexcharpoly_bench hexminpoly_bench hexgramschmidt_bench hexhermite_bench hexsmith_bench hexrealroots_bench hexrcf_bench hexroots_bench hexresultant_bench hexnumberfield_bench hexnumberfieldtower_bench hexinterval_decision_bench hexroots_demo hex_interval_representation_spike hex_interval_center_spike hex_interval_scale_spike hex_interval_scheduler_spike hex_interval_policy_frontier_spike hex_arith_floor hexlll_bench hexlll_gram_bench hexlll_external_reduction hexbz_factor_service" >> "$GITHUB_ENV"
+          echo "HEX_LIB_TARGETS=HexBasic HexTruncatedSeries HexTruncatedSeriesMathlib HexArith HexPoly HexPolyFast HexMvPoly HexModArith HexGF2 HexPolyZ HexRoots HexResultant HexInterval HexIntervalExperiment HexIntervalMathlib HexIntervalMathlibExperiment HexIntervalReplayProbe HexIntervalMathlibReplayProbe HexRealRootsMathlibReplayProbe HexRCFProofProbe HexPolyFp HexGFqRing HexGFqField HexBerlekamp HexHensel HexConway HexGFq HexPrimality HexIntFactor HexPrimalityKernelProbe HexIntFactorKernelProbe HexMvGcdKernelProbe HexBerlekampZassenhaus HexRealRoots HexMatrix HexRowReduce HexDeterminant HexBareiss HexHermite HexSmith HexCharPoly HexMinPoly HexPolySmith HexGramSchmidt HexLLL HexMatrixMathlib HexHermiteMathlib HexSmithMathlib HexSmithTests HexCharPolyMathlib HexMinPolyMathlib HexPolySmithMathlib HexGramSchmidtMathlib HexLLLMathlib HexBerlekampZassenhausMathlib HexBerlekampZassenhausMathlibProofProbe HexBerlekampMathlibProofProbe HexPrimalityMathlib HexIntFactorMathlib HexPolyFpMathlib HexFactorizationModules HexRealRootsMathlib HexGF2Mathlib HexGFqMathlib HexMvPolyMathlib HexMvPolyMathlibProofProbe HexRCF HexRCFTests HexRootsMathlib HexResultantMathlib HexNumberField HexNumberFieldMathlib HexNumberFieldTower HexNumberFieldTowerMathlib HexReleaseTests HexReleaseExamples HexAggregateCheck" >> "$GITHUB_ENV"
+          echo "HEX_EXE_TARGETS=hextruncatedseries_bench hexarith_bench hexpoly_bench hexpolyfast_bench hexpolysmith_bench hexmvpoly_bench hexmvgcd_bench hexsparsepoly_bench hexpolyz_bench hexpolyfp_bench hexmodarith_bench hexmodular_bench hexgf2_bench hexgfqring_bench hexgfqfield_bench hexgfq_bench hexhensel_bench hexprimality_bench hexintfactor_bench hexberlekamp_bench hexbz_bench hexconway_bench hexmatrix_bench hexstrassen_compare hexdeterminant_bench hexbareiss_bench hexcharpoly_bench hexminpoly_bench hexgramschmidt_bench hexhermite_bench hexsmith_bench hexrealroots_bench hexrcf_bench hexroots_bench hexresultant_bench hexnumberfield_bench hexnumberfieldtower_bench hexinterval_decision_bench hexroots_demo hex_interval_representation_spike hex_interval_center_spike hex_interval_scale_spike hex_interval_scheduler_spike hex_interval_policy_frontier_spike hex_arith_floor hexlll_bench hexlll_gram_bench hexlll_external_reduction hexbz_factor_service" >> "$GITHUB_ENV"
       # Shared build. The libraries, bench exes, conformance #guard drivers, and
       # emit-fixture exes are all elaborated here so the two verification tails
       # below only *run* things, never rebuild them -- which is what lets the
@@ -222,7 +222,7 @@ jobs:
         run: |
           bash scripts/ci/check_bench_verify_budget.sh \
             hextruncatedseries_bench \
-            hexarith_bench hexpoly_bench hexpolysmith_bench \
+            hexarith_bench hexpoly_bench hexpolyfast_bench hexpolysmith_bench \
             hexmvpoly_bench hexmvgcd_bench \
             hexpolyz_bench \
             hexpolyzgcd_bench hexsparsepoly_bench hexpolyfp_bench \
diff --git a/HexPolyFast.lean b/HexPolyFast.lean
new file mode 100644
index 000000000..119c6f789
--- /dev/null
+++ b/HexPolyFast.lean
@@ -0,0 +1,19 @@
+/-
+Copyright (c) 2026 Lean FRO, LLC. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Kim Morrison
+-/
+
+module
+
+public import HexPolyFast.Plan
+public import HexPolyFast.Karatsuba
+public import HexPolyFast.Cyclic
+public import HexPolyFast.Reverse
+public import HexPolyFast.Reciprocal
+public import HexPolyFast.Division
+public import HexPolyFast.Tree
+
+public section
+
+/-! Fast proof-carrying dense-polynomial algorithms. -/
diff --git a/HexPolyFast/Cyclic.lean b/HexPolyFast/Cyclic.lean
new file mode 100644
index 000000000..e15531ce2
--- /dev/null
+++ b/HexPolyFast/Cyclic.lean
@@ -0,0 +1,104 @@
+/-
+Copyright (c) 2026 Lean FRO, LLC. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Kim Morrison
+-/
+
+module
+
+public import HexPolyFast.Plan
+
+public section
+
+/-!
+# Cyclic and negacyclic products
+
+The reference kernels fold an ordinary planned product by congruence classes
+of exponents. They allocate exactly the requested `n` coefficient slots and
+therefore also serve as independent comparators for transform-based plans.
+-/
+
+namespace Hex.DensePoly
+
+universe u
+
+attribute [local instance 1000] Lean.Grind.Semiring.ofNat
+
+variable {R : Type u} [DecidableEq R] [Lean.Grind.CommRing R]
+
+/-- Sum coefficients of `p` whose exponents are congruent to `k` modulo `n`. -/
+def cyclicCoeff (n k : Nat) (p : DensePoly R) : R :=
+  (List.range ((p.size + n - 1) / n)).foldl
+    (fun acc t => acc + p.coeff (k + t * n)) 0
+
+/-- Alternating sum of coefficients in one congruence class, implementing
+reduction by `x^n = -1`. -/
+def negacyclicCoeff (n k : Nat) (p : DensePoly R) : R :=
+  (List.range ((p.size + n - 1) / n)).foldl
+    (fun acc t =>
+      if t % 2 = 0 then acc + p.coeff (k + t * n)
+      else acc - p.coeff (k + t * n)) 0
+
+/-- Cyclic convolution of positive length `n`. -/
+def mulCyclic (plan : MulPlan R) (n : Nat) (_hn : 0 < n)
+    (a b : DensePoly R) : DensePoly R :=
+  let p := mulWith plan a b
+  ofList ((List.range n).map fun k => cyclicCoeff n k p)
+
+/-- Negacyclic convolution of positive length `n`. -/
+def mulNegacyclic (plan : MulPlan R) (n : Nat) (_hn : 0 < n)
+    (a b : DensePoly R) : DensePoly R :=
+  let p := mulWith plan a b
+  ofList ((List.range n).map fun k => negacyclicCoeff n k p)
+
+/-- Checked cyclic convolution; length zero has no quotient-ring meaning. -/
+def mulCyclic? (plan : MulPlan R) (n : Nat) (a b : DensePoly R) :
+    Option (DensePoly R) :=
+  if hn : 0 < n then some (mulCyclic plan n hn a b) else none
+
+/-- Checked negacyclic convolution; length zero has no quotient-ring meaning. -/
+def mulNegacyclic? (plan : MulPlan R) (n : Nat) (a b : DensePoly R) :
+    Option (DensePoly R) :=
+  if hn : 0 < n then some (mulNegacyclic plan n hn a b) else none
+
+/-- Cyclic products contain at most `n` stored coefficients. -/
+theorem size_mulCyclic_le (plan : MulPlan R) (n : Nat) (hn : 0 < n)
+    (a b : DensePoly R) : (mulCyclic plan n hn a b).size ≤ n := by
+  unfold mulCyclic
+  exact Nat.le_trans (size_ofList_le _) (by simp)
+
+/-- Negacyclic products contain at most `n` stored coefficients. -/
+theorem size_mulNegacyclic_le (plan : MulPlan R) (n : Nat) (hn : 0 < n)
+    (a b : DensePoly R) : (mulNegacyclic plan n hn a b).size ≤ n := by
+  unfold mulNegacyclic
+  exact Nat.le_trans (size_ofList_le _) (by simp)
+
+/-- Coefficient description of cyclic folding. -/
+theorem coeff_mulCyclic (plan : MulPlan R) (n : Nat) (hn : 0 < n)
+    (a b : DensePoly R) (i : Nat) :
+    (mulCyclic plan n hn a b).coeff i =
+      if i < n then cyclicCoeff n i (a * b) else 0 := by
+  unfold mulCyclic
+  rw [coeff_ofList]
+  by_cases hi : i < n
+  · rw [_root_.ite_eq_left hi]
+    simp [List.getD, hi, mulWith_eq]
+  · rw [_root_.ite_eq_right hi, List.getD_eq_getElem?_getD]
+    simp [hi]
+    rfl
+
+/-- Coefficient description of negacyclic folding. -/
+theorem coeff_mulNegacyclic (plan : MulPlan R) (n : Nat) (hn : 0 < n)
+    (a b : DensePoly R) (i : Nat) :
+    (mulNegacyclic plan n hn a b).coeff i =
+      if i < n then negacyclicCoeff n i (a * b) else 0 := by
+  unfold mulNegacyclic
+  rw [coeff_ofList]
+  by_cases hi : i < n
+  · rw [_root_.ite_eq_left hi]
+    simp [List.getD, hi, mulWith_eq]
+  · rw [_root_.ite_eq_right hi, List.getD_eq_getElem?_getD]
+    simp [hi]
+    rfl
+
+end Hex.DensePoly
diff --git a/HexPolyFast/Division.lean b/HexPolyFast/Division.lean
new file mode 100644
index 000000000..54711dc28
--- /dev/null
+++ b/HexPolyFast/Division.lean
@@ -0,0 +1,570 @@
+/-
+Copyright (c) 2026 Lean FRO, LLC. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Kim Morrison
+-/
+
+module
+
+public import HexPolyFast.Reciprocal
+public import HexPoly.Euclid.MonicUnique
+
+public section
+
+/-!
+Cached reciprocal division.
+
+The executable quotient reverses the dividend, multiplies its low prefix by
+the cached reciprocal of the reversed divisor, and reverses the result back.
+The remainder uses the same multiplication plan.
+-/
+
+namespace Hex.DensePoly
+
+universe u
+
+attribute [local instance 1000] Lean.Grind.Semiring.ofNat
+
+variable {R : Type u} [DecidableEq R] [Lean.Grind.CommRing R]
+
+/-- Reverse the first `k` coefficients of a represented series back into a
+normalized polynomial.  The explicit bound is independent of any trailing
+zeros removed by polynomial normalization. -/
+def reversePrefix {n : Nat} (k : Nat) (a : TSeries R n) : DensePoly R :=
+  ofList ((List.range k).map fun i => a.coeff (k - 1 - i))
+
+/-- Coefficient law for fixed-bound reverse conversion. -/
+theorem coeff_reversePrefix {n : Nat} (k : Nat) (a : TSeries R n) (i : Nat) :
+    (reversePrefix k a).coeff i =
+      if i < k then a.coeff (k - 1 - i) else 0 := by
+  unfold reversePrefix
+  rw [coeff_ofList]
+  by_cases hi : i < k
+  · simp [List.getD, hi]
+  · rw [List.getD_eq_getElem?_getD]
+    simp [hi]
+    rfl
+
+omit [DecidableEq R] in
+private theorem foldRangeExtend (f : Nat → R) (m n : Nat) (hm : m ≤ n)
+    (hz : ∀ i, m ≤ i → i < n → f i = 0) :
+    (List.range n).foldl (fun acc i => acc + f i) 0 =
+      (List.range m).foldl (fun acc i => acc + f i) 0 := by
+  have hn : m + (n - m) = n := by omega
+  rw [← hn, List.range_add, List.foldl_append, List.foldl_map]
+  apply List.foldl_add_eq_self
+  intro j hj
+  apply hz (m + j) (by omega)
+  have := List.mem_range.mp hj
+  omega
+
+private theorem size_le_of_coeff_zero_above {p : DensePoly R} {N : Nat}
+    (h : ∀ i, N ≤ i → p.coeff i = 0) : p.size ≤ N := by
+  by_cases hle : p.size ≤ N
+  · exact hle
+  · have hpos : 0 < p.size := by omega
+    have hzero := h (p.size - 1) (by omega)
+    exact False.elim (coeff_last_ne_zero_of_pos_size p hpos hzero)
+
+/-- Reversing a `k`-coefficient series prefix converts its product's high
+coefficients into the corresponding low series convolution. -/
+theorem coeff_reversePrefix_mul {n : Nat} (a : TSeries R n)
+    (q : DensePoly R) (k t : Nat) (hk : k ≤ n) (ht : t < k)
+    (hq : 0 < q.size) :
+    (reversePrefix k a * q).coeff (k + q.size - 2 - t) =
+      (reverseSeries q n * a).coeff t := by
+  let base := k - 1 - t
+  have hbase : base + q.size = k + q.size - 1 - t := by
+    dsimp [base]
+    omega
+  have hdeg : k + q.size - 2 - t + 1 = base + q.size := by
+    dsimp [base]
+    omega
+  rw [DensePoly.coeff_mul, mulCoeffSum_eq_diagonal,
+    diagonal_eq_degree_bound, TSeries.coeff_mul _ _ t (by omega)]
+  unfold TSeries.convCoeff
+  rw [show k + q.size - 2 - t + 1 = base + q.size from hdeg,
+    List.range_add, List.foldl_append, List.foldl_map]
+  have hprefix :
+      (List.range base).foldl
+          (fun acc i => acc + diagonalMulCoeffTerm
+            (reversePrefix k a) q (k + q.size - 2 - t) i) 0 = 0 := by
+    apply List.foldl_add_eq_self
+    intro i hi
+    have hib : i < base := List.mem_range.mp hi
+    unfold diagonalMulCoeffTerm
+    have hid : i ≤ k + q.size - 2 - t := by omega
+    rw [_root_.if_neg (by omega)]
+    have hqi : q.size ≤ k + q.size - 2 - t - i := by
+      dsimp [base] at hib
+      omega
+    rw [coeff_eq_zero_of_size_le q hqi]
+    exact Lean.Grind.Semiring.mul_zero _
+  rw [hprefix]
+  have hleft :
+      (List.range q.size).foldl
+          (fun acc j => acc + diagonalMulCoeffTerm
+            (reversePrefix k a) q (k + q.size - 2 - t) (base + j)) 0 =
+        (List.range q.size).foldl
+          (fun acc j => acc +
+            if j < q.size ∧ j ≤ t then
+              q.coeff (q.size - 1 - j) * a.coeff (t - j)
+            else 0) 0 := by
+    apply List.foldl_add_congr
+    intro j hj
+    have hjs : j < q.size := List.mem_range.mp hj
+    unfold diagonalMulCoeffTerm
+    rw [_root_.if_neg (by
+      dsimp [base]
+      omega), coeff_reversePrefix]
+    by_cases hjt : j ≤ t
+    · rw [_root_.ite_eq_left (by
+        dsimp [base]; omega), _root_.ite_eq_left ⟨hjs, hjt⟩]
+      have hidx₁ : k + q.size - 2 - t - (base + j) = q.size - 1 - j := by
+        dsimp [base]
+        omega
+      have hidx₂ : k - 1 - (base + j) = t - j := by
+        dsimp [base]
+        omega
+      rw [hidx₁, hidx₂]
+      grind
+    · rw [_root_.ite_eq_right (by
+        dsimp [base]; omega), _root_.ite_eq_right (by omega)]
+      exact Lean.Grind.Semiring.zero_mul _
+  rw [hleft]
+  have hright :
+      (List.range (t + 1)).foldl
+          (fun acc j => acc + (reverseSeries q n).coeff j * a.coeff (t - j)) 0 =
+        (List.range (t + 1)).foldl
+          (fun acc j => acc +
+            if j < q.size ∧ j ≤ t then
+              q.coeff (q.size - 1 - j) * a.coeff (t - j)
+            else 0) 0 := by
+    apply List.foldl_add_congr
+    intro j hj
+    have hjt : j ≤ t := by
+      have := List.mem_range.mp hj
+      omega
+    rw [coeff_reverseSeries q n j (by omega)]
+    by_cases hjs : j < q.size
+    · rw [_root_.ite_eq_left hjs, _root_.ite_eq_left ⟨hjs, hjt⟩]
+    · rw [_root_.ite_eq_right hjs, _root_.ite_eq_right (by omega)]
+      exact Lean.Grind.Semiring.zero_mul _
+  rw [hright]
+  let f := fun j =>
+    if j < q.size ∧ j ≤ t then
+      q.coeff (q.size - 1 - j) * a.coeff (t - j) else 0
+  let bound := max q.size (t + 1)
+  have hqext := foldRangeExtend f q.size bound (Nat.le_max_left _ _)
+    (by
+      intro i hi _
+      dsimp [f]
+      rw [_root_.ite_eq_right (by omega)])
+  have htext := foldRangeExtend f (t + 1) bound (Nat.le_max_right _ _)
+    (by
+      intro i hi _
+      dsimp [f]
+      rw [_root_.ite_eq_right (by omega)])
+  dsimp [f] at hqext htext
+  exact hqext.symm.trans htext
+
+/-- Number of quotient coefficients required for division of `p` by `q`.
+It is zero for the zero divisor and for a divisor larger than the dividend. -/
+def quotientLength (p q : DensePoly R) : Nat :=
+  if q.size = 0 || p.size < q.size then 0 else p.size - q.size + 1
+
+/-- A reusable divisor with a reciprocal cached to a fixed capacity. -/
+structure DivPlan (R : Type u) [DecidableEq R] [Lean.Grind.CommRing R] where
+  mul : MulPlan R
+  divisor : DensePoly R
+  capacity : Nat
+  unitInv : R
+  reciprocal : TSeries R capacity
+  divisor_ne : divisor ≠ 0
+  unitInv_spec : 0 < capacity →
+    (reverseSeries divisor capacity).coeff 0 * unitInv = 1
+  reciprocal_spec :
+    reciprocal = TSeries.invOfUnit (reverseSeries divisor capacity)
+      unitInv
+
+/-- Build a cached plan for a nonzero monic divisor. -/
+def DivPlan.ofMonic (mul : MulPlan R) (q : DensePoly R)
+    (hq : Monic q) (hqne : q ≠ 0) (capacity : Nat) : DivPlan R :=
+  let reciprocal := reciprocalWith mul (reverseSeries q capacity) 1
+  { mul
+    divisor := q
+    capacity
+    unitInv := 1
+    reciprocal
+    divisor_ne := hqne
+    unitInv_spec := by
+      intro hcap
+      rw [coeff_reverseSeries_of_lt q capacity 0 hcap (by
+        have : q.size ≠ 0 := by
+          intro hs
+          exact hqne ((size_eq_zero_iff q).mp hs)
+        omega)]
+      have hqpos : 0 < q.size := by
+        apply Nat.pos_of_ne_zero
+        intro hs
+        exact hqne ((size_eq_zero_iff q).mp hs)
+      simp only [Nat.sub_zero]
+      have hlead : q.coeff (q.size - 1) = q.leadingCoeff := by
+        rfl
+      rw [hlead, leadingCoeff_eq_one_of_monic hq]
+      grind
+    reciprocal_spec := by
+      dsimp [reciprocal]
+      exact reciprocalWith_eq mul _ 1 }
+
+/-- Build a cached plan for an arbitrary nonzero divisor over a field. -/
+def DivPlan.ofNonzero {F : Type u} [DecidableEq F] [Lean.Grind.Field F]
+    (mul : MulPlan F) (q : DensePoly F) (hqne : q ≠ 0)
+    (capacity : Nat) : DivPlan F :=
+  let u := q.leadingCoeff⁻¹
+  let reciprocal := reciprocalWith mul (reverseSeries q capacity) u
+  { mul
+    divisor := q
+    capacity
+    unitInv := u
+    reciprocal
+    divisor_ne := hqne
+    unitInv_spec := by
+      intro hcap
+      rw [coeff_reverseSeries_of_lt q capacity 0 hcap (by
+        have : q.size ≠ 0 := by
+          intro hs
+          exact hqne ((size_eq_zero_iff q).mp hs)
+        omega)]
+      simp only [Nat.sub_zero]
+      have hlead : q.coeff (q.size - 1) = q.leadingCoeff := by
+        rfl
+      rw [hlead]
+      have hqpos : 0 < q.size := by
+        apply Nat.pos_of_ne_zero
+        intro hs
+        exact hqne ((size_eq_zero_iff q).mp hs)
+      have hne : q.leadingCoeff ≠ 0 := leadingCoeff_ne_zero_of_pos_size q hqpos
+      dsimp [u]
+      exact Lean.Grind.Field.mul_inv_cancel hne
+    reciprocal_spec := by
+      dsimp [reciprocal]
+      exact reciprocalWith_eq mul _ u }
+
+/-- Quotient obtained from a cached reversed reciprocal. -/
+def DivPlan.quotient (plan : DivPlan R) (p : DensePoly R)
+    (_hcap : quotientLength p plan.divisor ≤ plan.capacity) : DensePoly R :=
+  let k := quotientLength p plan.divisor
+  let top := reverseSeries p plan.capacity
+  let qrev := seriesMulUpTo plan.mul k top plan.reciprocal
+  reversePrefix k qrev
+
+/-- The reciprocal quotient cancels every coefficient in the high window of
+the dividend.  The index `t` counts down from the leading coefficient. -/
+theorem DivPlan.coeff_quotient_mul_high (plan : DivPlan R) (p : DensePoly R)
+    (hcap : quotientLength p plan.divisor ≤ plan.capacity) (t : Nat)
+    (ht : t < quotientLength p plan.divisor) :
+    (plan.quotient p hcap * plan.divisor).coeff (p.size - 1 - t) =
+      p.coeff (p.size - 1 - t) := by
+  let k := quotientLength p plan.divisor
+  let top := reverseSeries p plan.capacity
+  let qrev := seriesMulUpTo plan.mul k top plan.reciprocal
+  have hdpos : 0 < plan.divisor.size := by
+    apply Nat.pos_of_ne_zero
+    intro hs
+    exact plan.divisor_ne ((size_eq_zero_iff plan.divisor).mp hs)
+  have hdle : plan.divisor.size ≤ p.size := by
+    by_cases hle : plan.divisor.size ≤ p.size
+    · exact hle
+    · have hlt : p.size < plan.divisor.size := Nat.lt_of_not_ge hle
+      have hkzero : k = 0 := by simp [k, quotientLength, hlt]
+      omega
+  have hkform : k = p.size - plan.divisor.size + 1 := by
+    simp [k, quotientLength, plan.divisor_ne, Nat.not_lt.mpr hdle]
+  have hkpos : 0 < k := by omega
+  have hklep : k ≤ p.size := by
+    rw [hkform]
+    omega
+  have hkcap : k ≤ plan.capacity := hcap
+  have hunit : (reverseSeries plan.divisor plan.capacity).coeff 0 *
+      plan.unitInv = 1 := plan.unitInv_spec (by omega)
+  have hinv : reverseSeries plan.divisor plan.capacity * plan.reciprocal = 1 := by
+    rw [plan.reciprocal_spec]
+    exact TSeries.invOfUnit_mul _ _ hunit
+  have hqrev : TSeries.Agree k qrev (top * plan.reciprocal) := by
+    dsimp [qrev]
+    rw [seriesMulUpTo_eq]
+    exact TSeries.Agree.mulUpTo k top plan.reciprocal
+  have hprod : TSeries.Agree k
+      (reverseSeries plan.divisor plan.capacity * qrev) top := by
+    have hmul := TSeries.Agree.mul
+      (TSeries.Agree.refl k (reverseSeries plan.divisor plan.capacity)) hqrev
+    intro i hi hik
+    have h := hmul i hi hik
+    rw [← TSeries.mul_assoc,
+      TSeries.mul_comm (reverseSeries plan.divisor plan.capacity) top,
+      TSeries.mul_assoc, hinv, TSeries.mul_one] at h
+    exact h
+  have hidx : k + plan.divisor.size - 2 - t = p.size - 1 - t := by
+    rw [hkform]
+    omega
+  have hbridge :
+      (reversePrefix k qrev * plan.divisor).coeff (p.size - 1 - t) =
+        (reverseSeries plan.divisor plan.capacity * qrev).coeff t := by
+    rw [← hidx]
+    exact coeff_reversePrefix_mul qrev plan.divisor k t hkcap ht hdpos
+  change (reversePrefix k qrev * plan.divisor).coeff (p.size - 1 - t) = _
+  rw [hbridge, hprod t (by omega) ht]
+  dsimp [top]
+  rw [coeff_reverseSeries_of_lt p plan.capacity t (by omega) (by omega)]
+
+/-- The normalized reciprocal quotient has no more than the requested number
+of coefficients. -/
+theorem DivPlan.size_quotient_le (plan : DivPlan R) (p : DensePoly R)
+    (hcap : quotientLength p plan.divisor ≤ plan.capacity) :
+    (plan.quotient p hcap).size ≤ quotientLength p plan.divisor := by
+  unfold DivPlan.quotient reversePrefix
+  exact Nat.le_trans (size_ofList_le _) (by simp)
+
+/-- Subtracting the reciprocal quotient product leaves fewer coefficients than
+the divisor.  This is the executable division algorithm's termination law. -/
+theorem DivPlan.remainder_size_le (plan : DivPlan R) (p : DensePoly R)
+    (hcap : quotientLength p plan.divisor ≤ plan.capacity) :
+    (p - mulWith plan.mul (plan.quotient p hcap) plan.divisor).size ≤
+      plan.divisor.size - 1 := by
+  apply size_le_of_coeff_zero_above
+  intro i hi
+  rw [coeff_sub_ring, mulWith_eq]
+  by_cases hip : i < p.size
+  · have hdle : plan.divisor.size ≤ p.size := by omega
+    have hkform : quotientLength p plan.divisor =
+        p.size - plan.divisor.size + 1 := by
+      simp [quotientLength, plan.divisor_ne, Nat.not_lt.mpr hdle]
+    let t := p.size - 1 - i
+    have ht : t < quotientLength p plan.divisor := by
+      dsimp [t]
+      rw [hkform]
+      omega
+    have hcancel := plan.coeff_quotient_mul_high p hcap t ht
+    have hidx : p.size - 1 - t = i := by
+      dsimp [t]
+      omega
+    rw [hidx] at hcancel
+    rw [hcancel]
+    grind
+  · rw [coeff_eq_zero_of_size_le p (Nat.le_of_not_gt hip)]
+    have hqsize := plan.size_quotient_le p hcap
+    have hprod := size_mul_le (plan.quotient p hcap) plan.divisor
+    by_cases hdp : plan.divisor.size ≤ p.size
+    · have hkform : quotientLength p plan.divisor =
+          p.size - plan.divisor.size + 1 := by
+        simp [quotientLength, plan.divisor_ne, Nat.not_lt.mpr hdp]
+      rw [coeff_eq_zero_of_size_le _ (by rw [hkform] at hqsize; omega)]
+      grind
+    · have hdlt : p.size < plan.divisor.size := Nat.lt_of_not_ge hdp
+      have hkzero : quotientLength p plan.divisor = 0 := by
+        simp [quotientLength, hdlt]
+      rw [coeff_eq_zero_of_size_le _ (by rw [hkzero] at hqsize; omega)]
+      grind
+
+/-- The cached division remainder has degree strictly below every
+positive-degree divisor. -/
+theorem DivPlan.remainder_degree_lt (plan : DivPlan R) (p : DensePoly R)
+    (hcap : quotientLength p plan.divisor ≤ plan.capacity)
+    (hdeg : 0 < plan.divisor.degree?.getD 0) :
+    (p - mulWith plan.mul (plan.quotient p hcap) plan.divisor).degree?.getD 0 <
+      plan.divisor.degree?.getD 0 := by
+  have hdpos : 0 < plan.divisor.size := by
+    rcases Nat.eq_zero_or_pos plan.divisor.size with hz | hz
+    · rw [(degree?_eq_none_iff plan.divisor).mpr hz, Option.getD_none] at hdeg
+      omega
+    · exact hz
+  have hddegree : plan.divisor.degree?.getD 0 = plan.divisor.size - 1 := by
+    rw [degree?_eq_some_of_pos_size plan.divisor hdpos, Option.getD_some]
+  let r := p - mulWith plan.mul (plan.quotient p hcap) plan.divisor
+  have hrsize : r.size ≤ plan.divisor.size - 1 := plan.remainder_size_le p hcap
+  rcases Nat.eq_zero_or_pos r.size with hz | hz
+  · rw [(degree?_eq_none_iff r).mpr hz, Option.getD_none, hddegree]
+    omega
+  · rw [degree?_eq_some_of_pos_size r hz, Option.getD_some, hddegree]
+    omega
+
+/-- Divide using a cached reciprocal.  The proof ensures the cached precision
+covers the requested quotient; it is erased from executable code. -/
+def DivPlan.divMod (plan : DivPlan R) (p : DensePoly R)
+    (hcap : quotientLength p plan.divisor ≤ plan.capacity) :
+    DensePoly R × DensePoly R :=
+  let q := plan.quotient p hcap
+  (q, p - mulWith plan.mul q plan.divisor)
+
+/-- Cached reciprocal division agrees with the established long-division
+result whenever coefficient division satisfies the laws needed by that
+result.  The size-one divisor case is handled as an exact product rather than
+being hidden by the default-zero degree convention. -/
+theorem DivPlan.divMod_eq_divMod [Div R] (plan : DivPlan R) (p : DensePoly R)
+    (hcap : quotientLength p plan.divisor ≤ plan.capacity)
+    (hcancel : ∀ a : R,
+      a - (a / plan.divisor.leadingCoeff) * plan.divisor.leadingCoeff = 0)
+    (hexact : ∀ a : R,
+      (a * plan.divisor.leadingCoeff) / plan.divisor.leadingCoeff = a)
+    (h_top_ne : ∀ a : R, a ≠ 0 →
+      a * plan.divisor.leadingCoeff ≠ 0) :
+    plan.divMod p hcap = _root_.Hex.DensePoly.divMod p plan.divisor := by
+  let q := plan.quotient p hcap
+  let r := p - mulWith plan.mul q plan.divisor
+  have hrec : q * plan.divisor + r = p := by
+    apply ext_coeff
+    intro i
+    rw [coeff_add_semiring, coeff_sub_ring, mulWith_eq]
+    grind
+  have hdpos : 0 < plan.divisor.size := by
+    apply Nat.pos_of_ne_zero
+    intro hs
+    exact plan.divisor_ne ((size_eq_zero_iff plan.divisor).mp hs)
+  have heq : _root_.Hex.DensePoly.divMod p plan.divisor = (q, r) := by
+    by_cases hdeg : 0 < plan.divisor.degree?.getD 0
+    · exact divMod_eq_of_reconstruction p plan.divisor q r hdeg hcancel
+        hexact h_top_ne hrec (plan.remainder_degree_lt p hcap hdeg)
+    · have hddegree : plan.divisor.degree?.getD 0 =
+          plan.divisor.size - 1 := by
+        rw [degree?_eq_some_of_pos_size plan.divisor hdpos, Option.getD_some]
+      have hdsize : plan.divisor.size = 1 := by omega
+      have hrsize : r.size ≤ plan.divisor.size - 1 :=
+        plan.remainder_size_le p hcap
+      have hrzero : r = 0 := by
+        apply (size_eq_zero_iff r).mp
+        rw [hdsize] at hrsize
+        omega
+      have hmul : q * plan.divisor = p := by
+        rw [hrzero, add_zero_poly] at hrec
+        exact hrec
+      simpa [hrzero] using
+        (divMod_eq_of_polynomial_mul p plan.divisor q plan.divisor_ne
+          hexact h_top_ne hmul)
+  change (q, r) = _root_.Hex.DensePoly.divMod p plan.divisor
+  exact heq.symm
+
+/-- Remainder-only cached division. -/
+def DivPlan.mod (plan : DivPlan R) (p : DensePoly R)
+    (hcap : quotientLength p plan.divisor ≤ plan.capacity) : DensePoly R :=
+  (plan.divMod p hcap).2
+
+/-- One-shot reciprocal division by a monic polynomial. -/
+def divModMonicWith (mul : MulPlan R) (p q : DensePoly R) (hq : Monic q) :
+    DensePoly R × DensePoly R :=
+  if hqne : q = 0 then
+    (0, p)
+  else
+    let k := quotientLength p q
+    let plan := DivPlan.ofMonic mul q hq hqne k
+    plan.divMod p (by
+      dsimp [plan, DivPlan.ofMonic]
+      exact Nat.le_refl k)
+
+/-- One-shot monic reciprocal division agrees with the existing specialized
+monic long-division operation. -/
+theorem divModMonicWith_eq (mul : MulPlan R) (p q : DensePoly R)
+    (hq : Monic q) :
+    divModMonicWith mul p q hq = divModMonic p q hq := by
+  unfold divModMonicWith
+  split
+  · rename_i hqzero
+    subst q
+    have hzero_one : (0 : R) = 1 := by
+      simpa using (leadingCoeff_eq_one_of_monic hq)
+    have hpzero : p = 0 := by
+      apply ext_coeff
+      intro i
+      rw [coeff_zero]
+      calc
+        p.coeff i = p.coeff i * 1 := by grind
+        _ = p.coeff i * 0 := by rw [hzero_one]
+        _ = 0 := by grind
+    subst p
+    letI : Div R := ⟨fun a _ => a⟩
+    have hnot_lt : ¬(0 : DensePoly R).degree?.getD 0 <
+        (0 : DensePoly R).degree?.getD 0 := by omega
+    rw [divModMonic_eq_divMod_of_monic_of_scale 0 0 hq hnot_lt
+      (fun _ => rfl)]
+    exact (divMod_eq_zero_self_of_size_zero 0 0 size_zero).symm
+  · rename_i hqne
+    let k := quotientLength p q
+    let plan := DivPlan.ofMonic mul q hq hqne k
+    letI : Div R := ⟨fun a _ => a⟩
+    have hlead : q.leadingCoeff = 1 := leadingCoeff_eq_one_of_monic hq
+    have hplan : plan.divMod p (by
+        dsimp [plan, DivPlan.ofMonic]
+        exact Nat.le_refl k) = _root_.Hex.DensePoly.divMod p q := by
+      apply plan.divMod_eq_divMod
+      · intro a
+        change a - a * q.leadingCoeff = 0
+        rw [hlead]
+        grind
+      · intro a
+        change a * q.leadingCoeff = a
+        rw [hlead]
+        grind
+      · intro a ha
+        change a * q.leadingCoeff ≠ 0
+        rw [hlead]
+        simpa only [Lean.Grind.Semiring.mul_one] using ha
+    by_cases hlt : p.degree?.getD 0 < q.degree?.getD 0
+    · have hfast : _root_.Hex.DensePoly.divMod p q = (0, p) :=
+        divMod_eq_zero_self_of_degree_lt p q hlt
+      have hmonic : divModMonic p q hq = (0, p) := by
+        unfold divModMonic
+        exact divModArray_eq_zero_self_of_degree_lt p q id hlt
+      rw [hplan, hfast, hmonic]
+    · have hmonic := divModMonic_eq_divMod_of_monic_of_scale p q hq hlt
+        (fun _ => rfl)
+      exact hplan.trans hmonic.symm
+
+/-- One-shot reciprocal division over a field. -/
+def divModWith {F : Type u} [DecidableEq F] [Lean.Grind.Field F]
+    (mul : MulPlan F) (p q : DensePoly F) : DensePoly F × DensePoly F :=
+  if hqne : q = 0 then
+    (0, p)
+  else
+    let k := quotientLength p q
+    let plan := DivPlan.ofNonzero mul q hqne k
+    plan.divMod p (by
+      dsimp [plan, DivPlan.ofNonzero]
+      exact Nat.le_refl k)
+
+/-- One-shot field reciprocal division is extensionally the existing verified
+long-division operation. -/
+theorem divModWith_eq {F : Type u} [DecidableEq F] [Lean.Grind.Field F]
+    (mul : MulPlan F) (p q : DensePoly F) :
+    divModWith mul p q = _root_.Hex.DensePoly.divMod p q := by
+  unfold divModWith
+  split
+  · rename_i hq
+    have hqsize : q.size = 0 := (size_eq_zero_iff q).mpr hq
+    exact (divMod_eq_zero_self_of_size_zero p q hqsize).symm
+  · rename_i hqne
+    let k := quotientLength p q
+    let plan := DivPlan.ofNonzero mul q hqne k
+    have hlead : q.leadingCoeff ≠ 0 := by
+      apply leadingCoeff_ne_zero_of_pos_size q
+      apply Nat.pos_of_ne_zero
+      intro hs
+      exact hqne ((size_eq_zero_iff q).mp hs)
+    apply plan.divMod_eq_divMod
+    · intro a
+      change a - (a / q.leadingCoeff) * q.leadingCoeff = 0
+      rw [Lean.Grind.Field.div_eq_mul_inv, Lean.Grind.Semiring.mul_assoc,
+        Lean.Grind.Field.inv_mul_cancel hlead, Lean.Grind.Semiring.mul_one]
+      grind
+    · intro a
+      change (a * q.leadingCoeff) / q.leadingCoeff = a
+      rw [Lean.Grind.Field.div_eq_mul_inv, Lean.Grind.Semiring.mul_assoc,
+        Lean.Grind.Field.mul_inv_cancel hlead, Lean.Grind.Semiring.mul_one]
+    · intro a ha
+      change a * q.leadingCoeff ≠ 0
+      intro hz
+      have hz' := congrArg (fun y : F => y * q.leadingCoeff⁻¹) hz
+      rw [Lean.Grind.Semiring.zero_mul, Lean.Grind.Semiring.mul_assoc,
+        Lean.Grind.Field.mul_inv_cancel hlead, Lean.Grind.Semiring.mul_one] at hz'
+      exact ha hz'
+
+end Hex.DensePoly
diff --git a/HexPolyFast/Karatsuba.lean b/HexPolyFast/Karatsuba.lean
new file mode 100644
index 000000000..adf24a172
--- /dev/null
+++ b/HexPolyFast/Karatsuba.lean
@@ -0,0 +1,500 @@
+/-
+Copyright (c) 2026 Lean FRO, LLC. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Kim Morrison
+-/
+
+module
+
+public import HexPolyFast.Plan
+
+public section
+
+/-!
+Splitting and algebraic assembly lemmas for Karatsuba multiplication.
+-/
+
+namespace Hex.DensePoly
+
+universe u
+
+attribute [local instance 1000] Lean.Grind.Semiring.ofNat
+
+variable {R : Type u} [DecidableEq R] [Lean.Grind.CommRing R]
+
+/-- Coefficients strictly below `k`. -/
+@[expose]
+def low (k : Nat) (p : DensePoly R) : DensePoly R :=
+  ofList ((List.range k).map p.coeff)
+
+/-- Coefficients at and above `k`, shifted down by `k`. -/
+@[expose]
+def high (k : Nat) (p : DensePoly R) : DensePoly R :=
+  ofList ((List.range (p.size - k)).map fun i => p.coeff (k + i))
+
+/-- Coefficient law for the low half. -/
+theorem coeff_low (k : Nat) (p : DensePoly R) (i : Nat) :
+    (low k p).coeff i = if i < k then p.coeff i else 0 := by
+  unfold low
+  rw [coeff_ofList]
+  by_cases hi : i < k
+  · simp [List.getD, hi]
+  · rw [List.getD_eq_getElem?_getD]
+    simp [hi]
+    rfl
+
+/-- Coefficient law for the shifted-down high half. -/
+theorem coeff_high (k : Nat) (p : DensePoly R) (i : Nat) :
+    (high k p).coeff i =
+      if i < p.size - k then p.coeff (k + i) else 0 := by
+  unfold high
+  rw [coeff_ofList]
+  by_cases hi : i < p.size - k
+  · simp [List.getD, hi]
+  · rw [List.getD_eq_getElem?_getD]
+    simp [hi]
+    rfl
+
+/-- Splitting at `k` and shifting the high half back reconstructs the input. -/
+theorem low_add_shift_high (k : Nat) (p : DensePoly R) :
+    low k p + shift k (high k p) = p := by
+  apply ext_coeff
+  intro i
+  have hz : (0 : R) + 0 = 0 := by grind
+  rw [coeff_add _ _ _ hz, coeff_low, coeff_shift]
+  by_cases hik : i < k
+  · rw [_root_.ite_eq_left hik, _root_.ite_eq_left hik]
+    change p.coeff i + (0 : R) = p.coeff i
+    exact Lean.Grind.Semiring.add_zero _
+  · have hki : k ≤ i := Nat.le_of_not_gt hik
+    rw [_root_.ite_eq_right hik, _root_.ite_eq_right hik]
+    by_cases hip : i < p.size
+    · have hh : i - k < p.size - k := by omega
+      rw [coeff_high, _root_.ite_eq_left hh]
+      have hidx : k + (i - k) = i := by omega
+      rw [hidx]
+      grind
+    · have hsize : p.size ≤ i := Nat.le_of_not_gt hip
+      have hh : ¬ i - k < p.size - k := by omega
+      rw [coeff_high, _root_.ite_eq_right hh,
+        coeff_eq_zero_of_size_le p hsize]
+      exact hz
+
+/-- Nested shifts add their offsets. -/
+theorem shift_shift (k l : Nat) (p : DensePoly R) :
+    shift k (shift l p) = shift (k + l) p := by
+  apply ext_coeff
+  intro i
+  rw [coeff_shift, coeff_shift, coeff_shift]
+  by_cases hi : i < k
+  · have hil : i < k + l := by omega
+    simp [hi, hil]
+  · by_cases hil : i - k < l
+    · have hisum : i < k + l := by omega
+      simp [hi, hil, hisum]
+    · have hisum : ¬ i < k + l := by omega
+      have hidx : i - k - l = i - (k + l) := by omega
+      simp [hi, hil, hisum, hidx]
+
+/-- A shifted left factor shifts the product. -/
+theorem shift_mul (k : Nat) (p q : DensePoly R) :
+    shift k p * q = shift k (p * q) := by
+  calc
+    shift k p * q = (monomial k 1 * p) * q := by
+      rw [monomial_one_mul_poly_eq_shift]
+    _ = monomial k 1 * (p * q) := mul_assoc_poly _ _ _
+    _ = shift k (p * q) := monomial_one_mul_poly_eq_shift _ _
+
+/-- A shifted right factor shifts the product. -/
+theorem mul_shift (k : Nat) (p q : DensePoly R) :
+    p * shift k q = shift k (p * q) := by
+  calc
+    p * shift k q = shift k q * p := mul_comm_poly _ _
+    _ = shift k (q * p) := shift_mul _ _ _
+    _ = shift k (p * q) := congrArg (shift k) (mul_comm_poly q p)
+
+/-- Multiplying two shifted polynomials adds their offsets. -/
+theorem shift_mul_shift (k l : Nat) (p q : DensePoly R) :
+    shift k p * shift l q = shift (k + l) (p * q) := by
+  rw [shift_mul, mul_shift, shift_shift]
+
+/-- The three-product Karatsuba assembly identity. -/
+theorem karatsuba_combine (k : Nat) (a₀ a₁ b₀ b₁ : DensePoly R) :
+    let z₀ := a₀ * b₀
+    let z₂ := a₁ * b₁
+    let z₁ := (a₀ + a₁) * (b₀ + b₁) - z₀ - z₂
+    z₀ + shift k z₁ + shift (2 * k) z₂ =
+      (a₀ + shift k a₁) * (b₀ + shift k b₁) := by
+  dsimp
+  rw [mul_add_right_poly (a₀ + a₁) b₀ b₁,
+    mul_add_left_poly a₀ a₁ b₀,
+    mul_add_left_poly a₀ a₁ b₁,
+    mul_add_right_poly (a₀ + shift k a₁) b₀ (shift k b₁),
+    mul_add_left_poly a₀ (shift k a₁) b₀,
+    mul_add_left_poly a₀ (shift k a₁) (shift k b₁),
+    shift_mul, mul_shift, shift_mul_shift]
+  have hkk : k + k = 2 * k := by omega
+  rw [hkk]
+  apply ext_coeff
+  intro i
+  have hzadd : (0 : R) + 0 = 0 := by grind
+  have hzsub : (0 : R) - 0 = 0 := by grind
+  simp only [coeff_add _ _ _ hzadd, coeff_sub _ _ _ hzsub, coeff_shift]
+  have hzero : (Zero.zero : R) = 0 := rfl
+  simp only [hzero]
+  split <;> split <;> grind
+
+/-- Fuelled three-product Karatsuba recursion.  Fuel is separate from the
+cutoff so cutoff zero remains total; public callers provide at least the
+larger operand size. -/
+@[expose]
+def karatsubaAux (cutoff : Nat) : Nat → DensePoly R → DensePoly R → DensePoly R
+  | 0, a, b => mulImpl a b
+  | fuel + 1, a, b =>
+      if a.size ≤ max 1 cutoff || b.size ≤ max 1 cutoff then
+        mulImpl a b
+      else
+        let k := (max a.size b.size + 1) / 2
+        let a₀ := low k a
+        let a₁ := high k a
+        let b₀ := low k b
+        let b₁ := high k b
+        let z₀ := karatsubaAux cutoff fuel a₀ b₀
+        let z₂ := karatsubaAux cutoff fuel a₁ b₁
+        let z₁ := karatsubaAux cutoff fuel (a₀ + a₁) (b₀ + b₁) - z₀ - z₂
+        z₀ + shift k z₁ + shift (2 * k) z₂
+
+/-- Specialized three-square Karatsuba recursion. -/
+@[expose]
+def karatsubaSquareAux (cutoff : Nat) : Nat → DensePoly R → DensePoly R
+  | 0, a => mulImpl a a
+  | fuel + 1, a =>
+      if a.size ≤ max 1 cutoff then
+        mulImpl a a
+      else
+        let k := (a.size + 1) / 2
+        let a₀ := low k a
+        let a₁ := high k a
+        let z₀ := karatsubaSquareAux cutoff fuel a₀
+        let z₂ := karatsubaSquareAux cutoff fuel a₁
+        let z₁ := karatsubaSquareAux cutoff fuel (a₀ + a₁) - z₀ - z₂
+        z₀ + shift k z₁ + shift (2 * k) z₂
+
+/-- Every fuelled Karatsuba product agrees with schoolbook multiplication. -/
+theorem karatsubaAux_eq (cutoff fuel : Nat) (a b : DensePoly R) :
+    karatsubaAux cutoff fuel a b = a * b := by
+  induction fuel generalizing a b with
+  | zero =>
+      exact (mul_eq_mulImpl a b).symm
+  | succ fuel ih =>
+      rw [karatsubaAux]
+      split
+      · exact (mul_eq_mulImpl a b).symm
+      · dsimp only
+        rw [ih, ih, ih]
+        rw [karatsuba_combine]
+        rw [low_add_shift_high, low_add_shift_high]
+
+/-- Every fuelled Karatsuba square agrees with ordinary squaring. -/
+theorem karatsubaSquareAux_eq (cutoff fuel : Nat) (a : DensePoly R) :
+    karatsubaSquareAux cutoff fuel a = a * a := by
+  induction fuel generalizing a with
+  | zero =>
+      exact (mul_eq_mulImpl a a).symm
+  | succ fuel ih =>
+      rw [karatsubaSquareAux]
+      split
+      · exact (mul_eq_mulImpl a a).symm
+      · dsimp only
+        rw [ih, ih, ih]
+        rw [karatsuba_combine]
+        rw [low_add_shift_high]
+
+/-- Balanced Karatsuba multiplication with an explicit schoolbook cutoff. -/
+def mulKaratsubaBalanced (cutoff : Nat) (a b : DensePoly R) : DensePoly R :=
+  karatsubaAux cutoff (max a.size b.size) a b
+
+/-- Specialized Karatsuba squaring with an explicit schoolbook cutoff. -/
+def squareKaratsuba (cutoff : Nat) (a : DensePoly R) : DensePoly R :=
+  karatsubaSquareAux cutoff a.size a
+
+/-- Balanced Karatsuba multiplication agrees exactly with `DensePoly.mul`. -/
+theorem mulKaratsubaBalanced_eq (cutoff : Nat) (a b : DensePoly R) :
+    mulKaratsubaBalanced cutoff a b = a * b :=
+  karatsubaAux_eq cutoff _ a b
+
+/-- Specialized Karatsuba squaring agrees exactly with `DensePoly.mul`. -/
+theorem squareKaratsuba_eq (cutoff : Nat) (a : DensePoly R) :
+    squareKaratsuba cutoff a = a * a :=
+  karatsubaSquareAux_eq cutoff _ a
+
+/-- Multiply a long operand by blocks of `blockSize` coefficients.  Each block
+uses balanced Karatsuba; the shifted block products are accumulated without
+padding the short operand to the long size. -/
+@[expose]
+def karatsubaBlocks (cutoff blockSize : Nat) :
+    Nat → DensePoly R → DensePoly R → DensePoly R
+  | 0, long, short => mulKaratsubaBalanced cutoff long short
+  | fuel + 1, long, short =>
+      if long.size = 0 then 0
+      else
+        mulKaratsubaBalanced cutoff (low blockSize long) short +
+          shift blockSize
+            (karatsubaBlocks cutoff blockSize fuel (high blockSize long) short)
+
+/-- Blocked unbalanced multiplication agrees with the ordinary product for
+every fuel value. -/
+theorem karatsubaBlocks_eq (cutoff blockSize fuel : Nat)
+    (long short : DensePoly R) :
+    karatsubaBlocks cutoff blockSize fuel long short = long * short := by
+  induction fuel generalizing long with
+  | zero => exact mulKaratsubaBalanced_eq cutoff long short
+  | succ fuel ih =>
+      rw [karatsubaBlocks]
+      split <;> rename_i hzero
+      · have hlong : long = 0 := (size_eq_zero_iff long).mp hzero
+        subst long
+        exact (zero_mul short).symm
+      · rw [mulKaratsubaBalanced_eq, ih, ← shift_mul,
+          ← mul_add_left_poly, low_add_shift_high]
+
+/-- Full Karatsuba multiplication.  Strongly skewed operands are processed in
+blocks near the shorter size rather than padded to the longer size. -/
+def mulKaratsuba (cutoff : Nat) (a b : DensePoly R) : DensePoly R :=
+  if a.size = 0 || b.size = 0 then
+    mulImpl a b
+  else if 2 * b.size < a.size then
+    karatsubaBlocks cutoff b.size a.size a b
+  else if 2 * a.size < b.size then
+    karatsubaBlocks cutoff a.size b.size b a
+  else
+    mulKaratsubaBalanced cutoff a b
+
+/-- Full balanced-or-blocked Karatsuba multiplication agrees exactly with
+`DensePoly.mul`. -/
+theorem mulKaratsuba_eq (cutoff : Nat) (a b : DensePoly R) :
+    mulKaratsuba cutoff a b = a * b := by
+  unfold mulKaratsuba
+  split
+  · exact (mul_eq_mulImpl a b).symm
+  · split
+    · exact karatsubaBlocks_eq cutoff b.size a.size a b
+    · split
+      · rw [karatsubaBlocks_eq, mul_comm_poly]
+      · exact mulKaratsubaBalanced_eq cutoff a b
+
+/-- Extract `len` coefficients beginning at `lo` from an already-computed
+polynomial, shifting them down to degree zero. -/
+@[expose]
+def coeffSlice (lo len : Nat) (p : DensePoly R) : DensePoly R :=
+  ofList ((List.range len).map fun i => p.coeff (lo + i))
+
+/-- Coefficient law for extracting a polynomial interval. -/
+theorem coeff_coeffSlice (lo len : Nat) (p : DensePoly R) (i : Nat) :
+    (coeffSlice lo len p).coeff i =
+      if i < len then p.coeff (lo + i) else 0 := by
+  unfold coeffSlice
+  rw [coeff_ofList]
+  by_cases hi : i < len
+  · simp [List.getD, hi]
+  · rw [List.getD_eq_getElem?_getD]
+    simp [hi]
+    rfl
+
+/-- Truncating both inputs above `n` preserves every product coefficient
+strictly below `n`. -/
+theorem coeff_low_mul_low (n : Nat) (a b : DensePoly R) (i : Nat)
+    (hi : i < n) :
+    (low n a * low n b).coeff i = (a * b).coeff i := by
+  have hdecomp :
+      (low n a + shift n (high n a)) *
+          (low n b + shift n (high n b)) = a * b := by
+    rw [low_add_shift_high, low_add_shift_high]
+  rw [← hdecomp]
+  rw [mul_add_right_poly, mul_add_left_poly, mul_add_left_poly,
+    shift_mul, mul_shift, shift_mul_shift]
+  have hz : (0 : R) + 0 = 0 := by grind
+  simp only [coeff_add _ _ _ hz, coeff_shift]
+  have h2i : i < n + n := by omega
+  have hzero : (Zero.zero : R) = 0 := rfl
+  simp only [hzero]
+  simp [hi, h2i]
+  grind
+
+/-- Fuelled interval-pruned Karatsuba recursion.  At a split, each of the
+three mathematical subproducts is computed once over the smallest bounding
+interval containing all of its shifted contributions to `[lo, lo + len)`.
+Branches whose shifted output starts beyond the interval receive length zero. -/
+@[expose]
+def karatsubaSliceAux (cutoff : Nat) :
+    Nat → Nat → Nat → DensePoly R → DensePoly R → DensePoly R
+  | 0, lo, len, a, b => schoolbookSlice lo len a b
+  | fuel + 1, lo, len, a, b =>
+      if a.size ≤ max 1 cutoff || b.size ≤ max 1 cutoff then
+        schoolbookSlice lo len a b
+      else
+        let k := (max a.size b.size + 1) / 2
+        let a₀ := low k a
+        let a₁ := high k a
+        let b₀ := low k b
+        let b₁ := high k b
+        let hi := lo + len
+        let base₀ := lo - k
+        let base₁ := lo - k
+        let base₂ := lo - 2 * k
+        let z₀ := karatsubaSliceAux cutoff fuel base₀ (hi - base₀) a₀ b₀
+        let z₁ := karatsubaSliceAux cutoff fuel base₁
+          ((hi - k) - base₁) (a₀ + a₁) (b₀ + b₁)
+        let z₂ := karatsubaSliceAux cutoff fuel base₂
+          ((hi - k) - base₂) a₁ b₁
+        ofList ((List.range len).map fun i =>
+          let d := lo + i
+          z₀.coeff (d - base₀) +
+            (if k ≤ d then
+              z₁.coeff (d - k - base₁) -
+                z₀.coeff (d - k - base₀) -
+                z₂.coeff (d - k - base₂)
+            else 0) +
+            if 2 * k ≤ d then z₂.coeff (d - 2 * k - base₂) else 0)
+
+/-- Every fuelled interval recursion returns exactly the requested product
+coefficients. -/
+theorem coeff_karatsubaSliceAux (cutoff fuel lo len : Nat)
+    (a b : DensePoly R) (i : Nat) :
+    (karatsubaSliceAux cutoff fuel lo len a b).coeff i =
+      if i < len then (a * b).coeff (lo + i) else 0 := by
+  induction fuel generalizing lo len a b i with
+  | zero => exact coeff_schoolbookSlice lo len a b i
+  | succ fuel ih =>
+      rw [karatsubaSliceAux]
+      split
+      · exact coeff_schoolbookSlice lo len a b i
+      · dsimp only
+        rw [coeff_ofList]
+        by_cases hil : i < len
+        · rw [List.getD_eq_getElem?_getD]
+          simp only [List.getElem?_map, List.getElem?_range, hil,
+            Option.map_some, Option.getD_some]
+          let k := (max a.size b.size + 1) / 2
+          let a₀ := low k a
+          let a₁ := high k a
+          let b₀ := low k b
+          let b₁ := high k b
+          let hi := lo + len
+          let base₀ := lo - k
+          let base₁ := lo - k
+          let base₂ := lo - 2 * k
+          have hprod :
+              a₀ * b₀ + shift k
+                  ((a₀ + a₁) * (b₀ + b₁) - a₀ * b₀ - a₁ * b₁) +
+                    shift (2 * k) (a₁ * b₁) = a * b := by
+            rw [karatsuba_combine]
+            dsimp [a₀, a₁, b₀, b₁]
+            rw [low_add_shift_high, low_add_shift_high]
+          simp only [ih]
+          have hc := congrArg (fun p : DensePoly R => p.coeff (lo + i)) hprod
+          have hz : (0 : R) + 0 = 0 := by grind
+          have hzs : (0 : R) - 0 = 0 := by grind
+          simp only [coeff_add _ _ _ hz, coeff_sub _ _ _ hzs,
+            coeff_shift] at hc
+          have hzero : (Zero.zero : R) = 0 := rfl
+          simp only [hzero] at hc
+          have h₀ : lo + i - base₀ < hi - base₀ := by
+            dsimp [hi, base₀]
+            omega
+          have he₀ : base₀ + (lo + i - base₀) = lo + i := by
+            dsimp [base₀]
+            omega
+          rw [_root_.ite_eq_left h₀, he₀]
+          by_cases hk : k ≤ lo + i
+          · have h₁ : lo + i - k - base₁ < (hi - k) - base₁ := by
+              dsimp [hi, base₁]
+              omega
+            have h₀' : lo + i - k - base₀ < hi - base₀ := by
+              dsimp [hi, base₀]
+              omega
+            have h₂ : lo + i - k - base₂ < (hi - k) - base₂ := by
+              dsimp [hi, base₂]
+              omega
+            have he₁ : base₁ + (lo + i - k - base₁) = lo + i - k := by
+              dsimp [base₁]
+              omega
+            have he₀' : base₀ + (lo + i - k - base₀) = lo + i - k := by
+              dsimp [base₀]
+              omega
+            have he₂ : base₂ + (lo + i - k - base₂) = lo + i - k := by
+              dsimp [base₂]
+              omega
+            rw [_root_.ite_eq_left hk, _root_.ite_eq_left h₁,
+              _root_.ite_eq_left h₀', _root_.ite_eq_left h₂]
+            by_cases h2k : 2 * k ≤ lo + i
+            · have h₂' : lo + i - 2 * k - base₂ < (hi - k) - base₂ := by
+                dsimp [hi, base₂]
+                omega
+              have he₂' : base₂ + (lo + i - 2 * k - base₂) =
+                  lo + i - 2 * k := by
+                dsimp [base₂]
+                omega
+              rw [_root_.ite_eq_left h2k, _root_.ite_eq_left h₂', he₂']
+              simp only [_root_.ite_true]
+              rw [_root_.if_neg (by omega), _root_.if_neg (by omega)] at hc
+              grind
+            · rw [_root_.ite_eq_right h2k]
+              simp only [_root_.ite_true]
+              rw [_root_.if_neg (by omega), _root_.if_pos (by omega)] at hc
+              grind
+          · have h2k : ¬2 * k ≤ lo + i := by omega
+            rw [_root_.ite_eq_right hk, _root_.ite_eq_right h2k]
+            simp only [_root_.ite_true]
+            rw [_root_.if_pos (by omega), _root_.if_pos (by omega)] at hc
+            grind
+        · have hlen :
+              (List.map
+                (fun i =>
+                  let k := (max a.size b.size + 1) / 2
+                  let a₀ := low k a
+                  let a₁ := high k a
+                  let b₀ := low k b
+                  let b₁ := high k b
+                  let hi := lo + len
+                  let base₀ := lo - k
+                  let base₁ := lo - k
+                  let base₂ := lo - 2 * k
+                  let z₀ := karatsubaSliceAux cutoff fuel base₀ (hi - base₀) a₀ b₀
+                  let z₁ := karatsubaSliceAux cutoff fuel base₁
+                    ((hi - k) - base₁) (a₀ + a₁) (b₀ + b₁)
+                  let z₂ := karatsubaSliceAux cutoff fuel base₂
+                    ((hi - k) - base₂) a₁ b₁
+                  let d := lo + i
+                  z₀.coeff (d - base₀) +
+                    (if k ≤ d then
+                      z₁.coeff (d - k - base₁) - z₀.coeff (d - k - base₀) -
+                        z₂.coeff (d - k - base₂)
+                    else 0) +
+                    if 2 * k ≤ d then z₂.coeff (d - 2 * k - base₂) else 0)
+                (List.range len)).length ≤ i := by simp; omega
+          rw [List.getD_eq_getElem?_getD]
+          simp [hil]
+          rfl
+
+/-- Interval-pruned Karatsuba slicing. -/
+def karatsubaSlice (cutoff lo len : Nat) (a b : DensePoly R) : DensePoly R :=
+  karatsubaSliceAux cutoff (max a.size b.size) lo len a b
+
+/-- Interval-pruned Karatsuba slicing has the exact planned-slice semantics. -/
+theorem coeff_karatsubaSlice (cutoff lo len : Nat) (a b : DensePoly R)
+    (i : Nat) :
+    (karatsubaSlice cutoff lo len a b).coeff i =
+      if i < len then (a * b).coeff (lo + i) else 0 :=
+  coeff_karatsubaSliceAux cutoff _ lo len a b i
+
+/-- A lawful Karatsuba plan. -/
+def karatsubaPlan (cutoff : Nat) : MulPlan R where
+  mul := mulKaratsuba cutoff
+  square := squareKaratsuba cutoff
+  slice := karatsubaSlice cutoff
+  mul_eq := mulKaratsuba_eq cutoff
+  square_eq := squareKaratsuba_eq cutoff
+  coeff_slice := coeff_karatsubaSlice cutoff
+
+end Hex.DensePoly
diff --git a/HexPolyFast/Plan.lean b/HexPolyFast/Plan.lean
new file mode 100644
index 000000000..4d13ae63b
--- /dev/null
+++ b/HexPolyFast/Plan.lean
@@ -0,0 +1,171 @@
+/-
+Copyright (c) 2026 Lean FRO, LLC. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Kim Morrison
+-/
+
+module
+
+public import HexPoly
+
+public section
+
+/-!
+Explicit multiplication plans for dense polynomials.
+
+The plan laws keep optimized kernels behind the existing schoolbook
+`DensePoly` semantics.  Plans are values rather than typeclass instances, so
+callers can select a coefficient-specific implementation locally.
+-/
+
+namespace Hex.DensePoly
+
+universe u
+
+attribute [local instance 1000] Lean.Grind.Semiring.ofNat
+
+/-- A proof-carrying implementation of full, square, and clipped dense
+polynomial multiplication.  `slice lo len a b` stores coefficients beginning
+at degree `lo`, shifted down to degree zero. -/
+structure MulPlan (R : Type u) [DecidableEq R] [Lean.Grind.CommRing R] where
+  /-- A complete normalized product. -/
+  mul : DensePoly R → DensePoly R → DensePoly R
+  /-- A specialized square. -/
+  square : DensePoly R → DensePoly R
+  /-- `len` coefficients beginning at product degree `lo`, shifted down. -/
+  slice : Nat → Nat → DensePoly R → DensePoly R → DensePoly R
+  mul_eq : ∀ a b, mul a b = a * b
+  square_eq : ∀ a, square a = a * a
+  coeff_slice : ∀ lo len a b i,
+    (slice lo len a b).coeff i =
+      if i < len then (a * b).coeff (lo + i) else 0
+
+variable {R : Type u} [DecidableEq R] [Lean.Grind.CommRing R]
+
+/-- Kernel-facing clipped schoolbook product.  The definition is phrased by
+the semantic product so its coefficient law is immediate; compiled code uses
+`schoolbookSliceImpl`, which evaluates the requested schoolbook coefficient
+folds directly and never materializes the full product. -/
+@[expose]
+noncomputable def schoolbookSlice (lo len : Nat) (a b : DensePoly R) :
+    DensePoly R :=
+  ofList ((List.range len).map fun i => (a * b).coeff (lo + i))
+
+/-- Allocation-conscious runtime implementation of `schoolbookSlice`. -/
+@[expose]
+def schoolbookSliceImpl (lo len : Nat) (a b : DensePoly R) : DensePoly R :=
+  ofList ((List.range len).map fun i => mulCoeffSum a b (lo + i))
+
+/-- The direct coefficient-fold implementation agrees with the semantic
+clipped product. -/
+theorem schoolbookSlice_eq_impl (lo len : Nat) (a b : DensePoly R) :
+    schoolbookSlice lo len a b = schoolbookSliceImpl lo len a b := by
+  unfold schoolbookSlice schoolbookSliceImpl
+  apply congrArg ofList
+  apply List.map_congr_left
+  intro i hi
+  exact coeff_mul a b (lo + i)
+
+/-- Compiled clipped schoolbook products use direct coefficient folds. -/
+@[csimp]
+theorem schoolbookSlice_csimp : @schoolbookSlice = @schoolbookSliceImpl := by
+  funext R instDecEq instRing lo len a b
+  exact schoolbookSlice_eq_impl lo len a b
+
+/-- Coefficient law for a clipped schoolbook product. -/
+theorem coeff_schoolbookSlice (lo len : Nat) (a b : DensePoly R) (i : Nat) :
+    (schoolbookSlice lo len a b).coeff i =
+      if i < len then (a * b).coeff (lo + i) else 0 := by
+  unfold schoolbookSlice
+  rw [coeff_ofList]
+  by_cases hi : i < len
+  · simp [List.getD, hi]
+  · have hlen : (List.map (fun j => (a * b).coeff (lo + j))
+        (List.range len)).length ≤ i := by simp; omega
+    rw [List.getD_eq_getElem?_getD]
+    simp [hi]
+    change (Zero.zero : R) = (Zero.zero : R)
+    rfl
+
+/-- The reference plan backed by the existing allocation-conscious
+schoolbook multiplication. -/
+def schoolbookPlan : MulPlan R where
+  mul := mulImpl
+  square := fun a => mulImpl a a
+  slice := schoolbookSlice
+  mul_eq := fun a b => (mul_eq_mulImpl a b).symm
+  square_eq := fun a => (mul_eq_mulImpl a a).symm
+  coeff_slice := coeff_schoolbookSlice
+
+/-- Multiply using an explicit plan. -/
+@[inline]
+def mulWith (plan : MulPlan R) (a b : DensePoly R) : DensePoly R :=
+  plan.mul a b
+
+/-- Square using an explicit plan. -/
+@[inline]
+def squareWith (plan : MulPlan R) (a : DensePoly R) : DensePoly R :=
+  plan.square a
+
+/-- Keep the first `len` coefficients of a planned product. -/
+@[inline]
+def mulLow (plan : MulPlan R) (len : Nat) (a b : DensePoly R) : DensePoly R :=
+  plan.slice 0 len a b
+
+/-- Keep `len` coefficients beginning at product degree `lo`. -/
+@[inline]
+def mulSlice (plan : MulPlan R) (lo len : Nat) (a b : DensePoly R) :
+    DensePoly R :=
+  plan.slice lo len a b
+
+/-- The standard middle product for operands of sizes `m ≥ n > 0`.
+The result contains product degrees `n - 1` through `m - 1`, shifted down. -/
+@[inline]
+def mulMiddle (plan : MulPlan R) (a b : DensePoly R)
+    (_hsize : b.size ≤ a.size) (_hpos : 0 < b.size) : DensePoly R :=
+  plan.slice (b.size - 1) (a.size - b.size + 1) a b
+
+/-- Checked middle product.  Operands are ordered by size; an empty operand
+has zero middle product. -/
+def mulMiddleChecked (plan : MulPlan R) (a b : DensePoly R) : DensePoly R :=
+  if ha : a.size = 0 then 0
+  else if hb : b.size = 0 then 0
+  else
+    if hba : b.size ≤ a.size then
+      mulMiddle plan a b hba (by omega)
+    else
+      mulMiddle plan b a (Nat.le_of_lt (Nat.lt_of_not_ge hba)) (by omega)
+
+/-- Planned multiplication has the existing dense-polynomial semantics. -/
+theorem mulWith_eq (plan : MulPlan R) (a b : DensePoly R) :
+    mulWith plan a b = a * b :=
+  plan.mul_eq a b
+
+/-- Planned squaring has the existing dense-polynomial semantics. -/
+theorem squareWith_eq (plan : MulPlan R) (a : DensePoly R) :
+    squareWith plan a = a * a :=
+  plan.square_eq a
+
+/-- Coefficient law for a planned low product. -/
+theorem coeff_mulLow (plan : MulPlan R) (len i : Nat) (a b : DensePoly R) :
+    (mulLow plan len a b).coeff i =
+      if i < len then (a * b).coeff i else 0 := by
+  simpa [mulLow] using plan.coeff_slice 0 len a b i
+
+/-- Coefficient law for an arbitrary planned slice. -/
+theorem coeff_mulSlice (plan : MulPlan R) (lo len i : Nat)
+    (a b : DensePoly R) :
+    (mulSlice plan lo len a b).coeff i =
+      if i < len then (a * b).coeff (lo + i) else 0 :=
+  plan.coeff_slice lo len a b i
+
+/-- Coefficient law for the standard proof-taking middle product. -/
+theorem coeff_mulMiddle (plan : MulPlan R) (a b : DensePoly R)
+    (hsize : b.size ≤ a.size) (hpos : 0 < b.size) (i : Nat) :
+    (mulMiddle plan a b hsize hpos).coeff i =
+      if i < a.size - b.size + 1 then
+        (a * b).coeff (b.size - 1 + i)
+      else 0 :=
+  plan.coeff_slice (b.size - 1) (a.size - b.size + 1) a b i
+
+end Hex.DensePoly
diff --git a/HexPolyFast/Reciprocal.lean b/HexPolyFast/Reciprocal.lean
new file mode 100644
index 000000000..2261d8a5a
--- /dev/null
+++ b/HexPolyFast/Reciprocal.lean
@@ -0,0 +1,157 @@
+/-
+Copyright (c) 2026 Lean FRO, LLC. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Kim Morrison
+-/
+
+module
+
+public import HexPolyFast.Reverse
+public import HexTruncatedSeries.Inverse
+
+public section
+
+/-!
+Plan-driven Newton reciprocals for fixed-precision truncated series.
+
+Every bounded product in the Newton loop goes through `MulPlan.slice`.  The
+agreement theorem below connects that executable polynomial kernel to the
+existing `TSeries.mulUpTo` semantics.
+-/
+
+namespace Hex.DensePoly
+
+universe u
+
+attribute [local instance 1000] Lean.Grind.Semiring.ofNat
+
+variable {R : Type u} [DecidableEq R] [Lean.Grind.CommRing R]
+
+private theorem fold_diagonal_extend (p q : DensePoly R) (n d : Nat) :
+    (List.range (p.size + d)).foldl
+        (fun acc i => acc + diagonalMulCoeffTerm p q n i) 0 =
+      (List.range p.size).foldl
+        (fun acc i => acc + diagonalMulCoeffTerm p q n i) 0 := by
+  induction d with
+  | zero => simp
+  | succ d ih =>
+      rw [Nat.add_succ, List.range_succ, List.foldl_append]
+      simp only [List.foldl_cons, List.foldl_nil]
+      rw [ih]
+      have hcoeff : p.coeff (p.size + d) = 0 :=
+        coeff_eq_zero_of_size_le p (by omega)
+      simp [diagonalMulCoeffTerm, hcoeff]
+      grind
+
+private theorem fold_diagonal_truncate (p q : DensePoly R) (n d : Nat) :
+    (List.range (n + 1 + d)).foldl
+        (fun acc i => acc + diagonalMulCoeffTerm p q n i) 0 =
+      (List.range (n + 1)).foldl
+        (fun acc i => acc + diagonalMulCoeffTerm p q n i) 0 := by
+  induction d with
+  | zero => simp
+  | succ d ih =>
+      rw [Nat.add_succ, List.range_succ, List.foldl_append]
+      simp only [List.foldl_cons, List.foldl_nil]
+      rw [ih]
+      simp [diagonalMulCoeffTerm]
+      grind
+
+/-- A polynomial multiplication diagonal may always be normalized to the
+canonical degree-sized range, independently of the left operand's support. -/
+theorem diagonal_eq_degree_bound (p q : DensePoly R) (n : Nat) :
+    (List.range p.size).foldl
+        (fun acc i => acc + diagonalMulCoeffTerm p q n i) 0 =
+      (List.range (n + 1)).foldl
+        (fun acc i => acc + diagonalMulCoeffTerm p q n i) 0 := by
+  by_cases hp : p.size ≤ n + 1
+  · have hsum : p.size + (n + 1 - p.size) = n + 1 := by omega
+    rw [← hsum]
+    exact (fold_diagonal_extend p q n (n + 1 - p.size)).symm
+  · have hsum : n + 1 + (p.size - (n + 1)) = p.size := by omega
+    rw [← hsum]
+    exact fold_diagonal_truncate p q n (p.size - (n + 1))
+
+/-- Polynomial conversion preserves every represented product coefficient. -/
+theorem coeff_polyOfSeries_mul {n : Nat} (a b : TSeries R n)
+    (i : Nat) (hi : i < n) :
+    (polyOfSeries a * polyOfSeries b).coeff i = (a * b).coeff i := by
+  rw [DensePoly.coeff_mul, mulCoeffSum_eq_diagonal,
+    diagonal_eq_degree_bound, TSeries.coeff_mul a b i hi]
+  unfold TSeries.convCoeff
+  apply List.foldl_add_congr
+  intro j hj
+  have hjle : j ≤ i := by
+    have := List.mem_range.mp hj
+    omega
+  unfold diagonalMulCoeffTerm
+  rw [_root_.if_neg (by omega), coeff_polyOfSeries_of_lt,
+    coeff_polyOfSeries_of_lt]
+  all_goals omega
+
+/-- Multiply two series through degree `m - 1` with a polynomial
+multiplication plan. -/
+def seriesMulUpTo (plan : MulPlan R) {n : Nat} (m : Nat)
+    (a b : TSeries R n) : TSeries R n :=
+  let p := mulLow plan m (polyOfSeries a) (polyOfSeries b)
+  TSeries.ofFn fun i => p.coeff i
+
+/-- Planned bounded series multiplication has exactly the established
+`TSeries.mulUpTo` semantics. -/
+theorem seriesMulUpTo_eq (plan : MulPlan R) {n : Nat} (m : Nat)
+    (a b : TSeries R n) :
+    seriesMulUpTo plan m a b = TSeries.mulUpTo m a b := by
+  apply TSeries.ext
+  intro i hi
+  unfold seriesMulUpTo
+  rw [TSeries.coeff_ofFn _ i hi, coeff_mulLow,
+    TSeries.coeff_mulUpTo m a b i hi]
+  split
+  · rw [coeff_polyOfSeries_mul a b i hi]
+  · rfl
+
+/-- One Newton reciprocal update whose two bounded products both use the
+supplied multiplication plan. -/
+def reciprocalStep (plan : MulPlan R) {n : Nat} (g h : TSeries R n)
+    (m : Nat) : TSeries R n :=
+  seriesMulUpTo plan m h
+    (TSeries.C (1 + 1) - seriesMulUpTo plan m g h)
+
+/-- A planned Newton update agrees exactly with the reference series update. -/
+theorem reciprocalStep_eq (plan : MulPlan R) {n : Nat}
+    (g h : TSeries R n) (m : Nat) :
+    reciprocalStep plan g h m = TSeries.invStep g h m := by
+  unfold reciprocalStep TSeries.invStep
+  rw [seriesMulUpTo_eq, seriesMulUpTo_eq]
+
+/-- Newton reciprocal at the full represented precision.  Each doubling step
+uses `plan.slice 0 k`; no schoolbook `TSeries` product is executed by this
+definition. -/
+def reciprocalWith (plan : MulPlan R) {n : Nat}
+    (g : TSeries R n) (u : R) : TSeries R n :=
+  TSeries.newton (reciprocalStep plan g) (TSeries.C u) (TSeries.steps n)
+
+/-- The plan-driven reciprocal has the established truncated-series
+semantics.  The inverse-witness hypothesis is intentionally not needed for
+algorithmic agreement; it is needed by the defining inverse equation below. -/
+theorem reciprocalWith_eq (plan : MulPlan R) {n : Nat}
+    (g : TSeries R n) (u : R) :
+    reciprocalWith plan g u = TSeries.invOfUnit g u := by
+  have hstep : reciprocalStep plan g = TSeries.invStep g := by
+    funext h m
+    exact reciprocalStep_eq plan g h m
+  unfold reciprocalWith TSeries.invOfUnit TSeries.invUpTo
+  simp only [Nat.min_self, hstep]
+  apply TSeries.ext
+  intro i hi
+  rw [TSeries.coeff_ofFn _ i hi, _root_.ite_eq_left hi]
+
+/-- A planned reciprocal satisfies the multiplicative inverse equation when
+the supplied constant coefficient really is an inverse. -/
+theorem reciprocalWith_mul (plan : MulPlan R) {n : Nat}
+    (g : TSeries R n) (u : R) (hu : g.coeff 0 * u = 1) :
+    g * reciprocalWith plan g u = 1 := by
+  rw [reciprocalWith_eq]
+  exact TSeries.invOfUnit_mul g u hu
+
+end Hex.DensePoly
diff --git a/HexPolyFast/Reverse.lean b/HexPolyFast/Reverse.lean
new file mode 100644
index 000000000..f18d77019
--- /dev/null
+++ b/HexPolyFast/Reverse.lean
@@ -0,0 +1,89 @@
+/-
+Copyright (c) 2026 Lean FRO, LLC. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Kim Morrison
+-/
+
+module
+
+public import HexPolyFast.Plan
+public import HexTruncatedSeries
+
+public section
+
+/-!
+Coefficientwise bridges between normalized dense polynomials and fixed-length
+truncated series.  Reversal carries an explicit zero-extension guard so Nat
+subtraction cannot alias out-of-range positions to the constant coefficient.
+-/
+
+namespace Hex.DensePoly
+
+universe u
+
+variable {R : Type u} [Zero R] [DecidableEq R]
+
+/-- Reverse a dense polynomial from its leading end into an exactly `n`-term
+truncated series, padding with zeros when `n` exceeds the polynomial size. -/
+@[expose]
+def reverseSeries (f : DensePoly R) (n : Nat) : TSeries R n :=
+  TSeries.ofFn fun i =>
+    if i < f.size then f.coeff (f.size - 1 - i) else 0
+
+/-- Coefficient law for zero-extending polynomial reversal. -/
+theorem coeff_reverseSeries (f : DensePoly R) (n i : Nat) (hi : i < n) :
+    (reverseSeries f n).coeff i =
+      if i < f.size then f.coeff (f.size - 1 - i) else 0 := by
+  exact TSeries.coeff_ofFn _ i hi
+
+/-- Inside the polynomial support, reversal reads down from the leading
+coefficient. -/
+theorem coeff_reverseSeries_of_lt (f : DensePoly R) (n i : Nat)
+    (hi : i < n) (hf : i < f.size) :
+    (reverseSeries f n).coeff i = f.coeff (f.size - 1 - i) := by
+  rw [coeff_reverseSeries f n i hi, _root_.ite_eq_left hf]
+
+/-- Outside the polynomial support, reversal is zero rather than the constant
+coefficient selected by saturated Nat subtraction. -/
+theorem coeff_reverseSeries_of_size_le (f : DensePoly R) (n i : Nat)
+    (hi : i < n) (hf : f.size ≤ i) :
+    (reverseSeries f n).coeff i = 0 := by
+  rw [coeff_reverseSeries f n i hi, _root_.ite_eq_right (Nat.not_lt.mpr hf)]
+
+/-- Convert all represented coefficients of a truncated series to a normalized
+dense polynomial. -/
+@[expose]
+def polyOfSeries {n : Nat} (a : TSeries R n) : DensePoly R :=
+  ofList ((List.range n).map a.coeff)
+
+/-- Coefficient law for conversion from a fixed series prefix. -/
+theorem coeff_polyOfSeries {n : Nat} (a : TSeries R n) (i : Nat) :
+    (polyOfSeries a).coeff i = if i < n then a.coeff i else 0 := by
+  unfold polyOfSeries
+  rw [coeff_ofList]
+  by_cases hi : i < n
+  · simp [List.getD, hi]
+  · rw [List.getD_eq_getElem?_getD]
+    simp [hi]
+    rfl
+
+/-- Converting a series to a polynomial and reading a represented coefficient
+returns the original coefficient. -/
+theorem coeff_polyOfSeries_of_lt {n : Nat} (a : TSeries R n) (i : Nat)
+    (hi : i < n) :
+    (polyOfSeries a).coeff i = a.coeff i := by
+  rw [coeff_polyOfSeries, _root_.ite_eq_left hi]
+
+/-- Converting a reversed polynomial prefix back to a polynomial exposes the
+same guarded reversal coefficient law. -/
+theorem coeff_polyOf_reverseSeries (f : DensePoly R) (n i : Nat) :
+    (polyOfSeries (reverseSeries f n)).coeff i =
+      if i < n then
+        if i < f.size then f.coeff (f.size - 1 - i) else 0
+      else 0 := by
+  rw [coeff_polyOfSeries]
+  split <;> rename_i h
+  · exact coeff_reverseSeries f n i h
+  · rfl
+
+end Hex.DensePoly
diff --git a/HexPolyFast/Tree.lean b/HexPolyFast/Tree.lean
new file mode 100644
index 000000000..5178b43e0
--- /dev/null
+++ b/HexPolyFast/Tree.lean
@@ -0,0 +1,189 @@
+/-
+Copyright (c) 2026 Lean FRO, LLC. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Kim Morrison
+-/
+
+module
+
+public import HexPolyFast.Division
+
+public section
+set_option backward.proofsInPublic true
+
+/-!
+Balanced product trees.
+
+The representation is private: clients observe the original leaves, the
+balanced levels, and the root product through the accessors below. Adjacent
+nodes are multiplied with the supplied lawful plan; an unpaired final node is
+carried to the next level unchanged.
+-/
+
+namespace Hex.DensePoly
+
+universe u
+
+variable {R : Type u} [DecidableEq R] [Lean.Grind.CommRing R]
+
+/-- Multiply adjacent entries, carrying an unpaired final entry unchanged. -/
+private def pairProducts (plan : MulPlan R) : List (DensePoly R) → List (DensePoly R)
+  | [] => []
+  | [p] => [p]
+  | p :: q :: rest => mulWith plan p q :: pairProducts plan rest
+
+/-- Product of a list through a lawful multiplication plan. -/
+private def plannedProduct (plan : MulPlan R) (xs : List (DensePoly R)) : DensePoly R :=
+  xs.foldl (fun acc p => mulWith plan acc p) 1
+
+private theorem plannedProduct_cons (plan : MulPlan R) (p : DensePoly R)
+    (xs : List (DensePoly R)) :
+    plannedProduct plan (p :: xs) = p * plannedProduct plan xs := by
+  unfold plannedProduct
+  simp only [List.foldl_cons]
+  have hfun : (fun acc q => mulWith plan acc q) = (fun acc q => acc * q) := by
+    funext acc q
+    exact mulWith_eq plan acc q
+  rw [mulWith_eq, DensePoly.mul_comm_poly (1 : DensePoly R) p,
+    DensePoly.mul_one_right_poly, hfun]
+  simpa using List.foldl_mul_eq_mul_foldl xs id p
+
+private theorem plannedProduct_pairProducts (plan : MulPlan R) :
+    ∀ xs, plannedProduct plan (pairProducts plan xs) = plannedProduct plan xs := by
+  intro xs
+  induction hlen : xs.length using Nat.strongRecOn generalizing xs with
+  | ind n ih =>
+      cases xs with
+      | nil => rfl
+      | cons p rest =>
+          cases rest with
+          | nil => rfl
+          | cons q rest =>
+              rw [pairProducts, plannedProduct_cons, mulWith_eq,
+                plannedProduct_cons, plannedProduct_cons,
+                ih rest.length (by simp at hlen; omega) rest rfl,
+                DensePoly.mul_assoc_poly]
+
+/-- Build all nonempty balanced levels and return their root. Fuel is only a
+totality guard; construction supplies the leaf count. -/
+private def buildLevels (plan : MulPlan R) :
+    Nat → List (DensePoly R) → List (List (DensePoly R)) × DensePoly R
+  | 0, xs => ([xs], plannedProduct plan xs)
+  | _ + 1, [] => ([[1]], 1)
+  | _ + 1, [p] => ([[p]], p)
+  | fuel + 1, p :: q :: rest =>
+      let current := p :: q :: rest
+      let next := pairProducts plan current
+      let built := buildLevels plan fuel next
+      (current :: built.1, built.2)
+
+private theorem buildLevels_root (plan : MulPlan R) : ∀ fuel xs,
+    (buildLevels plan fuel xs).2 = plannedProduct plan xs := by
+  intro fuel
+  induction fuel with
+  | zero => intro xs; rfl
+  | succ fuel ih =>
+      intro xs
+      cases xs with
+      | nil => simp [buildLevels, plannedProduct, mulWith_eq]
+      | cons p rest =>
+          cases rest with
+          | nil =>
+              simp [buildLevels, plannedProduct, mulWith_eq,
+                DensePoly.mul_comm_poly (1 : DensePoly R) p,
+                DensePoly.mul_one_right_poly]
+          | cons q rest =>
+              rw [buildLevels]
+              dsimp only
+              rw [ih, plannedProduct_pairProducts]
+
+/-- An opaque balanced product tree. -/
+structure ProductTree (R : Type u) [DecidableEq R] [Lean.Grind.CommRing R] where
+  private planData : MulPlan R
+  private leafData : Array (DensePoly R)
+  private levelData : Array (Array (DensePoly R))
+  private rootData : DensePoly R
+
+namespace ProductTree
+
+/-- Build a balanced product tree. The empty tree has root `1` and one
+singleton internal level containing that root. -/
+def build (plan : MulPlan R) (leaves : Array (DensePoly R)) : ProductTree R :=
+  let built := buildLevels plan leaves.size leaves.toList
+  { planData := plan
+    leafData := leaves
+    levelData := (built.1.map List.toArray).toArray
+    rootData := built.2 }
+
+/-- The leaf sequence in its original order. -/
+def leaves (tree : ProductTree R) : Array (DensePoly R) := tree.leafData
+
+/-- Number of stored nonempty levels. -/
+def levelCount (tree : ProductTree R) : Nat := tree.levelData.size
+
+/-- Lawful multiplication plan used to build the tree. -/
+def plan (tree : ProductTree R) : MulPlan R := tree.planData
+
+/-- A stored balanced level, from leaves upward. -/
+def level? (tree : ProductTree R) (i : Nat) : Option (Array (DensePoly R)) :=
+  tree.levelData[i]?
+
+/-- The product represented by the root node. -/
+def root (tree : ProductTree R) : DensePoly R := tree.rootData
+
+/-- The ordered leaf block represented by a node. Invalid level/node pairs
+represent the empty block; use `nodeProduct?` when validity matters. -/
+def nodeLeaves (tree : ProductTree R) (level index : Nat) : Array (DensePoly R) :=
+  let width := 2 ^ level
+  let lo := index * width
+  tree.leafData.extract lo (min tree.leafData.size (lo + width))
+
+/-- Product represented by a valid stored node. The observation is semantic:
+it folds exactly that node's leaf block, independently of the internal level
+layout. -/
+def nodeProduct? (tree : ProductTree R) (level index : Nat) : Option (DensePoly R) :=
+  if level < tree.levelCount &&
+      ((tree.leafData.isEmpty && level = 0 && index = 0) ||
+        index * 2 ^ level < tree.leafData.size) then
+    some ((tree.nodeLeaves level index).foldl
+      (fun acc p => mulWith tree.plan acc p) 1)
+  else
+    none
+
+/-- Every observed node product is the planned product of precisely its
+represented leaf block. -/
+theorem nodeProduct?_eq (tree : ProductTree R) (level index : Nat)
+    (p : DensePoly R) (h : tree.nodeProduct? level index = some p) :
+    p = (tree.nodeLeaves level index).foldl
+      (fun acc q => mulWith tree.plan acc q) 1 := by
+  unfold nodeProduct? at h
+  split at h
+  · exact Option.some.inj h |>.symm
+  · contradiction
+
+/-- Building preserves the supplied leaf sequence. -/
+@[simp] theorem leaves_build (plan : MulPlan R) (leaves : Array (DensePoly R)) :
+    ProductTree.leaves (build plan leaves) = leaves := by
+  rfl
+
+/-- The root is the planned product of all leaves in order. -/
+theorem root_build (plan : MulPlan R) (leaves : Array (DensePoly R)) :
+    (build plan leaves).root =
+      leaves.toList.foldl (fun acc p => mulWith plan acc p) 1 := by
+  unfold build root
+  dsimp only
+  exact buildLevels_root plan leaves.size leaves.toList
+
+/-- The root is independent of the selected lawful multiplication kernel. -/
+theorem root_build_eq_foldl (plan : MulPlan R) (leaves : Array (DensePoly R)) :
+    (build plan leaves).root = leaves.foldl (fun acc p => acc * p) 1 := by
+  rw [root_build]
+  rw [← Array.foldl_toList]
+  have hfun : (fun acc p => mulWith plan acc p) = (fun acc p => acc * p) := by
+    funext acc p
+    exact mulWith_eq plan acc p
+  rw [hfun]
+
+end ProductTree
+
+end Hex.DensePoly
diff --git a/SPEC/Libraries/hex-poly-fast.md b/SPEC/Libraries/hex-poly-fast.md
index ea7188839..b00154246 100644
--- a/SPEC/Libraries/hex-poly-fast.md
+++ b/SPEC/Libraries/hex-poly-fast.md
@@ -122,7 +122,7 @@ There is deliberately no typeclass for the plan. Plans are passed explicitly:
 ```lean
 namespace Hex.DensePoly
 
-structure MulPlan (R : Type u) [Zero R] [DecidableEq R]
+structure MulPlan (R : Type u) [DecidableEq R]
     [Lean.Grind.CommRing R] where
   /-- A complete normalized product. -/
   mul : DensePoly R → DensePoly R → DensePoly R
@@ -150,6 +150,11 @@ def mulLow (plan : MulPlan R) (len : Nat)
   plan.slice 0 len a b
 ```
 
+`Zero R` is derived from the commutative-ring numeral structure rather than
+accepted as an independent parameter.  This matters for normalized dense
+polynomials: a separately supplied `Zero` could disagree with the ring's
+additive identity, invalidating both trimming and Karatsuba subtraction.
+
 The proof fields are erased. Passing a plan does not move correctness into a
 runtime checker, and using a structure rather than a typeclass prevents a
 global choice of integer or finite-field multiplication from leaking into
@@ -243,12 +248,14 @@ the leading end first and zero-extends to exactly `n` entries:
 def reverseSeries (f : DensePoly R) (n : Nat) : TSeries R n
 
 theorem coeff_reverseSeries (f : DensePoly R) (n i : Nat) (hi : i < n) :
-    (reverseSeries f n).coeff i = f.coeff (f.size - 1 - i)
+    (reverseSeries f n).coeff i =
+      if i < f.size then f.coeff (f.size - 1 - i) else 0
 ```
 
-The theorem is accompanied by a range-aware form, because subtraction on
-`Nat` makes the displayed convenience statement useful only while
-`i < f.size`. `polyOfSeries` converts a fixed prefix back through
+The theorem is accompanied by a range form that removes the conditional under
+`i < f.size`.  The guard is essential: subtraction on `Nat` saturates at zero,
+so the unguarded right-hand side would read the constant coefficient rather
+than zero once `i` passes the leading end. `polyOfSeries` converts a fixed prefix back through
 `DensePoly.ofCoeffs`; its coefficient theorem, not structural array equality,
 is the bridge used by division proofs.
 
@@ -270,7 +277,7 @@ A reciprocal is worth caching whenever a fixed modulus divides many values,
 as in a remainder tree or quotient ring:
 
 ```lean
-structure DivPlan (R : Type u) [Zero R] [DecidableEq R]
+structure DivPlan (R : Type u) [DecidableEq R]
     [Lean.Grind.CommRing R] where
   divisor : DensePoly R
   capacity : Nat
@@ -279,7 +286,7 @@ structure DivPlan (R : Type u) [Zero R] [DecidableEq R]
   reciprocal_spec : ...
 
 def DivPlan.ofMonic (mul : MulPlan R) (q : DensePoly R)
-    (hq : Monic q) (capacity : Nat) : DivPlan R
+    (hq : Monic q) (hqne : q ≠ 0) (capacity : Nat) : DivPlan R
 
 def DivPlan.ofNonzero [Lean.Grind.Field R] (mul : MulPlan R)
     (q : DensePoly R) (hq : q ≠ 0) (capacity : Nat) : DivPlan R
@@ -317,7 +324,7 @@ The half-gcd implementation uses a dedicated four-polynomial transformation
 record rather than depending on hex-matrix:
 
 ```lean
-structure GcdStep (R : Type u) [Zero R] [DecidableEq R] where
+structure GcdStep (R : Type u) [DecidableEq R] where
   a00 : DensePoly R
   a01 : DensePoly R
   a10 : DensePoly R
@@ -360,8 +367,14 @@ theorem EvalPlan.get_eval (plan : EvalPlan R) (f) (i) (hi : i < plan.size) :
     (plan.eval f)[i] = f.eval plan.points[i]
 ```
 
-Evaluation works over a commutative ring because every divisor in the point
-tree is monic. The empty point sequence produces an empty result.
+The cached remainder-tree path applies when `f.size <= plan.size`; this is the
+finite capacity determined when the plan is built. `EvalPlan.eval` remains
+total for larger inputs and uses direct pointwise evaluation in that case.
+This fallback is necessary because the signature accepts polynomials of
+unbounded size: no finite plan built from only the points can cache the
+unbounded reciprocal precision required to reduce every such input at the
+root. Evaluation works over a commutative ring because every divisor in the
+point tree is monic. The empty point sequence produces an empty result.
 
 Interpolation needs a field and distinct points. `InterpPlan.build?` returns
 `none` exactly when duplicate points are present. It reuses the point product,
@@ -372,9 +385,10 @@ value array to zero.
 
 Soundness states that the result has size at most the point count and evaluates
 to every supplied value. Uniqueness states that any polynomial of smaller
-degree with those values is equal to the result. Both construction and one
-evaluation/interpolation cost `O(M(n) log n)`; the reusable plan cost is
-reported separately.
+degree with those values is equal to the result. Construction, bounded-size
+evaluation, and interpolation cost `O(M(n) log n)`; the oversized
+direct-evaluation fallback costs `O(plan.size * f.size)`. The reusable plan
+cost is reported separately.
 
 ## Padé approximation
 
diff --git a/bench/HexPolyFast/Bench.lean b/bench/HexPolyFast/Bench.lean
new file mode 100644
index 000000000..a8b2ea7d8
--- /dev/null
+++ b/bench/HexPolyFast/Bench.lean
@@ -0,0 +1,196 @@
+/-
+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 HexPolyFast
+import LeanBench
+
+/-!
+Scientific benchmark registrations for the generic multiplication crossover.
+Input preparation is excluded from timing, and every target consumes its
+result through a coefficient hash.
+-/
+
+namespace Hex.PolyFastBench
+
+open Hex Hex.DensePoly
+
+structure Binary where
+  left : DensePoly Int
+  right : DensePoly Int
+
+instance : Hashable Binary where
+  hash input := mixHash (hash input.left.toArray) (hash input.right.toArray)
+
+private def coeff (i salt : Nat) : Int :=
+  Int.ofNat (((i + 3) * (salt + 11)) % 101 + 1) - 50
+
+def prepBalanced (n : Nat) : Binary :=
+  { left := ofList ((List.range n).map fun i => coeff i 3)
+    right := ofList ((List.range n).map fun i => coeff i 19) }
+
+def prepSkew (n : Nat) : Binary :=
+  { left := ofList ((List.range (64 * n)).map fun i => coeff i 5)
+    right := ofList ((List.range n).map fun i => coeff i 23) }
+
+private def checksum (p : DensePoly Int) : UInt64 :=
+  p.toArray.foldl (fun acc x => mixHash acc (hash x)) 0
+
+def runSchoolbook (input : Binary) : UInt64 :=
+  checksum (mulWith schoolbookPlan input.left input.right)
+
+def runKaratsuba (input : Binary) : UInt64 :=
+  checksum (mulWith (karatsubaPlan 32) input.left input.right)
+
+def runKaratsubaSkew (input : Binary) : UInt64 :=
+  checksum (mulWith (karatsubaPlan 32) input.left input.right)
+
+def runKaratsubaSquare (input : Binary) : UInt64 :=
+  checksum (squareWith (karatsubaPlan 32) input.left)
+
+structure DivisionInput where
+  dividend : DensePoly Rat
+  divisor : DensePoly Rat
+
+instance : Hashable DivisionInput where
+  hash input := mixHash (hash input.dividend.toArray) (hash input.divisor.toArray)
+
+def prepDivision (n : Nat) : DivisionInput :=
+  { dividend := ofList ((List.range (2 * n + 1)).map fun i => (coeff i 31 : Rat))
+    divisor := ofList (((List.range n).map fun i => (coeff i 47 : Rat)) ++ [1]) }
+
+private def checksumRat (p : DensePoly Rat) : UInt64 :=
+  p.toArray.foldl (fun acc x => mixHash acc (hash x)) 0
+
+private def checksumDiv (qr : DensePoly Rat × DensePoly Rat) : UInt64 :=
+  mixHash (checksumRat qr.1) (checksumRat qr.2)
+
+def runLongDivision (input : DivisionInput) : UInt64 :=
+  checksumDiv (divMod input.dividend input.divisor)
+
+def runNewtonDivision (input : DivisionInput) : UInt64 :=
+  checksumDiv (divModWith (karatsubaPlan 32) input.dividend input.divisor)
+
+structure ProductTreeInput where
+  leaves : Array (DensePoly Int)
+
+instance : Hashable ProductTreeInput where
+  hash input := input.leaves.foldl
+    (fun acc p => mixHash acc (hash p.toArray)) 0
+
+def prepProductTree (n : Nat) : ProductTreeInput :=
+  { leaves := (List.range n).map (fun i => ofList [-(coeff i 59), 1]) |>.toArray }
+
+def runProductTree (input : ProductTreeInput) : UInt64 :=
+  checksum (ProductTree.build (karatsubaPlan 32) input.leaves).root
+
+/- Cost model: a balanced length-`n` schoolbook convolution evaluates one
+coefficient product for each input pair, hence `n²` ring multiplications. -/
+setup_benchmark runSchoolbook n => n ^ 2
+  with prep := prepBalanced
+  where {
+    paramFloor := 4
+    paramCeiling := 16384
+    paramSchedule := .custom #[4, 16, 64, 256, 1024, 4096, 16384]
+    maxSecondsPerCall := 3.0
+    targetInnerNanos := 200000000
+    signalFloorMultiplier := 1.0
+    tags := #["multiplication", "schoolbook", "balanced"]
+  }
+
+/- Cost model: balanced Karatsuba satisfies `T(n) = 3T(n/2) + O(n)`, hence
+`T(n) = Θ(n^(log₂ 3))`; `n * sqrt n` is its integer-valued surrogate.
+The nearby 31/32/33 rungs expose the fixed cutoff. -/
+setup_benchmark runKaratsuba n => n * (Nat.sqrt n)
+  with prep := prepBalanced
+  where {
+    paramFloor := 4
+    paramCeiling := 16384
+    paramSchedule := .custom #[4, 16, 31, 32, 33, 64, 256, 1024, 4096, 16384]
+    maxSecondsPerCall := 3.0
+    targetInnerNanos := 200000000
+    signalFloorMultiplier := 1.0
+    tags := #["multiplication", "karatsuba", "balanced"]
+  }
+
+/- Cost model: specialized squaring performs three recursive squares plus
+linear combination work, so it obeys the same `Θ(n^(log₂ 3))` recurrence;
+`n * sqrt n` is the integer-valued Karatsuba-range surrogate. -/
+setup_benchmark runKaratsubaSquare n => n * (Nat.sqrt n)
+  with prep := prepBalanced
+  where {
+    paramFloor := 4
+    paramCeiling := 16384
+    paramSchedule := .custom #[4, 16, 31, 32, 33, 64, 256, 1024, 4096, 16384]
+    maxSecondsPerCall := 3.0
+    targetInnerNanos := 200000000
+    signalFloorMultiplier := 1.0
+    tags := #["multiplication", "karatsuba", "square"]
+  }
+
+/- Cost model: the 64:1 dispatcher partitions the long operand into 64
+length-`n` blocks.  That constant factor preserves the balanced Karatsuba
+bound `Θ(n^(log₂ 3))`, represented by the `n * sqrt n` surrogate. -/
+setup_benchmark runKaratsubaSkew n => n * (Nat.sqrt n)
+  with prep := prepSkew
+  where {
+    paramFloor := 4
+    paramCeiling := 256
+    paramSchedule := .custom #[4, 8, 16, 32, 64, 128, 256]
+    maxSecondsPerCall := 3.0
+    targetInnerNanos := 200000000
+    signalFloorMultiplier := 1.0
+    tags := #["multiplication", "karatsuba", "ratio-64"]
+  }
+
+/- The long-division comparator eliminates one leading coefficient at a time
+and updates a linear suffix, giving a quadratic coefficient-operation model. -/
+setup_benchmark runLongDivision n => n ^ 2
+  with prep := prepDivision
+  where {
+    paramFloor := 4
+    paramCeiling := 4096
+    paramSchedule := .custom #[4, 16, 31, 32, 33, 64, 256, 1024, 4096]
+    maxSecondsPerCall := 3.0
+    targetInnerNanos := 200000000
+    signalFloorMultiplier := 1.0
+    tags := #["division", "long", "cold"]
+  }
+
+/- Cost model: a cold Newton call includes reciprocal construction and three clipped
+Karatsuba products.  Doubling is geometric, so the balanced model remains
+`Θ(M(n))`, represented by the Karatsuba-range integer surrogate. -/
+setup_benchmark runNewtonDivision n => n * (Nat.sqrt n)
+  with prep := prepDivision
+  where {
+    paramFloor := 4
+    paramCeiling := 4096
+    paramSchedule := .custom #[4, 16, 31, 32, 33, 64, 256, 1024, 4096]
+    maxSecondsPerCall := 3.0
+    targetInnerNanos := 200000000
+    signalFloorMultiplier := 1.0
+    tags := #["division", "newton", "cold"]
+  }
+
+/- A balanced product tree performs one multiplication per internal node over
+geometrically growing degrees, for `O(M(n) log n)` total work. The integer
+surrogate retains the logarithmic level count without assuming a particular
+coefficient-kernel exponent. -/
+setup_benchmark runProductTree n => n * (Nat.log2 n + 1)
+  with prep := prepProductTree
+  where {
+    paramFloor := 4
+    paramCeiling := 16384
+    paramSchedule := .custom #[4, 16, 64, 256, 1024, 4096, 16384]
+    maxSecondsPerCall := 3.0
+    targetInnerNanos := 200000000
+    signalFloorMultiplier := 1.0
+    tags := #["product-tree", "karatsuba", "cold"]
+  }
+
+end Hex.PolyFastBench
+
+def main (args : List String) : IO UInt32 :=
+  LeanBench.Cli.dispatch args
diff --git a/conformance/HexPolyFast/Conformance.lean b/conformance/HexPolyFast/Conformance.lean
new file mode 100644
index 000000000..6c526a93d
--- /dev/null
+++ b/conformance/HexPolyFast/Conformance.lean
@@ -0,0 +1,124 @@
+/-
+Copyright (c) 2026 Lean FRO, LLC. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Kim Morrison
+-/
+
+module
+
+public import HexPolyFast
+public meta import HexPoly.Dense
+public meta import HexPolyFast.Plan
+public meta import HexPolyFast.Karatsuba
+public meta import HexPolyFast.Reverse
+public meta import HexPolyFast.Cyclic
+public meta import HexPolyFast.Reciprocal
+public meta import HexPolyFast.Division
+public meta import HexPolyFast.Tree
+
+public section
+
+/-!
+# Fast polynomial core conformance
+
+Oracle: the independent `DensePoly.mul` schoolbook implementation.
+Mode: always.
+
+The initial suite covers full and specialized-square Karatsuba recursion,
+cutoff zero, odd splits, normalized trailing zeros, strongly unbalanced
+blocking, arbitrary slices, and the guarded polynomial/series reversal bridge.
+-/
+
+namespace HexPolyFast.Conformance
+
+open Hex Hex.DensePoly
+
+private def a : DensePoly Int :=
+  ofList [3, -2, 0, 5, 1, -7, 4]
+
+private def b : DensePoly Int :=
+  ofList [-1, 6, 2, 0, -3]
+
+private def trailing : DensePoly Int :=
+  ofList [3, -2, 0, 5, 1, -7, 4, 0, 0]
+
+private def long : DensePoly Int :=
+  ofList ((List.range 65).map fun i => Int.ofNat i - 17)
+
+private def short : DensePoly Int := ofList [2, -3]
+
+#guard mulKaratsuba 0 a b = a * b
+#guard mulKaratsuba 2 a b = a * b
+#guard mulKaratsuba 3 trailing b = trailing * b
+#guard squareKaratsuba 0 a = a * a
+#guard squareKaratsuba 3 a = a * a
+#guard mulKaratsuba 2 long short = long * short
+#guard mulKaratsuba 2 short long = short * long
+
+private def plan : MulPlan Int := karatsubaPlan 2
+
+#guard mulWith plan a b = a * b
+#guard squareWith plan a = a * a
+#guard mulSlice plan 2 4 a b = schoolbookSlice 2 4 a b
+#guard mulSlice plan 50 4 a b = 0
+
+#guard (reverseSeries (C (1 : Int)) 2).coeffs.toArray.toList = [1, 0]
+#guard (polyOfSeries (reverseSeries a a.size)).coeff 0 = a.coeff 6
+#guard (polyOfSeries (reverseSeries a a.size)).coeff 3 = a.coeff 3
+#guard (polyOfSeries (reverseSeries a a.size)).coeff 6 = a.coeff 0
+#guard (polyOfSeries (reverseSeries a a.size)).coeff 7 = 0
+
+#guard mulCyclic (karatsubaPlan 2) 3 (by omega)
+    (ofList [1, 2, 3] : DensePoly Int) (ofList [4, 5]) = ofList [19, 13, 22]
+#guard mulNegacyclic (karatsubaPlan 2) 3 (by omega)
+    (ofList [1, 2, 3] : DensePoly Int) (ofList [4, 5]) = ofList [-11, 13, 22]
+#guard mulCyclic? (karatsubaPlan 2) 0 a b = none
+
+private def unitSeries : TSeries Int 8 :=
+  TSeries.ofFn fun i => [1, 2, -1, 3, 0, -2, 4, 1].getD i 0
+
+#guard seriesMulUpTo plan 5 unitSeries unitSeries =
+  TSeries.mulUpTo 5 unitSeries unitSeries
+#guard reciprocalWith plan unitSeries 1 = TSeries.invOfUnit unitSeries 1
+#guard unitSeries * reciprocalWith plan unitSeries 1 = 1
+
+private def monicDivisor : DensePoly Int := ofList [2, -3, 1]
+
+private theorem monicDivisor_monic : Monic monicDivisor := by
+  rfl
+
+#guard divModMonicWith plan a monicDivisor monicDivisor_monic =
+  divModMonic a monicDivisor monicDivisor_monic
+#guard divModMonicWith plan (monicDivisor * b) monicDivisor
+    monicDivisor_monic = (b, 0)
+#guard divModMonicWith plan (ofList [4, -1]) monicDivisor
+    monicDivisor_monic = (0, ofList [4, -1])
+#guard divModMonicWith plan a (C (1 : Int)) (by rfl) = (a, 0)
+
+private def ratA : DensePoly Rat :=
+  ofList [3, -2, 0, 5, 1, -7, 4]
+
+private def ratB : DensePoly Rat := ofList [2, -3, 5]
+
+#guard divModWith (karatsubaPlan 2) ratA ratB = divMod ratA ratB
+#guard divModWith (karatsubaPlan 2) ratA 0 = (0, ratA)
+#guard divModWith (karatsubaPlan 2) (ratA * ratB) ratB = (ratA, 0)
+
+private def treeLeaves : Array (DensePoly Int) :=
+  #[ofList [1, 1], ofList [2, 1], ofList [3, 1], ofList [4, 1], ofList [5, 1]]
+
+private def tree : ProductTree Int := ProductTree.build plan treeLeaves
+
+#guard tree.leaves = treeLeaves
+#guard tree.levelCount = 4
+#guard tree.root = treeLeaves.foldl (fun acc p => acc * p) 1
+#guard (tree.level? 0).map Array.size = some 5
+#guard (tree.level? 1).map Array.size = some 3
+#guard (tree.level? 2).map Array.size = some 2
+#guard (tree.level? 3).map Array.size = some 1
+#guard tree.nodeProduct? 1 0 = some (treeLeaves.getD 0 0 * treeLeaves.getD 1 0)
+#guard tree.nodeProduct? 1 2 = some (treeLeaves.getD 4 0)
+#guard tree.nodeProduct? 3 0 = some tree.root
+#guard tree.nodeProduct? 4 0 = none
+
+end HexPolyFast.Conformance
diff --git a/lakefile.lean b/lakefile.lean
index f382207e0..19b652898 100644
--- a/lakefile.lean
+++ b/lakefile.lean
@@ -90,6 +90,8 @@ lean_lib HexArith where
 
 lean_lib HexPoly where
 
+lean_lib HexPolyFast where
+
 lean_lib HexMvPoly where
 
 lean_lib HexMvGcd where
@@ -623,6 +625,8 @@ lean_lib HexConformance where
 
     ++ #[`HexTruncatedSeries.Conformance].map Glob.one
 
+    ++ #[`HexPolyFast.Conformance].map Glob.one
+
     ++ #[`HexMvHensel.Conformance, `HexMvFactor.Conformance].map Glob.one
 
     ++ #[`HexIntFactor.Conformance].map Glob.one
@@ -1063,6 +1067,10 @@ lean_exe hextruncatedseries_bench where
   srcDir := "bench"
   root := `HexTruncatedSeries.Bench
 
+lean_exe hexpolyfast_bench where
+  srcDir := "bench"
+  root := `HexPolyFast.Bench
+
 lean_exe hexpoly_emit_fixtures where
   srcDir := "conformance"
   root := `HexPoly.EmitFixtures
diff --git a/libraries.yml b/libraries.yml
index dcc69ed8b..0463f9651 100644
--- a/libraries.yml
+++ b/libraries.yml
@@ -126,6 +126,29 @@ libraries:
           description: Dense integer polynomials with nontrivial coefficient content for content and primitive-part operations.
         - name: polynomial-crt
           description: Coprime monic rational-polynomial moduli and residues for CRT witness construction.
+  HexPolyFast:
+    deps: [HexPoly, HexTruncatedSeries]
+    mathlib: false
+    done_through: 0
+    status: active
+    phase4:
+      comparators:
+        - tool: FLINT fmpz_poly and nmod_poly via python-flint
+          class: informational
+          rationale: "FLINT has independently tuned coefficient-specific dispatch; within-Lean agreement and crossover cells gate production selection."
+      input_families:
+        - name: full-and-clipped-multiplication
+          description: Balanced and unbalanced full, square, low, and middle products across crossover sizes.
+        - name: newton-division
+          description: One-shot and cached-divisor division against long division.
+        - name: half-gcd
+          description: Gcd, full xgcd, and one-sided xgcd across balanced and skewed degree pairs.
+        - name: multipoint
+          description: Cold and reused product/remainder trees for evaluation and interpolation.
+        - name: pade
+          description: Homogeneous and normalized Padé cases, including normalized failure.
+        - name: coefficient-kernels
+          description: Forced Kronecker and direct/CRT-NTT paths across degree, width, and modulus ladders.
   HexMvPoly:
     deps: [HexPoly, HexBasic]
     mathlib: false