fix(poly-fast): prove remainder capacities

HexManual/Chapters/HexPolyFast.lean

@@ -233,8 +233,12 @@ their sibling subtree.
 
 {docstring Hex.DensePoly.RemainderTree.build}
 
+{docstring Hex.DensePoly.RemainderTree.rootDegree}
+
 {docstring Hex.DensePoly.RemainderTree.remainders?}
 
+{docstring Hex.DensePoly.RemainderTree.remainders?_isSome_of_capacity}
+
 {docstring Hex.DensePoly.RemainderTree.remainders?_sound}
 
 {name}`Hex.DensePoly.EvalPlan` specializes the leaves to `x - a` and caches

HexPolyFast/Multipoint.lean

@@ -97,25 +97,42 @@ inductive PointNode (R : Type u) [DecidableEq R]
       PointNode R mul (leftPoints ++ rightPoints)
         (mulWith mul leftPoly rightPoly)
 
+/-- The cached child plans have the exact divisor sizes and sibling capacities
+needed by a bounded remainder-tree traversal. -/
+def PointNode.CapacitySafe {mul : MulPlan R} :
+    {points : List R} → {poly : DensePoly R} → PointNode R mul points poly → Prop
+  | _, _, .leaf _ => True
+  | _, _, @PointNode.branch _ _ _ _ leftPoints rightPoints _ _
+      left right leftPlan rightPlan _ _ =>
+      left.CapacitySafe ∧ right.CapacitySafe ∧
+        leftPlan.divisor.size = leftPoints.length + 1 ∧
+        rightPlan.divisor.size = rightPoints.length + 1 ∧
+        leftPlan.capacity = rightPoints.length ∧
+        rightPlan.capacity = leftPoints.length
+
 /-- A constructed node together with its monic product polynomial. -/
 private structure BuiltNode (mul : MulPlan R) (points : List R) where
   poly : DensePoly R
   node : PointNode R mul points poly
   monic : poly.Monic
   ne : poly ≠ 0
+  size_eq : poly.size = points.length + 1
   evalZero : ∀ a, a ∈ points → poly.eval a = 0
+  capacitySafe : node.CapacitySafe
 
 private def leafNode (mul : MulPlan R) (hone : (1 : R) ≠ 0) (a : R) :
     BuiltNode mul [a] :=
   { poly := pointFactor a
     node := .leaf a
     monic := pointFactor_monic a
     ne := (pointFactor_monic a).neOfOneNe hone
+    size_eq := pointFactor_size a hone
     evalZero := by
       intro x hx
       simp only [List.mem_singleton] at hx
       subst x
-      exact pointFactor_eval a }
+      exact pointFactor_eval a
+    capacitySafe := True.intro }
 
 private def branchNode (mul : MulPlan R) (hone : (1 : R) ≠ 0)
     {leftPoints rightPoints : List R}
@@ -124,6 +141,14 @@ private def branchNode (mul : MulPlan R) (hone : (1 : R) ≠ 0)
   let leftPlan := DivPlan.ofMonic mul left.poly left.monic left.ne rightPoints.length
   let rightPlan := DivPlan.ofMonic mul right.poly right.monic right.ne leftPoints.length
   let poly := mulWith mul left.poly right.poly
+  have hleftPos : 0 < left.poly.size := by rw [left.size_eq]; omega
+  have hrightPos : 0 < right.poly.size := by rw [right.size_eq]; omega
+  have hproduct := size_mul_of_top_ne left.poly right.poly hleftPos hrightPos (by
+    rw [leadingCoeff_eq_one_of_monic left.monic,
+      leadingCoeff_eq_one_of_monic right.monic]
+    have hzero : (Zero.zero : R) = 0 := rfl
+    rw [hzero, Lean.Grind.Semiring.one_mul]
+    exact hone)
   { poly
     node := .branch left.node right.node leftPlan rightPlan
       (by intro a ha; simpa [leftPlan] using left.evalZero a ha)
@@ -138,14 +163,29 @@ private def branchNode (mul : MulPlan R) (hone : (1 : R) ≠ 0)
         rw [mulWith_eq]
         exact left.monic.mul right.monic
       · exact hone
+    size_eq := by
+      dsimp [poly]
+      rw [mulWith_eq, hproduct, left.size_eq, right.size_eq, List.length_append]
+      omega
     evalZero := by
       intro a ha
       dsimp [poly]
       rw [mulWith_eq, eval_mul_commring]
       rw [List.mem_append] at ha
       cases ha with
       | inl hleft => rw [left.evalZero a hleft]; grind
-      | inr hright => rw [right.evalZero a hright]; grind }
+      | inr hright => rw [right.evalZero a hright]; grind
+    capacitySafe := by
+      change left.node.CapacitySafe ∧ right.node.CapacitySafe ∧
+        leftPlan.divisor.size = leftPoints.length + 1 ∧
+        rightPlan.divisor.size = rightPoints.length + 1 ∧
+        leftPlan.capacity = rightPoints.length ∧
+        rightPlan.capacity = leftPoints.length
+      exact ⟨left.capacitySafe, right.capacitySafe,
+        (by simp [leftPlan, left.size_eq]),
+        (by simp [rightPlan, right.size_eq]),
+        (by simp [leftPlan]),
+        (by simp [rightPlan])⟩ }
 
 /-- Recursively build a count-balanced cached remainder tree. -/
 private def buildNode (mul : MulPlan R) (hone : (1 : R) ≠ 0) :
