perf(univariate): optimize eval and pow backends

CompPoly/Bivariate/ToPoly.lean

@@ -375,14 +375,16 @@ theorem evalEval_toPoly {R : Type*} [BEq R] [LawfulBEq R] [Nontrivial R] [Semiri
       have h_toPoly : (toPoly f).eval (Polynomial.C y) = (f.val.eval (CPolynomial.C y)).toPoly := by
         unfold CBivariate.toPoly
         simp +decide [ Polynomial.eval_finset_sum, CPolynomial.Raw.eval ]
-        unfold CPolynomial.Raw.eval₂
+        rw [CPolynomial.Raw.eval₂_eq_eval₂_naive]
+        unfold CPolynomial.Raw.eval₂Naive
         simp +decide
         rw [ toPoly_foldl_zipIdx_eq_sum, Finset.sum_subset ]
         · exact fun i hi ↦ Finset.mem_range.mpr
             (Nat.lt_of_lt_of_le (Finset.mem_range.mp (Finset.mem_filter.mp hi |>.1)) (by simp))
         · simp +contextual [ CPolynomial.support ]
           simp +decide [ CPolynomial.toPoly, CPolynomial.Raw.toPoly ]
-          unfold CPolynomial.Raw.eval₂
+          rw [CPolynomial.Raw.eval₂_eq_eval₂_naive]
+          unfold CPolynomial.Raw.eval₂Naive
           erw [ Array.foldl_empty ]
           simp
       -- `toPoly (f.val.eval (C y))` equals the polynomial with coefficients `f.val.coeff i`.
@@ -641,14 +643,16 @@ theorem evalY_toPoly {R : Type*} [BEq R] [LawfulBEq R] [Nontrivial R] [Semiring
   have h_toPoly : (toPoly f).eval (Polynomial.C a) = (f.val.eval (CPolynomial.C a)).toPoly := by
     unfold CBivariate.toPoly
     simp +decide [ Polynomial.eval_finset_sum, CPolynomial.Raw.eval ]
-    unfold CPolynomial.Raw.eval₂
+    rw [CPolynomial.Raw.eval₂_eq_eval₂_naive]
+    unfold CPolynomial.Raw.eval₂Naive
     simp +decide
     rw [ toPoly_foldl_zipIdx_eq_sum, Finset.sum_subset ]
     · exact fun i hi ↦ Finset.mem_range.mpr
         (Nat.lt_of_lt_of_le (Finset.mem_range.mp (Finset.mem_filter.mp hi |>.1)) (by simp))
     · simp +contextual [ CPolynomial.support ]
       simp +decide [ CPolynomial.toPoly, CPolynomial.Raw.toPoly ]
-      unfold CPolynomial.Raw.eval₂
+      rw [CPolynomial.Raw.eval₂_eq_eval₂_naive]
+      unfold CPolynomial.Raw.eval₂Naive
       erw [ Array.foldl_empty ]
       simp
   exact h_toPoly.symm

CompPoly/Univariate/Basic.lean

@@ -196,15 +196,16 @@ def natDegree [Zero R] (p : CPolynomial R) : ℕ :=
 of the trimmed array, or `0` if the trimmed array is empty. -/
 def leadingCoeff [Zero R] (p : CPolynomial R) : R := p.val.getLastD 0
 
-/-- Evaluate a polynomial at a point. -/
+/-- Evaluate a polynomial at a point using Horner's method. -/
 def eval [Semiring R] (x : R) (p : CPolynomial R) : R :=
-  p.val.zipIdx.foldl (fun acc ⟨a, i⟩ => acc + a * x ^ i) 0
+  p.val.eval x
 
 /-- Evaluate at `x : S` via a ring hom
-`f : R →+* S`; `eval₂ f x p = f(a₀) + f(a₁)*x + f(a₂)*x² + ...`. -/
+`f : R →+* S`; `eval₂ f x p = f(a₀) + f(a₁)*x + f(a₂)*x² + ...`.
+Uses the optimized Horner backend. -/
 def eval₂ {S : Type*} [Semiring R] [Semiring S]
     (f : R →+* S) (x : S) (p : CPolynomial R) : S :=
-  p.val.zipIdx.foldl (fun acc ⟨a, i⟩ => acc + f a * x ^ i) 0
+  p.val.eval₂ f x
 
 /-- The support of a polynomial: indices with nonzero coefficients. -/
 def support [Zero R] [BEq R] (p : CPolynomial R) : Finset ℕ :=
@@ -316,19 +317,13 @@ theorem support_empty_iff [Zero R] [BEq R] [LawfulBEq R] (p : CPolynomial R) :
 /-- Evaluation equals the sum over support of coefficients times powers. -/
 theorem eval_eq_sum_support [Semiring R] [BEq R] [LawfulBEq R] (p : CPolynomial R) (x : R) :
     p.eval x = p.support.sum (fun i => p.coeff i * x ^ i) := by
-  have h_eval_def : p.eval x =
-      (p.val.zipIdx.toList.map (fun ⟨a, i⟩ => a * x ^ i)).sum := by
-    unfold CPolynomial.eval
-    simp +decide
-    induction p.val
-    simp +decide [ * ]
-    induction' ‹List R› using List.reverseRecOn with a l ih <;>
-      simp +decide [ *, List.zipIdx_append ]
-  have h_sum_range : (p.val.zipIdx.toList.map (fun ⟨a, i⟩ => a * x ^ i)).sum =
+  have h_eval_sum : p.eval x =
       (Finset.range p.val.size).sum (fun i => p.val.coeff i * x ^ i) := by
-    convert CPolynomial.Raw.sum_zipIdx_eq_sum_range p.val (fun a i => a * x ^ i)
-      using 1
-  convert h_eval_def.trans h_sum_range using 1
+    show p.val.eval x = _
+    show p.val.eval₂ (RingHom.id R) x = _
+    rw [Raw.eval₂_eq_sum]
+    simp
+  convert h_eval_sum using 1
   refine' Finset.sum_subset _ _ <;> intro i hi <;>
     simp_all +decide [ CPolynomial.Raw.coeff ]
   · exact Finset.mem_range.mp (Finset.mem_filter.mp hi |>.1)
@@ -341,19 +336,11 @@ Evaluation via a ring hom equals the sum over support of mapped coefficients tim
 theorem eval₂_eq_sum_support {S : Type*} [Semiring R] [BEq R] [LawfulBEq R] [Semiring S]
     (f : R →+* S) (p : CPolynomial R) (x : S) :
     p.eval₂ f x = p.support.sum (fun i => f (p.coeff i) * x ^ i) := by
-  have h_eval_def : p.eval₂ f x =
-      (p.val.zipIdx.toList.map (fun ⟨a, i⟩ => f a * x ^ i)).sum := by
-    unfold CPolynomial.eval₂
-    simp +decide
-    induction p.val
-    simp +decide [*]
-    induction' ‹List R› using List.reverseRecOn with a l ih <;>
-      simp +decide [*, List.zipIdx_append]
-  have h_sum_range : (p.val.zipIdx.toList.map (fun ⟨a, i⟩ => f a * x ^ i)).sum =
+  have h_eval_sum : p.eval₂ f x =
       (Finset.range p.val.size).sum (fun i => f (p.val.coeff i) * x ^ i) := by
-    convert CPolynomial.Raw.sum_zipIdx_eq_sum_range p.val (fun a i => f a * x ^ i)
-      using 1
-  convert h_eval_def.trans h_sum_range using 1
+    show p.val.eval₂ f x = _
+    exact Raw.eval₂_eq_sum f x p.val
+  convert h_eval_sum using 1
   refine' Finset.sum_subset _ _ <;> intro i hi <;>
     simp_all +decide [CPolynomial.Raw.coeff]
   · exact Finset.mem_range.mp (Finset.mem_filter.mp hi |>.1)
@@ -736,25 +723,18 @@ lemma add_mul [Semiring R] [BEq R] [LawfulBEq R]
 
 lemma pow_is_trimmed [Semiring R] [BEq R] [LawfulBEq R] [Nontrivial R]
     (p : CPolynomial.Raw R) (n : ℕ) : (p ^ n).trim = p ^ n := by
-      induction' n with n ih generalizing p;
-      · convert one_is_trimmed
-        · infer_instance
-        · infer_instance
-      · have h_exp : p ^ (n + 1) = p * p ^ n := by
-          exact pow_succ p n
-        rw [h_exp]
-        convert mul_is_trimmed p ( p ^ n ) using 1
+      induction n with
+      | zero =>
+        show (Raw.pow p 0).trim = Raw.pow p 0
+        unfold Raw.pow
+        exact one_is_trimmed
+      | succ n ih =>
+        rw [Raw.pow_succ]
+        exact mul_is_trimmed p (p ^ n)
 
 lemma pow_succ_right [Semiring R] [BEq R] [LawfulBEq R] [Nontrivial R]
-    (p : CPolynomial.Raw R) (n : ℕ) : p ^ (n + 1) = p ^ n * p := by
-      convert pow_succ p n using 1;
-      induction' n with n ih;
-      · have h_pow_zero : p ^ 0 = 1 := by
-          exact rfl
-        rw [h_pow_zero, mul_one_trim, one_mul_trim];
-      · simp_all +decide [Raw.pow_succ];
-        convert Raw.mul_assoc p ( p ^ n ) p using 1;
-        grind
+    (p : CPolynomial.Raw R) (n : ℕ) : p ^ (n + 1) = p ^ n * p :=
+  Raw.pow_succ_right p n
 
 /--
 `CPolynomial R` forms a commutative monoid when `R` is a semiring.
@@ -795,16 +775,23 @@ instance [Semiring R] [BEq R] [LawfulBEq R] [Nontrivial R] : Semiring (CPolynomi
   nsmul_zero := CPolynomial.nsmul_zero
   nsmul_succ := CPolynomial.nsmul_succ
   npow n p := ⟨p.val ^ n, Trim.isCanonical_of_trim_eq (CPolynomial.pow_is_trimmed p.val n)⟩
-  npow_zero := by intro x; apply Subtype.ext; rfl
-  npow_succ := by intro n p; apply Subtype.ext; exact
-      (CPolynomial.pow_succ_right p.val n)
+  npow_zero := by
+    intro x
+    apply Subtype.ext
+    show Raw.pow x.val 0 = _
+    unfold Raw.pow
+    rfl
+  npow_succ := by
+    intro n p
+    apply Subtype.ext
+    exact CPolynomial.pow_succ_right p.val n
   natCast_zero := by rfl
   natCast_succ := by intro n; rfl
 
 /-- `C r * X^n = monomial n r` as canonical polynomials. -/
 lemma C_mul_X_pow_eq_monomial [Semiring R] [BEq R] [LawfulBEq R] [DecidableEq R] [Nontrivial R]
     (r : R) (n : ℕ) :
-    (C r : CPolynomial R) * (X ^ n) = monomial n r := by
+    (C r : CPolynomial R) * X ^ n = monomial n r := by
   by_cases hr : r = 0
   · convert Subtype.ext ?_
     convert zero_mul _
@@ -818,13 +805,13 @@ lemma C_mul_X_pow_eq_monomial [Semiring R] [BEq R] [LawfulBEq R] [DecidableEq R]
         (Raw.X : CPolynomial.Raw R) ^ n = 0 := by
       convert Raw.zero_mul _ using 1
       convert rfl
-      · exact Eq.symm ( Raw.trim_replicate_zero 1 )
+      · exact Eq.symm (Raw.trim_replicate_zero 1)
       · infer_instance
     convert h_lhs using 1
     exact Eq.symm (by induction n <;> simp +decide [*, Raw.monomial])
   · convert Subtype.ext ?_
     have h_trim : (Raw.mk #[r]).trim = Raw.C r := by
-      exact Trim.canonical_iff.mpr fun hp => hr
+      exact Trim.canonical_iff.mpr fun hp ↦ hr
     generalize_proofs at *
     convert Raw.C_mul_eq_smul_trim r (Raw.X ^ n) using 1
     · exact h_trim.symm ▸ rfl

CompPoly/Univariate/Quotient/Core.lean

@@ -407,18 +407,14 @@ lemma pow_descends [Semiring R] [BEq R] [LawfulBEq R] (n : ℕ) (p₁ p₂ : CPo
   unfold powDescending
   rw [Quotient.eq]
   simp [Raw.instSetoidCPolynomial]
-  unfold pow
-  have mul_pow_succ_equiv (p : CPolynomial.Raw R) (n : ℕ):
-    p.mul^[n + 1] (C 1) ≈ p.mul (p.mul^[n] (C 1)) := by
-    rw [mul_pow_succ]
+  rw [Raw.pow_eq_powIterate, Raw.pow_eq_powIterate]
   induction n with
-  | zero => simp
+  | zero => exact equiv_refl _
   | succ n ih =>
+    rw [powIterate_succ, powIterate_succ]
     calc
-      p₁.mul^[n + 1] (C 1) ≈ p₁.mul (p₁.mul^[n] (C 1)) := mul_pow_succ_equiv p₁ n
-      _ ≈ p₁.mul (p₂.mul^[n] (C 1)) := mul_equiv₂ p₁ _ _ ih
-      _ ≈ p₂.mul (p₂.mul^[n] (C 1)) := mul_equiv _ _ (p₂.mul^[n] (C 1)) heq
-      _ ≈ p₂.mul^[n + 1] (C 1) := equiv_symm (mul_pow_succ_equiv p₂ n)
+      p₁.mul (p₁.powIterate n) ≈ p₁.mul (p₂.powIterate n) := mul_equiv₂ p₁ _ _ ih
+      _ ≈ p₂.mul (p₂.powIterate n) := mul_equiv _ _ _ heq
 
 /-- Exponentiation on the quotient. -/
 @[inline, specialize]
@@ -631,6 +627,7 @@ lemma npow_zero : ∀ (x : QuotientCPolynomial R), x.pow 0 = 1 := by
   refine Quotient.inductionOn x ?_
   intro p; clear x
   apply Quotient.sound
+  show CPolynomial.Raw.pow p 0 ≈ C 1
   unfold CPolynomial.Raw.pow
   simp
 
@@ -648,10 +645,9 @@ lemma pow_succ_left (n : ℕ) (x : QuotientCPolynomial R) :
   refine Quotient.inductionOn x ?_
   intro p
   apply Quotient.sound
-  -- p.pow (n+1) = p * p.pow n is true by definition of pow for CPolynomial.Raw
-  -- By definition of pow, we have p.pow (n + 1) = p.mul (p.pow n).
-  have h_pow : p.pow (n + 1) = p.mul (p.pow n) := by
-    exact Function.iterate_succ_apply' _ _ _
+  -- By Raw.pow_succ, we have p ^ (n + 1) = p * p ^ n.
+  have h_pow : CPolynomial.Raw.pow p (n + 1) = p.mul (CPolynomial.Raw.pow p n) := by
+    exact Raw.pow_succ p n
   exact congrFun (congrArg coeff h_pow)
 
 /-
@@ -671,9 +667,9 @@ lemma npow_succ : ∀ (n : ℕ) (x : QuotientCPolynomial R), x.pow (n + 1) = x.p
   refine Quotient.inductionOn x ?_
   intro p; clear x
   apply Quotient.sound
-  -- By definition of exponentiation, we have `p.pow (n + 1) = p * p.pow n` for any `p`.
-  rw [show p.pow (n + 1) = p.mul (p.pow n) from by
-        exact Function.iterate_succ_apply' _ _ _]
+  -- By Raw.pow_succ, we have `p ^ (n + 1) = p * p ^ n`.
+  rw [show CPolynomial.Raw.pow p (n + 1) = p.mul (CPolynomial.Raw.pow p n) from
+        Raw.pow_succ p n]
   convert commute_pow_self n ( Quotient.mk ( Raw.instSetoidCPolynomial ) p ) using 1
   erw [ Quotient.eq ]
   rfl

CompPoly/Univariate/Raw/Ops.lean

@@ -24,13 +24,20 @@ section Semiring
 
 variable {S : Type*}
 
-/-- Evaluates a `CPolynomial.Raw` at `x : S` using a ring homomorphism `f : R →+* S`.
+/-- Naive sum-of-powers evaluation (reference implementation).
 
   Computes `f(a₀) + f(a₁) * x + f(a₂) * x² + ...` where `aᵢ` are the coefficients.
+  Retained as a specification target for the optimized Horner backend. -/
+def eval₂Naive [Semiring R] [Semiring S] (f : R →+* S) (x : S) (p : CPolynomial.Raw R) : S :=
+  p.zipIdx.foldl (fun acc ⟨a, i⟩ ↦ acc + f a * x ^ i) 0
+
+/-- Evaluates a `CPolynomial.Raw` at `x : S` using a ring homomorphism `f : R →+* S`.
 
-  TODO: define an efficient version of this with caching -/
+  Uses Horner's method, processing coefficients from high degree to low degree:
+  `f(a₀) + x * (f(a₁) + x * (... + x * f(aₙ)))`, which avoids explicit exponentiation. -/
+@[inline, specialize]
 def eval₂ [Semiring R] [Semiring S] (f : R →+* S) (x : S) (p : CPolynomial.Raw R) : S :=
-  p.zipIdx.foldl (fun acc ⟨a, i⟩ => acc + f a * x ^ i) 0
+  p.foldr (fun a acc ↦ f a + acc * x) 0
 
 /-- Evaluates a `CPolynomial.Raw` at a given value -/
 @[inline, specialize]
@@ -55,12 +62,12 @@ section SMulDefs
 /-- Scalar multiplication: multiplies each coefficient by `r`. -/
 @[inline, specialize]
 def smul [Mul R] (r : R) (p : CPolynomial.Raw R) : CPolynomial.Raw R :=
-  .mk (Array.map (fun a => r * a) p)
+  .mk (Array.map (fun a ↦ r * a) p)
 
 /-- Raw scalar multiplication by a natural number (may have trailing zeros). -/
 @[inline, specialize]
 def nsmulRaw [Semiring R] (n : ℕ) (p : CPolynomial.Raw R) : CPolynomial.Raw R :=
-  .mk (Array.map (fun a => n * a) p)
+  .mk (Array.map (fun a ↦ n * a) p)
 
 /-- Scalar multiplication of `CPolynomial.Raw` by a natural number, with result trimmed. -/
 @[inline, specialize]
@@ -85,21 +92,30 @@ end MulPowXDefs
 /-- Multiplication using the naive `O(n²)` algorithm: `Σᵢ (aᵢ * q) * X^i`. -/
 @[inline, specialize]
 def mul [Semiring R] [BEq R] (p q : CPolynomial.Raw R) : CPolynomial.Raw R :=
-  p.zipIdx.foldl (fun acc ⟨a, i⟩ => acc.add <| (smul a q).mulPowX i) (mk #[])
+  p.zipIdx.foldl (fun acc ⟨a, i⟩ ↦ acc.add <| (smul a q).mulPowX i) (mk #[])
 
-/-- Exponentiation of a `CPolynomial.Raw` by a natural number `n` via repeated multiplication. -/
-@[inline, specialize]
-def pow [Semiring R] [BEq R] (p : CPolynomial.Raw R) (n : Nat) : CPolynomial.Raw R :=
+/-- Linear exponentiation of a `CPolynomial.Raw` by repeated multiplication (reference impl). -/
+def powIterate [Semiring R] [BEq R] (p : CPolynomial.Raw R) (n : Nat) : CPolynomial.Raw R :=
   (mul p)^[n] (C 1)
 
+/-- Exponentiation of a `CPolynomial.Raw` by a natural number `n` via squaring. -/
+@[inline, specialize]
+def pow [Semiring R] [BEq R] (p : CPolynomial.Raw R) : Nat → CPolynomial.Raw R
+  | 0 => C 1
+  | 1 => p.mul (C 1)
+  | n + 2 =>
+    let half := pow p ((n + 2) / 2)
+    let sq := mul half half
+    if (n + 2) % 2 == 0 then sq else mul p sq
+
 instance : Zero (CPolynomial.Raw R) := ⟨#[]⟩
 instance [One R] : One (CPolynomial.Raw R) := ⟨C 1⟩
 instance [Zero R] [Add R] [BEq R] : Add (CPolynomial.Raw R) := ⟨add⟩
 instance [Mul R] : SMul R (CPolynomial.Raw R) := ⟨smul⟩
 instance [Semiring R] [BEq R] : SMul ℕ (CPolynomial.Raw R) := ⟨nsmul⟩
 instance [Semiring R] [BEq R] : Mul (CPolynomial.Raw R) := ⟨mul⟩
 instance [Semiring R] [BEq R] : Pow (CPolynomial.Raw R) Nat := ⟨pow⟩
-instance [NatCast R] : NatCast (CPolynomial.Raw R) := ⟨fun n => C (n : R)⟩
+instance [NatCast R] : NatCast (CPolynomial.Raw R) := ⟨fun n ↦ C (n : R)⟩
 
 /-- Upper bound on degree: `size - 1` if non-empty, `⊥` if empty. -/
 def degreeBound (p : CPolynomial.Raw R) : WithBot Nat :=
@@ -123,7 +139,7 @@ section Ring
 
 /-- Negation of a `CPolynomial.Raw`. -/
 @[inline, specialize]
-def neg [Neg R] (p : CPolynomial.Raw R) : CPolynomial.Raw R := p.map (fun a => -a)
+def neg [Neg R] (p : CPolynomial.Raw R) : CPolynomial.Raw R := p.map (fun a ↦ -a)
 
 /-- Subtraction of two `CPolynomial.Raw`s. -/
 @[inline, specialize]
@@ -133,7 +149,7 @@ def sub [Zero R] [Add R] [Neg R] [BEq R]
 
 instance [Neg R] : Neg (CPolynomial.Raw R) := ⟨neg⟩
 instance [Zero R] [Add R] [Neg R] [BEq R] : Sub (CPolynomial.Raw R) := ⟨sub⟩
-instance [IntCast R] : IntCast (CPolynomial.Raw R) := ⟨fun n => C (n : R)⟩
+instance [IntCast R] : IntCast (CPolynomial.Raw R) := ⟨fun n ↦ C (n : R)⟩
 
 /-- Erase the coefficient at index `n`: same as `p` except `coeff n = 0`, then trimmed. -/
 def erase [Zero R] [Add R] [Neg R] [BEq R] [DecidableEq R]