diff --git a/CompPoly/Bivariate/ToPoly.lean b/CompPoly/Bivariate/ToPoly.lean index 0c4a86d8..469a84c6 100644 --- a/CompPoly/Bivariate/ToPoly.lean +++ b/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 diff --git a/CompPoly/Univariate/Basic.lean b/CompPoly/Univariate/Basic.lean index c6ecb5a9..90a65ab7 100644 --- a/CompPoly/Univariate/Basic.lean +++ b/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 diff --git a/CompPoly/Univariate/Quotient/Core.lean b/CompPoly/Univariate/Quotient/Core.lean index 9a35d2dc..c407ca7e 100644 --- a/CompPoly/Univariate/Quotient/Core.lean +++ b/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 diff --git a/CompPoly/Univariate/Raw/Ops.lean b/CompPoly/Univariate/Raw/Ops.lean index a28f7719..09276c29 100644 --- a/CompPoly/Univariate/Raw/Ops.lean +++ b/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,13 +92,22 @@ 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⟩ @@ -99,7 +115,7 @@ 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] diff --git a/CompPoly/Univariate/Raw/Proofs.lean b/CompPoly/Univariate/Raw/Proofs.lean index c5abdfc8..4f0bedcb 100644 --- a/CompPoly/Univariate/Raw/Proofs.lean +++ b/CompPoly/Univariate/Raw/Proofs.lean @@ -30,12 +30,98 @@ variable (p q r : CPolynomial.Raw R) lemma pow_zero (p : CPolynomial.Raw R) : p ^ 0 = C 1 := by - exact rfl + show pow p 0 = C 1 + unfold pow + rfl -lemma pow_succ (p : CPolynomial.Raw R) (n : ℕ) : - p ^ (n + 1) = p * (p ^ n) := by - convert ( Function.iterate_succ_apply' ( mul p ) n ( C 1 ) ) - using 1 +/-- `powIterate` unfolds one step: `powIterate p (n+1) = p * powIterate p n`. -/ +lemma powIterate_succ (p : CPolynomial.Raw R) (n : ℕ) : + powIterate p (n + 1) = mul p (powIterate p n) := + Function.iterate_succ_apply' (mul p) n (C 1) + +end Semiring + +section EvalBridge + +variable {R : Type*} [Semiring R] +variable {S : Type*} [Semiring S] + +/-- Horner eval₂ on an empty array is zero. -/ +@[simp] +lemma eval₂_empty (f : R →+* S) (x : S) : + eval₂ f x (#[] : CPolynomial.Raw R) = 0 := by + simp [eval₂] + +/-- Naive eval₂ on an empty array is zero. -/ +@[simp] +lemma eval₂Naive_empty (f : R →+* S) (x : S) : + eval₂Naive f x (#[] : CPolynomial.Raw R) = 0 := by + simp [eval₂Naive] + +/- +Helper: Horner foldr equals the naive sum-of-powers for lists. +-/ +private lemma horner_eq_naive_list (f : R →+* S) (x : S) : + ∀ (l : List R), + l.foldr (fun a acc ↦ f a + acc * x) 0 = + (l.zipIdx.map (fun ⟨a, i⟩ ↦ f a * x ^ i)).sum := by + intro l + induction' l using List.reverseRecOn with a l ih + · rfl + · simp +decide [*, List.zipIdx_append] + rw [← ih] + clear ih + induction a <;> simp +decide [*, pow_succ, mul_assoc, add_mul, add_assoc] + +/- +The Horner backend agrees with the naive sum-of-powers backend. +-/ +theorem eval₂_eq_eval₂_naive (f : R →+* S) (x : S) (p : CPolynomial.Raw R) : + eval₂ f x p = eval₂Naive f x p := by + convert horner_eq_naive_list f x p.toList using 1 + · unfold eval₂ + aesop + · unfold CPolynomial.Raw.eval₂Naive + induction p using Array.recOn + simp +decide [*] + induction ‹List R› using List.reverseRecOn <;> + simp +decide [*, List.zipIdx_append] + +/- +`eval₂` equals the Finset sum over `range p.size`. +-/ +theorem eval₂_eq_sum (f : R →+* S) (x : S) (p : CPolynomial.Raw R) : + eval₂ f x p = + (Finset.range p.size).sum (fun i ↦ f (p.coeff i) * x ^ i) := by + convert eval₂_eq_eval₂_naive f x p using 1 + · unfold CPolynomial.Raw.eval₂Naive + induction p using Array.recOn + simp_all +decide + induction ‹List R› using List.reverseRecOn <;> + simp_all +decide [Finset.sum_range_succ] + simp_all +decide [Finset.sum_range, List.zipIdx_append] + +/- +`eval₂` equals the list sum of mapped `zipIdx` pairs. +-/ +theorem eval₂_eq_list_sum (f : R →+* S) (x : S) (p : CPolynomial.Raw R) : + eval₂ f x p = + (p.zipIdx.toList.map (fun ⟨a, i⟩ ↦ f a * x ^ i)).sum := by + convert eval₂_eq_eval₂_naive f x p using 1 + · unfold CPolynomial.Raw.eval₂Naive + induction p using Array.recOn + simp +decide [*] + induction' ‹List R› using List.reverseRecOn with _ _ ih <;> + simp +decide [*, List.zipIdx_append] + +end EvalBridge + +section Semiring + +variable {R : Type*} [Semiring R] [BEq R] +variable {Q : Type*} [Semiring Q] +variable {S : Type*} +variable (p q r : CPolynomial.Raw R) section AddDefs @@ -659,18 +745,6 @@ lemma mulX_monomial_one [DecidableEq R] [LawfulBEq R] [Nontrivial R] (n : ℕ) : exact Nat.recOn n (by simp +decide) fun n ih => by simp +decide [ List.replicate ] at ih ⊢; tauto -lemma X_pow_eq_monomial_one [DecidableEq R] [LawfulBEq R] [Nontrivial R] (n : ℕ) : - (X : CPolynomial.Raw R) ^ n = monomial n 1 := by - have h_monomial : ∀ n : ℕ, - (monomial n (1 : R)).trim = - monomial n (1 : R) := by - exact fun n => monomial_canonical n 1 - induction' n with n ih; - · unfold X monomial - aesop - · rw [ pow_succ, ih ]; - rw [ X_mul_eq_mulX_trim ]; - rw [ mulX_monomial_one, h_monomial ] lemma smul_monomial_one_trim [DecidableEq R] [LawfulBEq R] [Nontrivial R] (n : ℕ) (r : R) : @@ -905,6 +979,103 @@ protected theorem mul_assoc [LawfulBEq R] (p q r : CPolynomial.Raw R) : end MulTheorems +section PowTheorems + +variable [LawfulBEq R] + +/-- Multiplication absorbs a left `trim`: `p.trim * q = p * q`. -/ +lemma mul_trim_left (p q : CPolynomial.Raw R) : mul p.trim q = mul p q := by + apply Trim.canonical_ext (mul_is_trimmed _ _) (mul_is_trimmed _ _) + intro i; rw [mul_coeff, mul_coeff]; congr 1; ext j; congr 1 + exact Trim.coeff_eq_coeff p j + +/-- Multiplication absorbs a right `trim`: `p * q.trim = p * q`. -/ +lemma mul_trim_right (p q : CPolynomial.Raw R) : mul p q.trim = mul p q := by + apply Trim.canonical_ext (mul_is_trimmed _ _) (mul_is_trimmed _ _) + intro i; rw [mul_coeff, mul_coeff]; congr 1; ext j; congr 1 + exact Trim.coeff_eq_coeff q _ + +/-- Product of iterate-powers adds exponents (for positive first exponent). -/ +lemma powIterate_mul_add (p : CPolynomial.Raw R) (m n : ℕ) : + mul (powIterate p (m + 1)) (powIterate p n) = powIterate p (m + 1 + n) := by + induction m with + | zero => + simp only [Nat.zero_add] + show mul (powIterate p 1) (powIterate p n) = powIterate p (1 + n) + have h1 : powIterate p 1 = p.trim := by + simp [powIterate_succ, show powIterate p 0 = C 1 from rfl] + exact mul_one_trim p + rw [h1, mul_trim_left, show 1 + n = n + 1 from by omega, powIterate_succ] + | succ m ih => + have h_unfold : powIterate p (m + 1 + 1) = mul p (powIterate p (m + 1)) := by + exact powIterate_succ p (m + 1) + rw [h_unfold] + have h_assoc : mul (mul p (powIterate p (m + 1))) (powIterate p n) = + mul p (mul (powIterate p (m + 1)) (powIterate p n)) := by + exact Raw.mul_assoc p (powIterate p (m + 1)) (powIterate p n) + rw [h_assoc, ih, ← powIterate_succ] + congr 1; omega + +/-- The fast squaring backend agrees with linear iteration under `LawfulBEq`. -/ +theorem pow_eq_powIterate (p : CPolynomial.Raw R) : ∀ n : ℕ, + pow p n = powIterate p n + | 0 => by unfold pow; rfl + | 1 => by unfold pow; rfl + | n + 2 => by + unfold pow + have ih := pow_eq_powIterate p ((n + 2) / 2) + have h_half : (n + 2) / 2 = n / 2 + 1 := by omega + rw [ih, h_half] + by_cases h : (n + 2) % 2 = 0 + · simp [h] + rw [powIterate_mul_add] + have : n / 2 + 1 + (n / 2 + 1) = n + 2 := by omega + rw [this] + · simp [show (n + 2) % 2 = 1 from by omega] + rw [powIterate_mul_add] + have h_sum : n / 2 + 1 + (n / 2 + 1) = n + 1 := by omega + rw [h_sum, ← powIterate_succ] + +/-- `p ^ (n + 1) = p * p ^ n` under `LawfulBEq`. -/ +lemma pow_succ (p : CPolynomial.Raw R) (n : ℕ) : + p ^ (n + 1) = p * (p ^ n) := by + change pow p (n + 1) = mul p (pow p n) + rw [pow_eq_powIterate, pow_eq_powIterate] + exact powIterate_succ p n + +/-- `p` commutes with `p ^ n` under multiplication (self-commutativity). -/ +lemma pow_mul_comm (p : CPolynomial.Raw R) : ∀ n : ℕ, + p * (p ^ n) = p ^ n * p + | 0 => by + show mul p (pow p 0) = mul (pow p 0) p + have h0 : pow p 0 = C 1 := by unfold pow; rfl + rw [h0] + show p * 1 = 1 * p + rw [mul_one_trim, one_mul_trim] + | n + 1 => by + have ih := pow_mul_comm p n + rw [Raw.pow_succ] + calc p * (p * p ^ n) = p * (p ^ n * p) := by rw [ih] + _ = (p * p ^ n) * p := by rw [Raw.mul_assoc] + +/-- `p ^ (n + 1) = p ^ n * p` under `LawfulBEq`. -/ +lemma pow_succ_right (p : CPolynomial.Raw R) (n : ℕ) : + p ^ (n + 1) = p ^ n * p := by + rw [pow_succ] + exact pow_mul_comm p n + +lemma X_pow_eq_monomial_one [DecidableEq R] [Nontrivial R] (n : ℕ) : + (X : CPolynomial.Raw R) ^ n = monomial n 1 := by + have h_monomial : ∀ n : ℕ, + (monomial n (1 : R)).trim = + monomial n (1 : R) := by + exact fun n => monomial_canonical n 1 + induction' n with n ih; + · unfold X monomial; show pow _ 0 = _; unfold pow; aesop + · rw [ pow_succ X n, ih, X_mul_eq_mulX_trim, mulX_monomial_one, h_monomial ] + +end PowTheorems + end Semiring section CommutativeSemiring diff --git a/CompPoly/Univariate/ToPoly/Core.lean b/CompPoly/Univariate/ToPoly/Core.lean index 11edab61..daf4d04e 100644 --- a/CompPoly/Univariate/ToPoly/Core.lean +++ b/CompPoly/Univariate/ToPoly/Core.lean @@ -67,7 +67,9 @@ alias ofPoly := Polynomial.toImpl /-- Evaluation is preserved by `toPoly`. -/ theorem eval_toPoly_eq_eval (x : Q) (p : CPolynomial.Raw Q) : p.toPoly.eval x = p.eval x := by - unfold Raw.toPoly Raw.eval eval₂ + unfold Raw.toPoly Raw.eval + rw [eval₂_eq_eval₂_naive, eval₂_eq_eval₂_naive] + unfold eval₂Naive rw [← Array.foldl_hom (Polynomial.eval x) (g₁ := fun acc (t : Q × ℕ) ↦ acc + Polynomial.C t.1 * Polynomial.X ^ t.2) (g₂ := fun acc (a, i) ↦ acc + a * x ^ i) ] @@ -76,7 +78,9 @@ theorem eval_toPoly_eq_eval (x : Q) (p : CPolynomial.Raw Q) : p.toPoly.eval x = /-- The coefficients of `p.toPoly` match those of `p`. -/ lemma coeff_toPoly {p : CPolynomial.Raw Q} {n : ℕ} : p.toPoly.coeff n = p.coeff n := by - unfold Raw.toPoly Raw.eval₂ + unfold Raw.toPoly + rw [eval₂_eq_eval₂_naive] + unfold eval₂Naive let f := fun (acc: Q[X]) ((a,i): Q × ℕ) ↦ acc + Polynomial.C a * Polynomial.X ^ i change (Array.foldl f 0 p.zipIdx).coeff n = p.coeff n diff --git a/CompPoly/Univariate/ToPoly/Impl.lean b/CompPoly/Univariate/ToPoly/Impl.lean index 341a0019..c755f713 100644 --- a/CompPoly/Univariate/ToPoly/Impl.lean +++ b/CompPoly/Univariate/ToPoly/Impl.lean @@ -54,10 +54,12 @@ theorem eval_toPoly [BEq R] [LawfulBEq R] (x : R) (p : CPolynomial R) : theorem Raw.eval₂_toPoly {S : Type*} [Semiring S] (f : R →+* S) (x : S) (p : CPolynomial.Raw R) : p.eval₂ f x = p.toPoly.eval₂ f x := by - unfold CompPoly.CPolynomial.Raw.toPoly CompPoly.CPolynomial.Raw.eval₂ - rw [← Array.foldl_hom (fun q : R[X] => q.eval₂ f x) - (g₁ := fun acc (t : R × ℕ) => acc + Polynomial.C t.1 * Polynomial.X ^ t.2) - (g₂ := fun acc (a, i) => acc + f a * x ^ i)] + unfold CompPoly.CPolynomial.Raw.toPoly + rw [CPolynomial.Raw.eval₂_eq_eval₂_naive, CPolynomial.Raw.eval₂_eq_eval₂_naive] + unfold CompPoly.CPolynomial.Raw.eval₂Naive + rw [← Array.foldl_hom (fun q : R[X] ↦ q.eval₂ f x) + (g₁ := fun acc (t : R × ℕ) ↦ acc + Polynomial.C t.1 * Polynomial.X ^ t.2) + (g₂ := fun acc (a, i) ↦ acc + f a * x ^ i)] · simp · intro acc t rcases t with ⟨a, i⟩ @@ -67,9 +69,7 @@ theorem Raw.eval₂_toPoly {S : Type*} [Semiring S] theorem eval₂_toPoly {S : Type*} [Semiring S] (f : R →+* S) (x : S) (p : CPolynomial R) : eval₂ f x p = p.toPoly.eval₂ f x := by - simpa [CompPoly.CPolynomial.eval₂, CompPoly.CPolynomial.toPoly, CompPoly.CPolynomial.Raw.eval₂] - using - (Raw.eval₂_toPoly (f := f) (x := x) (p := p.val)) + exact Raw.eval₂_toPoly f x p.val /-- CPolynomial.coeff is correct wrt the Mathlib spec. -/ theorem coeff_toPoly [BEq R] [LawfulBEq R] (p : CPolynomial R) (i : ℕ) : @@ -187,11 +187,12 @@ theorem erase_toPoly {R : Type*} [Ring R] [BEq R] [LawfulBEq R] [DecidableEq R] intros p q; have h_erase_toPoly : ∀ (p q : CPolynomial.Raw R), (p + -q).toPoly = p.toPoly + (-q).toPoly := by - exact fun p q => Raw.toPoly_add p (-q) + exact fun p q ↦ Raw.toPoly_add p (-q) convert h_erase_toPoly p q using 1 - simp +decide [ Raw.toPoly ] - rw [ show ( -q : CompPoly.CPolynomial.Raw R ) = q.map ( fun x => -x ) from ?_ ] - · simp +decide [ Raw.eval₂ ] + simp +decide [Raw.toPoly] + rw [show (-q : CompPoly.CPolynomial.Raw R) = q.map (fun x ↦ -x) from ?_] + · simp only [CPolynomial.Raw.eval₂_eq_eval₂_naive] + simp +decide [ Raw.eval₂Naive ] induction' q using Array.recOn with q ih; simp +decide [ *, Array.zipIdx ] induction' q using List.reverseRecOn with q ih <;> simp +decide [ *, List.mapIdx_append ] diff --git a/tests/CompPolyTests.lean b/tests/CompPolyTests.lean index e2185125..9e87ed16 100644 --- a/tests/CompPolyTests.lean +++ b/tests/CompPolyTests.lean @@ -1,3 +1,4 @@ +import CompPolyTests.Univariate.Eval import CompPolyTests.Univariate.Raw import CompPolyTests.Univariate.Basic import CompPolyTests.Univariate.Linear diff --git a/tests/CompPolyTests/Univariate/Basic.lean b/tests/CompPolyTests/Univariate/Basic.lean index 142c70e1..2f59e9f1 100644 --- a/tests/CompPolyTests/Univariate/Basic.lean +++ b/tests/CompPolyTests/Univariate/Basic.lean @@ -41,4 +41,41 @@ example (p : CPolynomial ℚ) : 1 * p = p := by simp end CPolynomial + +-- Power regression tests: canonical powers agree with repeated multiplication +section PowerTests + +open CPolynomial.Raw + +private abbrev p1x : CPolynomial.Raw ℤ := CPolynomial.Raw.mk #[1, 1] + +-- Zero polynomial: 0^n = 0 for n ≥ 1 +#guard (0 : CPolynomial.Raw ℤ) ^ 0 == C 1 +#guard (0 : CPolynomial.Raw ℤ) ^ 1 == (0 : CPolynomial.Raw ℤ) +#guard (0 : CPolynomial.Raw ℤ) ^ 5 == (0 : CPolynomial.Raw ℤ) + +-- Constant: (C 3)^n = C (3^n) +#guard (C 3 : CPolynomial.Raw ℤ) ^ 0 == C 1 +#guard (C 3 : CPolynomial.Raw ℤ) ^ 1 == C 3 +#guard (C 3 : CPolynomial.Raw ℤ) ^ 2 == C 9 +#guard (C 3 : CPolynomial.Raw ℤ) ^ 3 == C 27 + +-- Monomial: X^n +#guard (X : CPolynomial.Raw ℤ) ^ 0 == C 1 +#guard (X : CPolynomial.Raw ℤ) ^ 1 == (X : CPolynomial.Raw ℤ) +#guard (X : CPolynomial.Raw ℤ) ^ 3 == CPolynomial.Raw.mk #[(0 : ℤ), 0, 0, 1] + +-- Nontrivial: (1 + X)^n = binomial coefficients +#guard p1x ^ 0 == C 1 +#guard p1x ^ 1 == p1x +#guard p1x ^ 2 == CPolynomial.Raw.mk #[(1 : ℤ), 2, 1] +#guard p1x ^ 3 == CPolynomial.Raw.mk #[(1 : ℤ), 3, 3, 1] +#guard p1x ^ 4 == CPolynomial.Raw.mk #[(1 : ℤ), 4, 6, 4, 1] + +-- Agrees with repeated multiplication +#guard p1x ^ 3 == p1x * (p1x * p1x) +#guard p1x ^ 4 == p1x * (p1x * (p1x * p1x)) + +end PowerTests + end CompPoly diff --git a/tests/CompPolyTests/Univariate/Eval.lean b/tests/CompPolyTests/Univariate/Eval.lean new file mode 100644 index 00000000..646e2741 --- /dev/null +++ b/tests/CompPolyTests/Univariate/Eval.lean @@ -0,0 +1,123 @@ +/- +Copyright (c) 2026 CompPoly, Elias Judin, Harmonic. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Elias Judin, Aristotle (Harmonic) +-/ +import CompPoly.Univariate.ToPoly.Impl + +/-! + # Univariate Evaluation Regression Tests + + Regression tests for the Horner-style evaluation backend, covering: + - Zero polynomial + - Constant polynomials + - Monomial + - Trailing-zero polynomials + - Nontrivial polynomials + - Agreement with `toPoly` evaluation (Mathlib bridge) +-/ + +namespace CompPoly + +open CPolynomial.Raw + +-- ============================================================ +-- Raw evaluation: eval₂ with RingHom.id ≡ eval +-- ============================================================ + +section ZeroPoly + +-- eval of the zero polynomial should always be 0 +#guard (eval 42 (#[] : CPolynomial.Raw ℤ)) == 0 +#guard (eval 0 (#[] : CPolynomial.Raw ℤ)) == 0 +#guard (eval (-7) (#[] : CPolynomial.Raw ℤ)) == 0 + +end ZeroPoly + +section ConstPoly + +-- eval of a constant polynomial should return that constant +#guard (eval 100 (C 5 : CPolynomial.Raw ℤ)) == 5 +#guard (eval 0 (C 5 : CPolynomial.Raw ℤ)) == 5 +#guard (eval (-3) (C 0 : CPolynomial.Raw ℤ)) == 0 +#guard (eval 7 (C (-2) : CPolynomial.Raw ℤ)) == (-2) + +end ConstPoly + +section MonomialPoly + +-- eval of X at various points +#guard (eval 0 (X : CPolynomial.Raw ℤ)) == 0 +#guard (eval 1 (X : CPolynomial.Raw ℤ)) == 1 +#guard (eval 5 (X : CPolynomial.Raw ℤ)) == 5 +#guard (eval (-3) (X : CPolynomial.Raw ℤ)) == (-3) + +-- eval of X^2 = monomial 2 1 = [0, 0, 1] +#guard (eval 3 (monomial 2 1 : CPolynomial.Raw ℤ)) == 9 +#guard (eval (-2) (monomial 2 1 : CPolynomial.Raw ℤ)) == 4 + +end MonomialPoly + +section TrailingZeros + +-- Polynomial with trailing zeros: [1, 0, 0] represents constant 1 +-- Evaluation should still give 1 +#guard (eval 42 (mk #[(1 : ℤ), 0, 0])) == 1 +#guard (eval 0 (mk #[(1 : ℤ), 0, 0])) == 1 + +-- [0, 1, 0] represents X (with trailing zero) +#guard (eval 5 (mk #[(0 : ℤ), 1, 0])) == 5 +#guard (eval (-3) (mk #[(0 : ℤ), 1, 0])) == (-3) + +end TrailingZeros + +section NontrivialPoly + +-- 1 + X at x=3 should be 4 +private abbrev p1x : CPolynomial.Raw ℤ := mk #[1, 1] +#guard (eval 3 p1x) == 4 +#guard (eval 0 p1x) == 1 +#guard (eval (-1) p1x) == 0 + +-- 1 + 2X + 3X^2 at x=2 should be 1 + 4 + 12 = 17 +private abbrev p123 : CPolynomial.Raw ℤ := mk #[1, 2, 3] +#guard (eval 2 p123) == 17 +#guard (eval 0 p123) == 1 +#guard (eval 1 p123) == 6 +#guard (eval (-1) p123) == 2 + +-- 5 + 0*X + 3X^2 at x=2 should be 5 + 0 + 12 = 17 +private abbrev p503 : CPolynomial.Raw ℤ := mk #[5, 0, 3] +#guard (eval 2 p503) == 17 + +end NontrivialPoly + +section HornerNaiveAgreement + +-- Verify Horner eval₂ agrees with naive eval₂ on concrete examples +#guard (eval₂ (RingHom.id ℤ) 3 p1x) == (eval₂Naive (RingHom.id ℤ) 3 p1x) +#guard (eval₂ (RingHom.id ℤ) 2 p123) == (eval₂Naive (RingHom.id ℤ) 2 p123) +#guard (eval₂ (RingHom.id ℤ) (-1) p123) == (eval₂Naive (RingHom.id ℤ) (-1) p123) +#guard (eval₂ (RingHom.id ℤ) 0 (#[] : CPolynomial.Raw ℤ)) == + (eval₂Naive (RingHom.id ℤ) 0 (#[] : CPolynomial.Raw ℤ)) + +end HornerNaiveAgreement + +section ToPolyAgreement + +-- Verify eval agrees with toPoly evaluation via proved bridge theorems + +example : CPolynomial.Raw.eval 3 p1x = p1x.toPoly.eval 3 := by + rw [CPolynomial.Raw.eval_toPoly_eq_eval] + +example : CPolynomial.Raw.eval 2 p123 = p123.toPoly.eval 2 := by + rw [CPolynomial.Raw.eval_toPoly_eq_eval] + +-- eval₂ bridge +example : CPolynomial.Raw.eval₂ (RingHom.id ℤ) 2 p123 = + p123.toPoly.eval₂ (RingHom.id ℤ) 2 := by + exact CPolynomial.Raw.eval₂_toPoly (RingHom.id ℤ) 2 p123 + +end ToPolyAgreement + +end CompPoly