@@ -194,7 +234,8 @@ structure EvalPlan (R : Type u) [DecidableEq R] [Lean.Grind.CommRing R] where
   /-- Empty and trivial-ring plans use `none`; otherwise the type index ties
   every cached divisor to the exact public point sequence. -/
   private nodeData : Option
-    (Sigma fun poly => PointNode R mulData pointsData.toList poly)
+    (Sigma fun poly => { node : PointNode R mulData pointsData.toList poly //
+      node.CapacitySafe })
 
 namespace EvalPlan
 
@@ -206,7 +247,7 @@ def build (mul : MulPlan R) (points : Array R) : EvalPlan R :=
       else if hempty : points.toList = [] then none
       else
         let built := buildNode mul hone points.toList hempty
-        some ⟨built.poly, built.node⟩ }
+        some ⟨built.poly, ⟨built.node, built.capacitySafe⟩⟩ }
 
 /-- The planned point sequence. -/
 def points (plan : EvalPlan R) : Array R := plan.pointsData
@@ -226,14 +267,13 @@ def mulPlan (plan : EvalPlan R) : MulPlan R := plan.mulData
 nondegenerate. -/
 def cachedNode (plan : EvalPlan R) : Option
     (Sigma fun poly => PointNode R plan.mulPlan plan.points.toList poly) :=
-  plan.nodeData
+  plan.nodeData.map fun built => ⟨built.1, built.2.1⟩
 
 /-- A nonempty point sequence over a nontrivial ring has a cached node. -/
 theorem cachedNode_build_isSome (mul : MulPlan R) (points : Array R)
     (hone : (1 : R) ≠ 0) (hne : points.toList ≠ []) :
     (build mul points).cachedNode.isSome := by
   simp [cachedNode, mulPlan, build, hone, hne]
-  rfl
 
 /-- Rebuild the observational product-tree view. Constructing an evaluation
 plan does not eagerly build this redundant level representation. -/
@@ -260,64 +300,99 @@ theorem eval_eq_map (plan : EvalPlan R) (f : DensePoly R) :
     plan.eval f = plan.points.map (f.eval ·) := by
   rfl
 
-/-- Reduce one parent remainder into a cached child node.  The capacity guard
-is executable defensive checking; plans built by `build` satisfy it throughout
-the bounded traversal. -/
-private def reduceNode (parent : DensePoly R) (node : DivPlan R) : DensePoly R :=
-  if hcap : quotientLength parent node.divisor ≤ node.capacity then
-    node.mod parent hcap
-  else
-    parent
+/-- A divisor of size `degree + 1` needs at most `capacity` quotient
+coefficients when the parent has size at most `degree + capacity`. -/
+private theorem quotientLength_le_capacity (parent : DensePoly R)
+    (node : DivPlan R) {degree capacity : Nat}
+    (hdivisor : node.divisor.size = degree + 1)
+    (hparent : parent.size ≤ degree + capacity) :
+    quotientLength parent node.divisor ≤ capacity := by
+  have hne : node.divisor.size ≠ 0 := by rw [hdivisor]; omega
+  by_cases hlt : parent.size < node.divisor.size
+  · simp [quotientLength_eq, hlt]
+  · rw [quotientLength_eq]
+    simp [hne, hlt]
+    rw [hdivisor] at hlt
+    omega
 
 /-- Reduction by a divisor vanishing at `a` preserves evaluation at `a`. -/
-private theorem eval_reduceNode (parent : DensePoly R) (node : DivPlan R) (a : R)
+private theorem eval_mod (parent : DensePoly R) (node : DivPlan R)
+    (hcap : quotientLength parent node.divisor ≤ node.capacity) (a : R)
     (hzero : node.divisor.eval a = 0) :
-    (reduceNode parent node).eval a = parent.eval a := by
-  unfold reduceNode
-  split
-  · rename_i hcap
-    rw [node.mod_eq parent hcap]
-    simp only [eval_sub_ring, mulWith_eq, eval_mul_commring, hzero,
-      Lean.Grind.Semiring.mul_zero]
-    grind
-  · rfl
+    (node.mod parent hcap).eval a = parent.eval a := by
+  rw [node.mod_eq parent hcap]
+  simp only [eval_sub_ring, mulWith_eq, eval_mul_commring, hzero,
+    Lean.Grind.Semiring.mul_zero]
+  grind
+
+private theorem map_eval_mod (parent : DensePoly R) (node : DivPlan R)
+    (hcap : quotientLength parent node.divisor ≤ node.capacity)
+    (points : List R) (hzero : ∀ a, a ∈ points → node.divisor.eval a = 0) :
+    points.map ((node.mod parent hcap).eval ·) =
+      points.map (parent.eval ·) := by
+  apply List.map_congr_left
+  intro a ha
+  exact eval_mod parent node hcap a (hzero a ha)
 
 /-- Execute the balanced remainder tree, preserving left-to-right point order. -/
 private def evalNode {points : List R} {poly : DensePoly R} :
-    PointNode R plan points poly → DensePoly R → List R
-  | .leaf a, f => [f.eval a]
-  | .branch left right leftPlan rightPlan _ _, f =>
-      evalNode left (reduceNode f leftPlan) ++
-        evalNode right (reduceNode f rightPlan)
+    (node : PointNode R plan points poly) →
+      node.CapacitySafe → (f : DensePoly R) → f.size ≤ points.length → List R
+  | .leaf a, _, f, _ => [f.eval a]
+  | @PointNode.branch _ _ _ _ leftPoints rightPoints _ _
+      left right leftPlan rightPlan leftZero rightZero, safe, f, hsize =>
+      have leftSafe := safe.1
+      have rightSafe := safe.2.1
+      have leftSize := safe.2.2.1
+      have rightSize := safe.2.2.2.1
+      have leftCapacity := safe.2.2.2.2.1
+      have rightCapacity := safe.2.2.2.2.2
+      have hleftCap : quotientLength f leftPlan.divisor ≤ leftPlan.capacity := by
+        rw [leftCapacity]
+        exact quotientLength_le_capacity f leftPlan leftSize (by
+          simp [List.length_append] at hsize ⊢
+          omega)
+      have hrightCap : quotientLength f rightPlan.divisor ≤ rightPlan.capacity := by
+        rw [rightCapacity]
+        exact quotientLength_le_capacity f rightPlan rightSize (by
+          simp [List.length_append] at hsize ⊢
+          omega)
+      let leftRemainder := leftPlan.mod f hleftCap
+      let rightRemainder := rightPlan.mod f hrightCap
+      have hleftSize : leftRemainder.size ≤ leftPoints.length := by
+        have h := leftPlan.size_mod_le f hleftCap
+        rw [leftSize] at h
+        simp [leftRemainder] at h ⊢
+        omega
+      have hrightSize : rightRemainder.size ≤ rightPoints.length := by
+        have h := rightPlan.size_mod_le f hrightCap
+        rw [rightSize] at h
+        simp [rightRemainder] at h ⊢
+        omega
+      evalNode left leftSafe leftRemainder hleftSize ++
+        evalNode right rightSafe rightRemainder hrightSize
 
 private theorem evalNode_eq {points : List R} {poly : DensePoly R}
-    (node : PointNode R plan points poly)
-    (f : DensePoly R) : evalNode node f = points.map (f.eval ·) := by
+    {node : PointNode R plan points poly} (safe : node.CapacitySafe)
+    (f : DensePoly R) (hsize : f.size ≤ points.length) :
+    evalNode node safe f hsize = points.map (f.eval ·) := by
   induction node generalizing f with
   | leaf a => rfl
   | @branch leftPoints rightPoints leftPoly rightPoly left right leftPlan rightPlan
       leftZero rightZero leftIH rightIH =>
       simp only [evalNode, List.map_append]
       rw [leftIH, rightIH]
-      have hleft : leftPoints.map ((reduceNode f leftPlan).eval ·) =
-          leftPoints.map (f.eval ·) := by
-        apply List.map_congr_left
-        intro a ha
-        exact eval_reduceNode f leftPlan a (leftZero a ha)
-      have hright : rightPoints.map ((reduceNode f rightPlan).eval ·) =
-          rightPoints.map (f.eval ·) := by
-        apply List.map_congr_left
-        intro a ha
-        exact eval_reduceNode f rightPlan a (rightZero a ha)
-      rw [hleft, hright]
+      rw [map_eval_mod (hzero := leftZero), map_eval_mod (hzero := rightZero)]
 
 /-- Executable cached remainder-tree evaluation.  Oversized inputs and the
 degenerate ring use direct Horner evaluation, exactly as specified. -/
 def evalImpl (plan : EvalPlan R) (f : DensePoly R) : Array R :=
-  if f.size ≤ plan.size then
+  if hsize : f.size ≤ plan.size then
     match plan.nodeData with
     | none => plan.points.map (f.eval ·)
-    | some built => (evalNode built.2 f).toArray
+    | some built =>
+        (evalNode built.2.1 built.2.2 f (by
+          simpa [size, points] using hsize)).toArray
   else
     plan.points.map (f.eval ·)
 
@@ -331,7 +406,7 @@ theorem eval_eq_impl (plan : EvalPlan R) (f : DensePoly R) :
     | none => rfl
     | some built =>
         simp only
-        rw [evalNode_eq built.2 f]
+        rw [evalNode_eq built.2.2 f]
         unfold points
         rw [← Array.toList_map, Array.toArray_toList]
   · rfl