diff --git a/CompPoly.lean b/CompPoly.lean index d8c0bff9..aa13c9a6 100644 --- a/CompPoly.lean +++ b/CompPoly.lean @@ -99,6 +99,7 @@ import CompPoly.Univariate.BatchEval.Correctness import CompPoly.Univariate.BatchEval.Naive import CompPoly.Univariate.BatchEval.SubproductTree import CompPoly.Univariate.CMvEquiv +import CompPoly.Univariate.Context import CompPoly.Univariate.Deriv import CompPoly.Univariate.DivisionCorrectness import CompPoly.Univariate.EuclideanAlgorithm @@ -107,6 +108,7 @@ import CompPoly.Univariate.Linear import CompPoly.Univariate.ManyEval import CompPoly.Univariate.ManyEval.Basic import CompPoly.Univariate.ManyEval.Correctness +import CompPoly.Univariate.Modular import CompPoly.Univariate.NTT.BabyBear import CompPoly.Univariate.NTT.Domain import CompPoly.Univariate.NTT.Evaluation @@ -133,10 +135,23 @@ import CompPoly.Univariate.NTTFast.Plan import CompPoly.Univariate.Quotient.Core import CompPoly.Univariate.Quotient.Equiv import CompPoly.Univariate.Raw +import CompPoly.Univariate.Raw.Context import CompPoly.Univariate.Raw.Core import CompPoly.Univariate.Raw.Division +import CompPoly.Univariate.Raw.Modular import CompPoly.Univariate.Raw.Ops import CompPoly.Univariate.Raw.Proofs +import CompPoly.Univariate.Roots +import CompPoly.Univariate.Roots.Backend +import CompPoly.Univariate.Roots.Context +import CompPoly.Univariate.Roots.Correctness +import CompPoly.Univariate.Roots.Enumeration +import CompPoly.Univariate.Roots.Extraction +import CompPoly.Univariate.Roots.RootProduct +import CompPoly.Univariate.Roots.SmoothSubgroup +import CompPoly.Univariate.Roots.SmoothSubgroup.Basic +import CompPoly.Univariate.Roots.SmoothSubgroup.Correctness +import CompPoly.Univariate.Roots.Splitter import CompPoly.Univariate.ToPoly import CompPoly.Univariate.ToPoly.Core import CompPoly.Univariate.ToPoly.Degree diff --git a/CompPoly/Bivariate/ToPoly.lean b/CompPoly/Bivariate/ToPoly.lean index e031cb63..413f7ec6 100644 --- a/CompPoly/Bivariate/ToPoly.lean +++ b/CompPoly/Bivariate/ToPoly.lean @@ -210,7 +210,7 @@ theorem ofPoly_coeff {R : Type*} [BEq R] [LawfulBEq R] [Nontrivial R] [Semiring /-- The outer coefficient of `toPoly p` is `CPolynomial.coeff p n` converted to `R[X]`. -/ theorem toPoly_coeff {R : Type*} [BEq R] [LawfulBEq R] [Semiring R] (p : CBivariate R) (n : ℕ) : (toPoly p).coeff n = (CPolynomial.coeff p n).toPoly := by - rw [ CBivariate.toPoly, Polynomial.finset_sum_coeff ] + rw [ CBivariate.toPoly, Polynomial.finsetSum_coeff ] rw [ Finset.sum_eq_single n ] <;> simp +contextual [ Polynomial.coeff_monomial ] simp_all +decide [ CPolynomial.mem_support_iff ] aesop @@ -378,7 +378,7 @@ theorem evalY_toPoly {R : Type*} [BEq R] [LawfulBEq R] [Nontrivial R] [Semiring simp only [Array.length_toList, Array.getElem?_toList] symm unfold CBivariate.toPoly - simp +decide [ Polynomial.eval_finset_sum ] + simp +decide [ Polynomial.eval_finsetSum ] rw [ 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)) @@ -468,7 +468,7 @@ theorem natDegreeY_toPoly {R : Type*} [BEq R] [LawfulBEq R] [Nontrivial R] [Semi theorem coeff_toPoly_Y {R : Type*} [BEq R] [LawfulBEq R] [Nontrivial R] [Semiring R] (f : CBivariate R) (j : ℕ) : (toPoly f).coeff j = CPolynomial.toPoly (f.val.coeff j) := by - erw [ Polynomial.finset_sum_coeff ] + erw [ Polynomial.finsetSum_coeff ] rw [ Finset.sum_eq_single j ] <;> simp +contextual [ Polynomial.coeff_monomial ] intro hj rw [ CPolynomial.support ] at hj @@ -628,7 +628,7 @@ theorem evalX_toPoly_eval_commute {R : Type*} [BEq R] [LawfulBEq R] [Nontrivial Polynomial.eval a (∑ j ∈ s, g j * Polynomial.C (y ^ j)) = ∑ j ∈ s, Polynomial.eval a (g j) * y ^ j := by exact fun s g hg => by - rw [Polynomial.eval_finset_sum, + rw [Polynomial.eval_finsetSum, Finset.sum_congr rfl (fun j hj => h_eval_mul_C _ _ (hg j hj))] convert h_sum _ _ _ · simp +decide [Polynomial.eval_eq_sum, Polynomial.sum_def] diff --git a/CompPoly/Data/Array/Lemmas.lean b/CompPoly/Data/Array/Lemmas.lean index f62f4f62..b7659143 100644 --- a/CompPoly/Data/Array/Lemmas.lean +++ b/CompPoly/Data/Array/Lemmas.lean @@ -17,6 +17,16 @@ namespace Array variable {α : Type*} {unit : α} +/-- Remove duplicates from an array while preserving first occurrences. -/ +def eraseDups [BEq α] (xs : Array α) : Array α := + xs.foldl + (fun out x ↦ + if out.contains x then + out + else + out.push x) + #[] + /-- Checks if an array of elements from a type `R` is a boolean array, i.e., if every element is either `0` or `1`. -/ def isBoolean {R : Type _} [Zero R] [One R] (a : Array R) : Prop := @@ -65,6 +75,68 @@ lemma getElem?_eq_toList {a : Array α} {i : ℕ} : a.toList[i]? = a[i]? := by rw (occs := .pos [2]) [← Array.toArray_toList (xs := a)] rw [List.getElem?_toArray] +/-- `Array.map` and `getD` agree with indexed access on in-bounds indices. -/ +theorem getD_map_of_lt {β : Type*} (xs : Array α) (f : α → β) + (d : β) {i : Nat} (hi : i < xs.size) : + (xs.map f).getD i d = f xs[i] := by + rw [Array.getD_eq_getD_getElem?, Array.getElem?_map, Array.getElem?_eq_getElem hi] + simp + +theorem foldl_zipIdx_eq_foldl_toList_zipIdx {β : Type*} + (f : β → α × Nat → β) (init : β) (a : Array α) : + a.zipIdx.foldl f init = a.toList.zipIdx.foldl f init := by + cases a + simp + +theorem foldl_zipIdx_eq_foldl_toList_zipIdx_size {β : Type*} + (f : β → α × Nat → β) (init : β) (a : Array α) : + Array.foldl f init a.zipIdx 0 a.size = a.toList.zipIdx.foldl f init := by + cases a + simp + +theorem mem_foldl_append_of_mem {β : Type*} + (xs : Array α) (f : α → Array β) {x : α} {y : β} + (hx : x ∈ xs.toList) (hy : y ∈ (f x).toList) : + y ∈ (xs.foldl (fun out x ↦ out ++ f x) #[]).toList := by + cases xs with + | mk data => + simp at hx ⊢ + have haux : ∀ (data : List α) (acc : Array β), + y ∈ acc.toList ∨ (∃ x, x ∈ data ∧ y ∈ (f x).toList) → + y ∈ (data.foldl (fun out x ↦ out ++ f x) acc).toList := by + intro data + induction data with + | nil => + intro acc h + simp at h ⊢ + exact h + | cons z zs ih => + intro acc h + simp only [List.foldl_cons] + apply ih + rcases h with hacc | ⟨x, hx, hyx⟩ + · left + simp [hacc] + · simp only [List.mem_cons] at hx + cases hx with + | inl hxz => + subst x + left + simp [hyx] + | inr hxzs => + right + exact ⟨x, hxzs, hyx⟩ + simpa using haux data #[] (Or.inr ⟨x, hx, hy⟩) + +theorem mem_flatten_map_of_mem {β : Type*} + (xs : Array α) (f : α → Array β) {x : α} {y : β} + (hx : x ∈ xs.toList) (hy : y ∈ (f x).toList) : + y ∈ (xs.map f).flatten := by + cases xs with + | mk data => + simp at hx hy ⊢ + exact ⟨x, hx, hy⟩ + attribute [simp] Array.getElem?_eq_getElem -- @[simp] theorem matchSize_comm (a : Array α) (b : Array α) (unit : α) : diff --git a/CompPoly/Data/List/Lemmas.lean b/CompPoly/Data/List/Lemmas.lean index 25d325a0..9f22c96d 100644 --- a/CompPoly/Data/List/Lemmas.lean +++ b/CompPoly/Data/List/Lemmas.lean @@ -29,6 +29,19 @@ theorem append_getLast_dropLast {α : Type u} (l : List α) (h : l ≠ []) : simp only [cons_append, cons.injEq, true_and] apply ih +theorem foldl_congr_of_mem {α : Type u} {β : Type v} {f g : α → β → α} + (xs : List β) (acc : α) (h : ∀ acc' x, x ∈ xs → f acc' x = g acc' x) : + xs.foldl f acc = xs.foldl g acc := by + induction xs generalizing acc with + | nil => + simp + | cons x xs ih => + simp only [foldl_cons] + rw [h acc x (by simp)] + apply ih + intro acc' y hy + exact h acc' y (by simp [hy]) + theorem foldl_split_outer {α : Type u} {β : Type v} (f : α → β → α) (init : α) (l : List β) (h : l ≠ []): foldl (f:=f) (init:=init) (l) = f (foldl (f:=f) (init:=init) (l.dropLast)) (l.getLast (by omega)) := by diff --git a/CompPoly/Data/Nat/Bitwise.lean b/CompPoly/Data/Nat/Bitwise.lean index d2325dc0..a82b3bba 100644 --- a/CompPoly/Data/Nat/Bitwise.lean +++ b/CompPoly/Data/Nat/Bitwise.lean @@ -547,9 +547,8 @@ lemma sum_eq_xor_plus_twice_and (n : Nat) : ∀ m : ℕ, n + m = (n ^^^ m) + 2 * rw [zero_add, Nat.zero_and, mul_zero, add_zero, Nat.zero_xor] | bit bn n2 ih => intro m - let resDiv2M := Nat.boddDiv2 m - let bm := resDiv2M.fst - let m2 := resDiv2M.snd + let bm := Nat.bodd m + let m2 := Nat.div2 m have h_m2 : m2 = Nat.div2 m := by rfl have h_bm : bm = Nat.bodd m := by @@ -1295,7 +1294,7 @@ lemma getBit_of_binaryFinMapToNat {n : ℕ} (m : Fin n → ℕ) (h_binary: ∀ j induction n with | zero => intro k; - simp only [Nat.pow_zero, Fin.val_eq_zero, not_lt_zero', ↓reduceDIte] + simp only [Nat.pow_zero, Fin.val_eq_zero, Nat.not_lt_zero, ↓reduceDIte] exact getBit_zero_eq_zero | succ n ih => -- Inductive step: Assume the property holds for `n`, prove it for `n+1`. diff --git a/CompPoly/Data/Polynomial/MonomialBasis.lean b/CompPoly/Data/Polynomial/MonomialBasis.lean index a48928e1..adb8002c 100644 --- a/CompPoly/Data/Polynomial/MonomialBasis.lean +++ b/CompPoly/Data/Polynomial/MonomialBasis.lean @@ -64,7 +64,7 @@ noncomputable def monomialBasisOfDegreeLT {n : ℕ} : Basis (Fin n) L (L⦃< n simp only [SetLike.mk_smul_mk] conv => lhs - simp only [Finsupp.sum, AddSubmonoidClass.coe_finset_sum, finset_sum_coeff, coeff_smul, + simp only [Finsupp.sum, AddSubmonoidClass.coe_finsetSum, finsetSum_coeff, coeff_smul, coeff_X_pow, smul_eq_mul, mul_ite, mul_one, mul_zero, monomials_in_submodule, v] -- ⊢ (∑ x ∈ l.support, if ↑i = ↑x then l x else 0) = l i simp_rw [Fin.val_eq_val, eq_comm] @@ -110,7 +110,7 @@ noncomputable def monomialBasisOfDegreeLT {n : ℕ} : Basis (Fin n) L (L⦃< n -- we have to convert (↑p).sum into Fin n → L form using Polynomial.sum_fin simp only [monomial_zero_right, implies_true, ←Polynomial.sum_fin (hn := h_deg_p)] -- ⊢ ↑(∑ i, c i • monomials_in_submodule i) = ∑ i, (monomial ↑i) ((↑p).coeff ↑i) - rw [AddSubmonoidClass.coe_finset_sum] -- bring both sides back to L[X] + rw [AddSubmonoidClass.coe_finsetSum] -- bring both sides back to L[X] apply Finset.sum_congr rfl intro ⟨i, hi_finN⟩ hi simp only [SetLike.mk_smul_mk, c, monomials_in_submodule] diff --git a/CompPoly/Data/RingTheory/AlgebraTower.lean b/CompPoly/Data/RingTheory/AlgebraTower.lean index 4195a7b3..17a2cfeb 100644 --- a/CompPoly/Data/RingTheory/AlgebraTower.lean +++ b/CompPoly/Data/RingTheory/AlgebraTower.lean @@ -61,7 +61,7 @@ instance AlgebraTower.toIsScalarTower (a : AlgebraTower C) {i j k : ι} simp only [map_mul] rw [←RingHom.comp_apply] unfold instIJ instJK instIK AlgebraTower.toAlgebra - simp_rw [algebraMap, Algebra.algebraMap] + simp_rw [RingHom.algebraMap_toAlgebra] have h_assoc := a.coherence' (i:=i) (j:=j) (k:=k) (h1:=h1) (h2:=h2) rw [h_assoc] rw [mul_assoc] @@ -123,8 +123,7 @@ def AlgebraTowerEquiv.toAlgEquivOverLeft (e : AlgebraTowerEquiv A B) (i j : ι) commutes' := fun r => by simp only [RingEquiv.toEquiv_eq_coe, Equiv.toFun_as_coe, EquivLike.coe_coe] unfold instAij instAiBij - rw [algebraMap, algebraMap, Algebra.algebraMap, Algebra.algebraMap,AlgebraTower.toAlgebra, - AlgebraTowerEquiv.toAlgebraOverLeft, AlgebraTowerEquiv.algebraMapRightUp] + simp only [AlgebraTowerEquiv.algebraMapRightUp, RingHom.algebraMap_toAlgebra] simp only [RingEquiv.toRingHom_eq_coe, RingHom.coe_comp, RingHom.coe_coe, Function.comp_apply] exact Eq.symm (e.commutesLeft' i j h r) map_mul' := fun x y => by diff --git a/CompPoly/Fields/Binary/AdditiveNTT/Intermediate.lean b/CompPoly/Fields/Binary/AdditiveNTT/Intermediate.lean index 4cc5370f..cf8edf65 100644 --- a/CompPoly/Fields/Binary/AdditiveNTT/Intermediate.lean +++ b/CompPoly/Fields/Binary/AdditiveNTT/Intermediate.lean @@ -277,7 +277,7 @@ lemma getSDomainBasisCoeff_of_sum_repr [NeZero R_rate] (i : Fin (ℓ + 1)) simp only [Basis.repr_self] have hx_at_j_simplified : (∑ j_x, x_coeffs j_x • (b.repr (b j_x))) j = x_coeffs j := by - simp only [h_repr_basis, Finsupp.smul_single, smul_eq_mul, mul_one, Finsupp.coe_finset_sum, + simp only [h_repr_basis, Finsupp.smul_single, smul_eq_mul, mul_one, Finsupp.coe_finsetSum, Finset.sum_apply, Finsupp.single_apply, Finset.sum_ite_eq', Finset.mem_univ, ↓reduceIte] let x_coeffs_fs := Finsupp.equivFunOnFinite.symm x_coeffs let rhs_sum := ∑ j_x, (x_coeffs_fs j_x) • (b j_x) diff --git a/CompPoly/Fields/Binary/AdditiveNTT/NovelPolynomialBasis.lean b/CompPoly/Fields/Binary/AdditiveNTT/NovelPolynomialBasis.lean index f504e7f7..f47a7d04 100644 --- a/CompPoly/Fields/Binary/AdditiveNTT/NovelPolynomialBasis.lean +++ b/CompPoly/Fields/Binary/AdditiveNTT/NovelPolynomialBasis.lean @@ -189,7 +189,7 @@ lemma βᵢ_not_in_Uᵢ (i : Fin r) : omega else push Not at h_i - have h_i_eq_0: i = 0 := by exact Fin.le_zero_iff'.mp h_i + have h_i_eq_0: i = 0 := by exact nonpos_iff_eq_zero.mp h_i have set_empty: Set.Ico 0 i = ∅ := by rw [h_i_eq_0] simp only [Set.Ico_eq_empty_iff] diff --git a/CompPoly/Fields/Binary/Common.lean b/CompPoly/Fields/Binary/Common.lean index b4c61bc5..6e1189ee 100644 --- a/CompPoly/Fields/Binary/Common.lean +++ b/CompPoly/Fields/Binary/Common.lean @@ -360,7 +360,7 @@ lemma toPoly_ne_zero_iff_ne_zero {w : Nat} (v : BitVec w) : have h_coeff_formula : ∀ i : Fin w, (toPoly v).coeff i.val = if v.getLsb i then 1 else 0 := by intro i unfold toPoly - rw [finset_sum_coeff] + rw [finsetSum_coeff] have h_ite: ∀ b : Fin w, (if v.getLsb b then (X : (ZMod 2)[X])^b.val else 0).coeff i.val = if v.getLsb b then (if b = i then 1 else 0) else 0 := by intro b @@ -470,7 +470,7 @@ lemma BitVec_lt_two_pow_of_toPoly_degree_lt {w d : ℕ} (v : BitVec w) have h_i_lt_w : i < w := by omega have h_coeff_one : (toPoly v).coeff i = 1 := by unfold toPoly - rw [finset_sum_coeff] + rw [finsetSum_coeff] rw [Finset.sum_eq_single ⟨i, h_i_lt_w⟩] · have h_getLsb_eq_testBit : v.getLsb ⟨i, h_i_lt_w⟩ = v.toNat.testBit i := by rfl @@ -627,7 +627,7 @@ theorem toPoly_coeff {w : ℕ} (v : BitVec w) (n : ℕ) : (if h : n < w then (if v.getLsb ⟨n, h⟩ = true then (1 : ZMod 2) else 0) else 0) := by classical unfold toPoly - rw [Polynomial.finset_sum_coeff] + rw [Polynomial.finsetSum_coeff] by_cases h : n < w · -- case n < w simp only [dif_pos h] @@ -667,7 +667,7 @@ theorem toPoly_coeff {w : ℕ} (v : BitVec w) (n : ℕ) : intro hi0_not simp [Finset.mem_univ] at hi0_not -- conclude - simpa [Polynomial.finset_sum_coeff, i0] using hmain + simpa [Polynomial.finsetSum_coeff, i0] using hmain · -- case ¬ n < w simp only [dif_neg h] -- show every summand has coefficient 0 diff --git a/CompPoly/Fields/Binary/Tower/Abstract/Algebra.lean b/CompPoly/Fields/Binary/Tower/Abstract/Algebra.lean index e9a35b6a..9c35e455 100644 --- a/CompPoly/Fields/Binary/Tower/Abstract/Algebra.lean +++ b/CompPoly/Fields/Binary/Tower/Abstract/Algebra.lean @@ -240,7 +240,6 @@ lemma algebra_adjacent_tower_eq_AdjoinRoot_algebra (k : ℕ) : rw [algebra_adjacent_tower_def] unfold canonicalEmbedding rw [←AdjoinRoot.algebraMap_eq] - rw [algebraMap, Algebra.algebraMap] exact Algebra.algebra_ext (AdjoinRoot.instAlgebra (poly k)).2.toAlgebra (AdjoinRoot.instAlgebra (poly k)) (congrFun rfl) diff --git a/CompPoly/Fields/Binary/Tower/Abstract/Basis.lean b/CompPoly/Fields/Binary/Tower/Abstract/Basis.lean index 3e358a74..1c6ce968 100644 --- a/CompPoly/Fields/Binary/Tower/Abstract/Basis.lean +++ b/CompPoly/Fields/Binary/Tower/Abstract/Basis.lean @@ -174,7 +174,7 @@ lemma algebraMap_𝕏_eq_of_index_eq (r k m : ℕ) (h_k_le : k + 1 ≤ r) (h_m_l (h_eq : k = m) : letI := binaryAlgebraTower (l := k + 1) (r := r) (h_le := h_k_le) letI := binaryAlgebraTower (l := m + 1) (r := r) (h_le := h_m_le) - (Algebra.algebraMap (𝕏 k) : BTField r) = (Algebra.algebraMap (𝕏 m) : BTField r) := by + (algebraMap _ _ (𝕏 k) : BTField r) = (algebraMap _ _ (𝕏 m) : BTField r) := by subst h_eq rfl @@ -302,13 +302,11 @@ theorem multilinearBasis_apply (r : ℕ) : ∀ l : ℕ, (h_le : l ≤ r) → ∀ simp_rw [algebraMap.coe_pow] -- rhs simp_rw [algebraMap.coe_prod] -- lhs unfold Algebra.cast - rw! (castMode:=.all) [←algebraMap] conv_lhs => rw [←Fin.prod_congr' (b:=r1-l) (a:=prevDiff) (h:=by omega)] simp only [Fin.val_cast] - simp (config := { failIfUnchanged := false }) only [algebraMap, instAlgebraSucc] - erw [RingHom.map_pow] - simp (config := { failIfUnchanged := false }) only [←binaryTowerAlgebra_apply_assoc] + simp only [RingHom.map_pow] + simp only [←binaryTowerAlgebra_apply_assoc] ------------------ Equality of bit-based powers of generators ----------------- --- The outtermost term have hfinProd_msb := bit_revFinProdFinEquiv_symm_2_pow_succ (n:=prevDiff) diff --git a/CompPoly/Fields/Binary/Tower/Abstract/Split.lean b/CompPoly/Fields/Binary/Tower/Abstract/Split.lean index 50934803..8e6946e0 100644 --- a/CompPoly/Fields/Binary/Tower/Abstract/Split.lean +++ b/CompPoly/Fields/Binary/Tower/Abstract/Split.lean @@ -149,10 +149,10 @@ theorem mul_join_via_add_smul (k : ℕ) (h_pos : k > 0) (a₁ a₀ b₁ b₀ : B have h_a₁_b₀_Z_k : (algebraMap (BTField (k - 1)) (BTField k)) b₀ * a₁ • Z k = (a₁ * b₀) • Z k := by - rw [Algebra.smul_def', ←algebraMap, ←mul_assoc, ←map_mul, ←Algebra.smul_def, mul_comm] + rw [Algebra.smul_def', ←mul_assoc, ←map_mul, ←Algebra.smul_def, mul_comm] have h_a₀_b₁_Z_k : (algebraMap (BTField (k - 1)) (BTField k)) a₀ * b₁ • Z k = (a₀ * b₁) • Z k := by - rw [Algebra.smul_def', ←algebraMap, ←mul_assoc, ←map_mul, ←Algebra.smul_def, mul_comm] + rw [Algebra.smul_def', ←mul_assoc, ←map_mul, ←Algebra.smul_def, mul_comm] have h_Z_k_pow_2 : (Z k) ^ 2 = Z (k - 1) • Z k + 1 := by rw [sumZeroIffEq (x:=(Z k)^2) (y:=Z (k - 1) • Z k + 1).mp] rw [←add_assoc] @@ -161,7 +161,6 @@ theorem mul_join_via_add_smul (k : ℕ) (h_pos : k > 0) (a₁ a₀ b₁ b₀ : B rw! (castMode:=.all) [Nat.sub_one_add_one (by omega)] at h simp only [eq_mp_eq_cast] at h convert h - rw [Algebra.algebraMap] conv_lhs => simp only [instAlgebra]; change (towerAlgebraMap (l:=k-1) (r:=k) (h_le:=by omega)) (Z (k - 1)) @@ -307,9 +306,9 @@ lemma split_algebraMap_eq_zero_x {k : ℕ} (h_pos : k > 0) (x : BTField (k - 1)) apply h.mp -- ⊢ mappedVal = join_via_add_smul h_pos 0 x unfold mappedVal - rw [algebraMap, Algebra.algebraMap] - unfold instAlgebra binaryAlgebraTower - rw [AlgebraTower.toAlgebra, AlgebraTower.algebraMap, instAlgebraTowerNatBTField] + unfold instAlgebra binaryAlgebraTower AlgebraTower.toAlgebra + simp only [RingHom.algebraMap_toAlgebra] + rw [AlgebraTower.algebraMap, instAlgebraTowerNatBTField] simp only have h_concrete_embedding_succ_1 := towerAlgebraMap_succ_1 (k:=k-1) rw! (castMode:=.all) [Nat.sub_one_add_one (by omega)] at h_concrete_embedding_succ_1 @@ -327,12 +326,11 @@ lemma split_algebraMap_eq_zero_x {k : ℕ} (h_pos : k > 0) (x : BTField (k - 1)) have h := algebraMap_adjacent_tower_def (l:=k-1) rw! (castMode:=.all) [Nat.sub_one_add_one (by omega)] at h simp only [eqRec_eq_cast] at h - rw [algebraMap, Algebra.algebraMap] at ⊢ h rw! (castMode:=.all) [Nat.sub_one_add_one (by omega)] at h simp only [cast_eq] at h unfold binaryAlgebraTower AlgebraTower.toAlgebra AlgebraTower.algebraMap instAlgebraTowerNatBTField - simp only [] -- normalize eqRec before rewrite + simp only [RingHom.algebraMap_toAlgebra] -- unfold algebraMap (v4.30: no longer rw-unfoldable) -- Both sides reduce to (cast ⋯ (canonicalEmbedding (k-1))) x through different paths erw [h_concrete_embedding_succ_1]; simp only [eqRec_eq_cast] diff --git a/CompPoly/Fields/Binary/Tower/Concrete/Algebra.lean b/CompPoly/Fields/Binary/Tower/Concrete/Algebra.lean index 4baac08e..f250965a 100644 --- a/CompPoly/Fields/Binary/Tower/Concrete/Algebra.lean +++ b/CompPoly/Fields/Binary/Tower/Concrete/Algebra.lean @@ -234,9 +234,9 @@ lemma split_algebraMap_eq_zero_x {k : ℕ} (h_pos : k > 0) (x : ConcreteBTField apply h -- ⊢ mappedVal = join h_pos zero x unfold mappedVal - rw [algebraMap, Algebra.algebraMap] - unfold instAlgebra ConcreteBTFieldAlgebra - rw [AlgebraTower.toAlgebra, AlgebraTower.algebraMap, instAlgebraTowerConcreteBTF] + unfold instAlgebra ConcreteBTFieldAlgebra AlgebraTower.toAlgebra + simp only [RingHom.algebraMap_toAlgebra] + rw [AlgebraTower.algebraMap, instAlgebraTowerConcreteBTF] simp only have h_concrete_embedding_succ_1 := concreteTowerAlgebraMap_succ_1 (k:=k-1) rw! (castMode:=.all) [Nat.sub_one_add_one (by omega)] at h_concrete_embedding_succ_1 @@ -439,8 +439,7 @@ lemma aeval_definingPoly_at_Z_succ (k : ℕ) : -- Z_square_mul_form uses instAlgebraLiftConcreteBTField internally rw [Z_square_mul_form (k:=k) (prev:=(getBTFResult (k:=k)))] rw [add_assoc] - rw [algebraMap, Algebra.algebraMap, instAlgebraLiftConcreteBTField] - simp only + simp only [RingHom.algebraMap_toAlgebra] -- f uses ConcreteBTFieldAlgebra, it's same as instAlgebraLiftConcreteBTField at step = 1 rw [h_f_is_canonical_embedding, concreteTowerAlgebraMap_succ_1] simp only [canonicalAlgMap]; rw [mul_comm] diff --git a/CompPoly/Fields/Binary/Tower/Concrete/Basis.lean b/CompPoly/Fields/Binary/Tower/Concrete/Basis.lean index 70750465..7c189ab9 100644 --- a/CompPoly/Fields/Binary/Tower/Concrete/Basis.lean +++ b/CompPoly/Fields/Binary/Tower/Concrete/Basis.lean @@ -136,8 +136,7 @@ def basisSucc (k : ℕ) : Basis (Fin 2) (ConcreteBTField k) (ConcreteBTField (k rw [←h_add_smul] unfold join_via_add_smul simp only [Nat.add_one_sub_one] - rw [algebraMap, Algebra.algebraMap, ConcreteBTFieldAlgebra_def] - simp only + simp only [RingHom.algebraMap_toAlgebra] simp only [generator] rw [add_comm] congr -- .Q.E.D @@ -422,8 +421,8 @@ lemma algebraMap_𝕏_eq_of_index_eq (r k m : ℕ) (h_k_le : k + 1 ≤ r) (h_m_l (h_eq : k = m) : letI := ConcreteBTFieldAlgebra (l := k + 1) (r := r) (h_le := h_k_le) letI := ConcreteBTFieldAlgebra (l := m + 1) (r := r) (h_le := h_m_le) - (Algebra.algebraMap (𝕏 k) : ConcreteBTField r) = - (Algebra.algebraMap (𝕏 m) : ConcreteBTField r) := by + (algebraMap _ _ (𝕏 k) : ConcreteBTField r) = + (algebraMap _ _ (𝕏 m) : ConcreteBTField r) := by subst h_eq rfl @@ -548,13 +547,11 @@ theorem multilinearBasis_apply (r : ℕ) : ∀ l : ℕ, (h_le : l ≤ r) → ∀ simp_rw [algebraMap.coe_pow] simp_rw [algebraMap.coe_prod] unfold Algebra.cast - rw! (castMode:=.all) [←algebraMap] conv_lhs => rw [←Fin.prod_congr' (b:=r1 - l) (a:=prevDiff) (h:=by omega)] simp only [Fin.val_cast] - simp (config := { failIfUnchanged := false }) only [algebraMap, instAlgebraSucc] - erw [RingHom.map_pow] - simp (config := { failIfUnchanged := false }) only [←ConcreteBTFieldAlgebra_apply_assoc] + simp only [RingHom.map_pow] + simp only [←ConcreteBTFieldAlgebra_apply_assoc] ------------------ Equality of bit-based powers of generators ----------------- have hfinProd_msb := bit_revFinProdFinEquiv_symm_2_pow_succ (n:=prevDiff) (i:=⟨prevDiff, by omega⟩) (j:=⟨j, by omega⟩) diff --git a/CompPoly/Fields/KoalaBear/Basic.lean b/CompPoly/Fields/KoalaBear/Basic.lean index 44288729..e32de202 100644 --- a/CompPoly/Fields/KoalaBear/Basic.lean +++ b/CompPoly/Fields/KoalaBear/Basic.lean @@ -164,7 +164,7 @@ lemma twoAdicity_maximal : ¬ (2 ^ (twoAdicity + 1)) ∣ (fieldSize - 1) := by decide /-- Repeated squaring: `sqChain g n = g ^ (2^n)`. - Does `n` multiplications instead of `2^n`, making it kernel-friendly. -/ + Uses `n` multiplications and avoids expanding the exponent into `2^n` steps. -/ private def sqChain (g : Field) : Nat → Field | 0 => g | n + 1 => let h := sqChain g n; h * h @@ -286,4 +286,68 @@ lemma twoAdicGenerators_order (bits : Fin (twoAdicity + 1)) : (twoAdicGenerators_pow_twoPow_ne_one_of_lt (bits := ⟨n + 1, hb⟩) (m := n) (by simp)) (twoAdicGenerators_pow_twoPow_eq_one ⟨n + 1, hb⟩) +/-- Primitive generator used by the smooth field-root splitter. -/ +def primitiveRoot : Field := (3 : Field) + +set_option maxRecDepth 100000 in +/-- `primitiveRoot ^ 127` is the maximal two-adic generator. -/ +private lemma primitiveRoot_pow_127_eq_twoAdicGenerator : + primitiveRoot ^ 127 = + twoAdicGenerators[(⟨twoAdicity, by omega⟩ : Fin (twoAdicity + 1))] := by + unfold primitiveRoot twoAdicity + decide + +/-- `primitiveRoot ^ 2^twoAdicity` is nontrivial. -/ +private lemma primitiveRoot_pow_twoAdicity_ne_one : + primitiveRoot ^ (2 ^ twoAdicity) ≠ (1 : Field) := by + rw [← sqChain_eq_pow_two_pow] + unfold primitiveRoot twoAdicity + decide + +/-- Prime divisors of `fieldSize - 1` are exactly `2` and `127`. -/ +private lemma prime_dvd_fieldSize_sub_one_cases {p : Nat} + (hp : p.Prime) (hdvd : p ∣ fieldSize - 1) : + p = 2 ∨ p = 127 := by + have hpdvd : p ∣ 2 ^ twoAdicity * 127 := by + rw [← fieldSize_sub_one_factorization] + exact hdvd + rcases hp.dvd_mul.mp hpdvd with h2pow | h127 + · left + exact (Nat.prime_dvd_prime_iff_eq hp Nat.prime_two).mp + (hp.dvd_of_dvd_pow h2pow) + · right + exact (Nat.prime_dvd_prime_iff_eq hp (by decide : Nat.Prime 127)).mp h127 + +/-- The smooth field-root splitter generator has full multiplicative order. -/ +lemma primitiveRoot_order : orderOf primitiveRoot = fieldSize - 1 := by + refine orderOf_eq_of_pow_and_pow_div_prime (n := fieldSize - 1) ?_ ?_ ?_ + · unfold fieldSize + omega + · exact ZMod.pow_card_sub_one_eq_one (a := primitiveRoot) (by + unfold primitiveRoot + decide) + · intro p hp hdvd + rcases prime_dvd_fieldSize_sub_one_cases hp hdvd with rfl | rfl + · rw [fieldSize_sub_one_factorization, twoAdicity] + have hdiv : (2 ^ 24 * 127) / 2 = 127 * 2 ^ 23 := by decide + rw [hdiv, pow_mul, primitiveRoot_pow_127_eq_twoAdicGenerator] + exact twoAdicGenerators_pow_twoPow_ne_one_of_lt + (bits := (⟨twoAdicity, by omega⟩ : Fin (twoAdicity + 1))) (m := 23) + (by simp [twoAdicity]) + · rw [fieldSize_sub_one_factorization] + have hdiv : (2 ^ twoAdicity * 127) / 127 = 2 ^ twoAdicity := + by decide + rw [hdiv] + exact primitiveRoot_pow_twoAdicity_ne_one + +/-- Smooth subgroup refinement schedule for `fieldSize - 1 = 2^24 * 127`. -/ +def smoothRootSchedule : Array Nat := + (Array.replicate twoAdicity 2).push 127 + +/-- The KoalaBear smooth schedule refines the multiplicative group down to singleton cosets. -/ +lemma smoothRootSchedule_fold_eq_one : + smoothRootSchedule.toList.foldl (fun order ell ↦ order / ell) (fieldSize - 1) = 1 := by + unfold smoothRootSchedule fieldSize twoAdicity + decide + end KoalaBear diff --git a/CompPoly/Multivariate/CMvPolynomialEvalLemmas.lean b/CompPoly/Multivariate/CMvPolynomialEvalLemmas.lean index 6593d711..56fd6134 100644 --- a/CompPoly/Multivariate/CMvPolynomialEvalLemmas.lean +++ b/CompPoly/Multivariate/CMvPolynomialEvalLemmas.lean @@ -106,11 +106,10 @@ theorem CMvPolynomial.eval_ext_univariate have hagreeUni : d < (S.filter (fun r ↦ pUni.eval r = qUni.eval r)).card := by convert hagree using 3 - · ext r - rw [show pUni = (CompPoly.CPolynomial.cmvEquiv (R := R)).symm p from rfl, - show qUni = (CompPoly.CPolynomial.cmvEquiv (R := R)).symm q from rfl, - CompPoly.CPolynomial.eval_cmvEquiv_symm, - CompPoly.CPolynomial.eval_cmvEquiv_symm] + rw [show pUni = (CompPoly.CPolynomial.cmvEquiv (R := R)).symm p from rfl, + show qUni = (CompPoly.CPolynomial.cmvEquiv (R := R)).symm q from rfl, + CompPoly.CPolynomial.eval_cmvEquiv_symm, + CompPoly.CPolynomial.eval_cmvEquiv_symm] have hUni : pUni = qUni := CompPoly.CPolynomial.eval_ext (p := pUni) (q := qUni) hdegUni hagreeUni exact (CompPoly.CPolynomial.cmvEquiv (R := R)).symm.injective hUni diff --git a/CompPoly/ToMathlib/MvPolynomial/Equiv.lean b/CompPoly/ToMathlib/MvPolynomial/Equiv.lean index ab4527e9..79278f42 100644 --- a/CompPoly/ToMathlib/MvPolynomial/Equiv.lean +++ b/CompPoly/ToMathlib/MvPolynomial/Equiv.lean @@ -111,11 +111,11 @@ theorem eval_eq_eval_mv_eval_finSuccEquivNth (s : Fin n → R) (y : R) congr 2 apply MvPolynomial.algHom_ext simp only [Fin.forall_iff_succAbove p, aeval_X, Fin.insertNth_apply_same, Polynomial.mapAlgHom, - AlgHom.toRingHom_eq_coe, coe_aeval_eq_eval, AlgEquiv.toAlgHom_eq_coe, AlgHom.coe_comp, - Polynomial.coe_aeval_eq_eval, AlgHom.coe_mk, coe_mapRingHom, AlgHom.coe_coe, comp_apply, - finSuccEquivNth_apply, eval₂Hom_X', Polynomial.map_X, Polynomial.eval_X, - Fin.insertNth_apply_succAbove, Polynomial.map_C, eval_X, Polynomial.eval_C, implies_true, - and_self] + AlgHom.toRingHom_eq_coe, coe_aeval_eq_eval, AlgHom.coe_comp, + Polynomial.coe_aeval_eq_eval, AlgHom.coe_mk, coe_mapRingHom, comp_apply, + AlgEquiv.toAlgHom_apply, finSuccEquivNth_apply, eval₂Hom_X', Polynomial.map_X, + Polynomial.eval_X, Fin.insertNth_apply_succAbove, Polynomial.map_C, eval_X, + Polynomial.eval_C, implies_true, and_self] /-- A monomial index `m` is in the support of the `i`-th coefficient of `finSuccEquivNth R p f` if and only if `m.insertNth p i` is in the support of `f`. -/ diff --git a/CompPoly/ToMathlib/Polynomial/BivariateDegree.lean b/CompPoly/ToMathlib/Polynomial/BivariateDegree.lean index 09c8aae2..7e6526df 100644 --- a/CompPoly/ToMathlib/Polynomial/BivariateDegree.lean +++ b/CompPoly/ToMathlib/Polynomial/BivariateDegree.lean @@ -321,7 +321,7 @@ theorem natDegree_sum_eq_of_unique {α : Type*} {s : Finset α} {f : α → F[X] intro y hy hym exact hcoeff_others y hy hym have hcoeff_eq : (∑ x ∈ s, f x).coeff deg = (f mx).coeff deg := by - rw [Polynomial.finset_sum_coeff (s := s) (f := f) (n := deg)] + rw [Polynomial.finsetSum_coeff (s := s) (f := f) (n := deg)] exact hsum_coeff have hcoeff_ne0 : (∑ x ∈ s, f x).coeff deg ≠ 0 := by simpa [hcoeff_eq] using hmx_coeff_ne0 diff --git a/CompPoly/Univariate/Basic.lean b/CompPoly/Univariate/Basic.lean index b1b5d471..93446282 100644 --- a/CompPoly/Univariate/Basic.lean +++ b/CompPoly/Univariate/Basic.lean @@ -173,6 +173,19 @@ instance [Semiring R] [BEq R] [LawfulBEq R] [Nontrivial R] : Nontrivial (CPolyno @[reducible] def coeff [Zero R] (p : CPolynomial R) (i : ℕ) : R := p.val.coeff i +/-- Build a canonical polynomial from a dense coefficient array. -/ +def ofArray [Zero R] [BEq R] [LawfulBEq R] + (coeffs : Array R) : CPolynomial R := + let raw : CPolynomial.Raw R := coeffs + ⟨raw.trim, Trim.isCanonical_trim raw⟩ + +/-- Coefficients of `ofArray` are exactly the source array entries, with zero default. -/ +theorem coeff_ofArray [Zero R] [BEq R] [LawfulBEq R] + (coeffs : Array R) (i : Nat) : + (CPolynomial.ofArray coeffs).coeff i = coeffs.getD i 0 := by + unfold CPolynomial.coeff CPolynomial.ofArray + rw [CPolynomial.Raw.Trim.coeff_eq_coeff] + /-- The constant polynomial `C r`. -/ def C [Zero R] [BEq R] [LawfulBEq R] (r : R) : CPolynomial R := ⟨(Raw.C r).trim, Trim.isCanonical_trim (Raw.C r)⟩ @@ -181,6 +194,10 @@ def C [Zero R] [BEq R] [LawfulBEq R] (r : R) : CPolynomial R := def X [Semiring R] [BEq R] [LawfulBEq R] [Nontrivial R] : CPolynomial R := ⟨Raw.X, Trim.isCanonical_of_trim_eq X_canonical⟩ +/-- The monic linear factor `X - C x`. -/ +def linearFactor [Field R] [BEq R] [LawfulBEq R] (x : R) : CPolynomial R := + (C (-x) : CPolynomial R) + (X : CPolynomial R) + /-- Construct a canonical monomial `c * X^n` as a `CPolynomial R`. The result is canonical (no trailing zeros) when `c ≠ 0`. @@ -842,6 +859,21 @@ theorem mod_zero [Field R] [BEq R] [LawfulBEq R] (p : CPolynomial R) : p.mod 0 = instance [Field R] [BEq R] [LawfulBEq R] : Div (CPolynomial R) := ⟨div⟩ instance [Field R] [BEq R] [LawfulBEq R] : Mod (CPolynomial R) := ⟨mod⟩ +/-- Normalize a nonzero polynomial to monic form. The zero polynomial stays zero. -/ +def monicNormalize [Field R] [BEq R] [LawfulBEq R] + (p : CPolynomial R) : CPolynomial R := + CPolynomial.ofArray (Raw.monicNormalize p.val) + +/-- Euclidean gcd with explicit fuel, normalized to a monic result. -/ +def gcdMonicWithFuel [Field R] [BEq R] [LawfulBEq R] : + Nat → CPolynomial R → CPolynomial R → CPolynomial R + | fuel, p, q => CPolynomial.ofArray (Raw.gcdMonicWithFuel fuel p.val q.val) + +/-- Monic Euclidean gcd for canonical univariate polynomials. -/ +def gcdMonic [Field R] [BEq R] [LawfulBEq R] + (p q : CPolynomial R) : CPolynomial R := + CPolynomial.ofArray (Raw.gcdMonic p.val q.val) + end Division end Operations @@ -1047,8 +1079,8 @@ instance [Ring R] [BEq R] [LawfulBEq R] : Sub (CPolynomial R) where /-- Erase the coefficient at index `n` (same as `p` except `coeff n = 0`, then trimmed). - Uses an in-place `Array.setIfInBounds` rather than subtracting a monomial, avoiding - the allocation of a length-`n` monomial array plus the padding/zip passes of `sub`. -/ + Uses an in-place `Array.setIfInBounds` and avoids allocating a length-`n` monomial + array plus the padding/zip passes of `sub`. -/ def erase [Zero R] [BEq R] [LawfulBEq R] (n : ℕ) (p : CPolynomial R) : CPolynomial R := let arr : CPolynomial.Raw R := p.val.setIfInBounds n 0 ⟨arr.trim, Trim.isCanonical_trim arr⟩ diff --git a/CompPoly/Univariate/BatchEval.lean b/CompPoly/Univariate/BatchEval.lean index b79d287d..b07b2cb3 100644 --- a/CompPoly/Univariate/BatchEval.lean +++ b/CompPoly/Univariate/BatchEval.lean @@ -3,10 +3,11 @@ Copyright (c) 2026 CompPoly. All rights reserved. Released under Apache 2.0 license as described in the file LICENSE. Authors: Valerii Huhnin -/ -import CompPoly.Univariate.BatchEval.Context import CompPoly.Univariate.BatchEval.Naive import CompPoly.Univariate.BatchEval.SubproductTree import CompPoly.Univariate.BatchEval.Correctness +import CompPoly.Univariate.BatchEval.Context +import CompPoly.Univariate.Context /-! # Univariate Batch Evaluation diff --git a/CompPoly/Univariate/BatchEval/Context.lean b/CompPoly/Univariate/BatchEval/Context.lean index 38d083f5..072193c3 100644 --- a/CompPoly/Univariate/BatchEval/Context.lean +++ b/CompPoly/Univariate/BatchEval/Context.lean @@ -3,91 +3,44 @@ Copyright (c) 2026 CompPoly. All rights reserved. Released under Apache 2.0 license as described in the file LICENSE. Authors: Valerii Huhnin -/ -import CompPoly.Univariate.DivisionCorrectness -import CompPoly.Univariate.NTT.FastMul -import CompPoly.Univariate.NTTFast.Correctness.Pipeline + +import CompPoly.Univariate.BatchEval.Correctness /-! -# Batch Evaluation Contexts +# Batch-Evaluation Contexts -Algorithm dictionaries for univariate batch-evaluation implementations. +Explicit context wrappers for reusable univariate batch-evaluation backends. -/ namespace CompPoly -namespace CPolynomial - -variable {R : Type*} - -/-- Explicit multiplication backend for batch-evaluation algorithms. -/ -structure MulContext (R : Type*) [Semiring R] [BEq R] [LawfulBEq R] where - /-- Multiply two canonical polynomials. -/ - mul : CPolynomial R → CPolynomial R → CPolynomial R - /-- The backend agrees with canonical polynomial multiplication. -/ - mul_eq_mul : ∀ p q, mul p q = p * q - -/-- Explicit remainder backend for algorithms that only need reduction modulo monic divisors. -/ -structure ModContext (R : Type*) [Field R] [BEq R] [LawfulBEq R] where - /-- Reduce the first polynomial modulo the second, assuming the divisor is monic. -/ - modByMonic : CPolynomial R → CPolynomial R → CPolynomial R - /-- The backend agrees with the canonical monic-remainder operation. -/ - modByMonic_eq_modByMonic : ∀ p q, modByMonic p q = CPolynomial.modByMonic p q - -namespace MulContext - -/-- The default multiplication context, backed by canonical `CPolynomial` multiplication. -/ -def naive [Semiring R] [BEq R] [LawfulBEq R] : MulContext R where - mul p q := p * q - mul_eq_mul _ _ := rfl - -/-- -NTT-backed multiplication context with canonical multiplication as a fallback. -The context asks the selector for a domain that fits the current operands. If no -supported domain is available, it falls back to ordinary `CPolynomial` -multiplication. --/ -def ntt [Field R] [BEq R] [LawfulBEq R] - (bestDomainForLength? : (requiredLen : Nat) → - Option (NTT.FittingDomain R requiredLen)) : - MulContext R where - mul := NTT.FastMul.withFallback bestDomainForLength? - mul_eq_mul := NTT.FastMul.withFallback_eq_mul bestDomainForLength? - -/-- -NTTFast-backed multiplication context with canonical multiplication as a fallback. - -The context asks the selector for a domain that fits the current operands. If no -supported domain is available, it falls back to ordinary `CPolynomial` -multiplication. --/ -def nttFast [Field R] [BEq R] [LawfulBEq R] - (bestDomainForLength? : (requiredLen : Nat) → - Option (NTT.FittingDomain R requiredLen)) : - MulContext R where - mul := NTTFast.withFallback bestDomainForLength? - mul_eq_mul := NTTFast.withFallback_eq_mul bestDomainForLength? - -end MulContext +namespace CPolynomial -namespace ModContext +/-- Contract wrapper for batch evaluation backends. -/ +structure BatchEvalContext (R : Type*) [Semiring R] where + evalBatchWith : CPolynomial R → Array R → Array R + correct : ∀ p xs, evalBatchWith p xs = CPolynomial.evalBatch p xs -/-- The default monic-remainder context, backed by `CPolynomial.modByMonic`. -/ -def naive [Field R] [BEq R] [LawfulBEq R] : ModContext R where - modByMonic p q := CPolynomial.modByMonic p q - modByMonic_eq_modByMonic _ _ := rfl +namespace BatchEvalContext -/-- A remainder-only backend for monic remainders. -/ -def remainderOnly [Field R] [BEq R] [LawfulBEq R] : ModContext R where - modByMonic p q := CPolynomial.modByMonicRemainderOnly p q - modByMonic_eq_modByMonic p q := CPolynomial.modByMonicRemainderOnly_eq_modByMonic p q +/-- Horner-backed batch evaluation context. -/ +def horner (R : Type*) [Semiring R] : BatchEvalContext R where + evalBatchWith := CPolynomial.evalBatchHorner + correct := by + intro p xs + exact CPolynomial.evalBatchHorner_eq_evalBatch p xs -/-- A reversal-based monic-remainder backend parameterized by low-product multiplication. -/ -def reversal [Field R] [BEq R] [LawfulBEq R] - (M : Raw.MulLowContext R) : ModContext R where - modByMonic p q := CPolynomial.modByMonicByReversal M p q - modByMonic_eq_modByMonic p q := CPolynomial.modByMonicByReversal_eq_modByMonic M p q +/-- Subproduct-tree-backed batch evaluation context. -/ +def subproduct (R : Type*) [Field R] [BEq R] [LawfulBEq R] + (M : CPolynomial.MulContext R) (D : CPolynomial.ModContext R) : + BatchEvalContext R where + evalBatchWith := CPolynomial.evalBatchSubproduct M D + correct := by + intro p xs + exact CPolynomial.evalBatchSubproduct_eq_evalBatch M D p xs -end ModContext +end BatchEvalContext end CPolynomial + end CompPoly diff --git a/CompPoly/Univariate/BatchEval/Correctness.lean b/CompPoly/Univariate/BatchEval/Correctness.lean index 51267229..c84bc85b 100644 --- a/CompPoly/Univariate/BatchEval/Correctness.lean +++ b/CompPoly/Univariate/BatchEval/Correctness.lean @@ -20,7 +20,7 @@ variable {R : Type*} private theorem eval_linearFactor [Field R] [BEq R] [LawfulBEq R] (x : R) : (linearFactor x).eval x = 0 := by rw [eval_toPoly] - simp [linearFactor, toPoly_sub, X_toPoly, C_toPoly] + simp [linearFactor, toPoly_add, X_toPoly, C_toPoly] private theorem eval_modContext_eq_self_of_eval_eq_zero [Field R] [BEq R] [LawfulBEq R] (D : ModContext R) (p q : CPolynomial R) diff --git a/CompPoly/Univariate/BatchEval/SubproductTree.lean b/CompPoly/Univariate/BatchEval/SubproductTree.lean index 5052ec3f..3db4df70 100644 --- a/CompPoly/Univariate/BatchEval/SubproductTree.lean +++ b/CompPoly/Univariate/BatchEval/SubproductTree.lean @@ -3,8 +3,8 @@ Copyright (c) 2026 CompPoly. All rights reserved. Released under Apache 2.0 license as described in the file LICENSE. Authors: Valerii Huhnin -/ -import CompPoly.Univariate.BatchEval.Context import CompPoly.Univariate.BatchEval.Naive +import CompPoly.Univariate.Context /-! # Subproduct-Tree Batch Evaluation @@ -18,10 +18,6 @@ namespace CPolynomial variable {R : Type*} -/-- The linear factor `X - C x` used as a subproduct-tree leaf. -/ -def linearFactor [Field R] [BEq R] [LawfulBEq R] (x : R) : CPolynomial R := - X - C x - /-- A subproduct tree stores the product polynomial at every node. -/ inductive SubproductTree (R : Type*) [Zero R] where /-- A leaf for one evaluation point. The polynomial is `X - C x`. -/ diff --git a/CompPoly/Univariate/Context.lean b/CompPoly/Univariate/Context.lean new file mode 100644 index 00000000..028812b1 --- /dev/null +++ b/CompPoly/Univariate/Context.lean @@ -0,0 +1,109 @@ +/- +Copyright (c) 2026 CompPoly. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ +import CompPoly.Univariate.DivisionCorrectness +import CompPoly.Univariate.NTT.FastMul +import CompPoly.Univariate.NTT.FastMulLow +import CompPoly.Univariate.NTTFast.Correctness.Pipeline +import CompPoly.Univariate.NTTFast.FastMulLow +import CompPoly.Univariate.Raw.Context + +/-! +# Univariate Algorithm Contexts + +Algorithm dictionaries for reusable univariate polynomial operations, including +canonical, NTT, and NTTFast-backed implementations. +-/ + +namespace CompPoly +namespace CPolynomial + +variable {R : Type*} + +/-- Explicit multiplication backend for univariate polynomial algorithms. -/ +structure MulContext (R : Type*) [Semiring R] [BEq R] [LawfulBEq R] where + /-- Multiply two canonical polynomials. -/ + mul : CPolynomial R → CPolynomial R → CPolynomial R + /-- The backend agrees with canonical polynomial multiplication. -/ + mul_eq_mul : ∀ p q, mul p q = p * q + +/-- Explicit remainder backend for algorithms that only need reduction modulo monic divisors. -/ +structure ModContext (R : Type*) [Field R] [BEq R] [LawfulBEq R] where + /-- Reduce the first polynomial modulo the second, assuming the divisor is monic. -/ + modByMonic : CPolynomial R → CPolynomial R → CPolynomial R + /-- The backend agrees with the canonical monic-remainder operation. -/ + modByMonic_eq_modByMonic : ∀ p q, modByMonic p q = CPolynomial.modByMonic p q + +namespace MulContext + +/-- The default multiplication context, backed by canonical `CPolynomial` multiplication. -/ +def naive [Semiring R] [BEq R] [LawfulBEq R] : MulContext R where + mul p q := p * q + mul_eq_mul _ _ := rfl + +/-- +NTT-backed multiplication context with canonical multiplication for unsupported lengths. + +The context asks the selector for a domain that fits the current operands. If no +supported domain is available, it uses ordinary `CPolynomial` multiplication. +-/ +def ntt [Field R] [BEq R] [LawfulBEq R] + (bestDomainForLength? : (requiredLen : Nat) → + Option (NTT.FittingDomain R requiredLen)) : + MulContext R where + mul := NTT.FastMul.withFallback bestDomainForLength? + mul_eq_mul := NTT.FastMul.withFallback_eq_mul bestDomainForLength? + +/-- +NTTFast-backed multiplication context with canonical multiplication for unsupported lengths. + +The context asks the selector for a domain that fits the current operands. If no +supported domain is available, it uses ordinary `CPolynomial` multiplication. +-/ +def nttFast [Field R] [BEq R] [LawfulBEq R] + (bestDomainForLength? : (requiredLen : Nat) → + Option (NTT.FittingDomain R requiredLen)) : + MulContext R where + mul := NTTFast.withFallback bestDomainForLength? + mul_eq_mul := NTTFast.withFallback_eq_mul bestDomainForLength? + +end MulContext + +namespace ModContext + +/-- The default monic-remainder context, backed by `CPolynomial.modByMonic`. -/ +def naive [Field R] [BEq R] [LawfulBEq R] : ModContext R where + modByMonic p q := CPolynomial.modByMonic p q + modByMonic_eq_modByMonic _ _ := rfl + +/-- A remainder-only backend for monic remainders. -/ +def remainderOnly [Field R] [BEq R] [LawfulBEq R] : ModContext R where + modByMonic p q := CPolynomial.modByMonicRemainderOnly p q + modByMonic_eq_modByMonic p q := CPolynomial.modByMonicRemainderOnly_eq_modByMonic p q + +/-- A reversal-based monic-remainder backend parameterized by low-product multiplication. -/ +def reversal [Field R] [BEq R] [LawfulBEq R] + (M : Raw.MulLowContext R) : ModContext R where + modByMonic p q := CPolynomial.modByMonicByReversal M p q + modByMonic_eq_modByMonic p q := CPolynomial.modByMonicByReversal_eq_modByMonic M p q + +/-- Monic remainders by reversal, using an NTT low-product backend. -/ +def reversalNtt [Field R] [BEq R] [LawfulBEq R] + (bestDomainForLength? : (requiredLen : Nat) → + Option (NTT.FittingDomain R requiredLen)) : + ModContext R := + reversal (NTT.FastMulLow.withFallback bestDomainForLength?) + +/-- Monic remainders by reversal, using an NTTFast low-product backend. -/ +def reversalNttFast [Field R] [BEq R] [LawfulBEq R] + (bestDomainForLength? : (requiredLen : Nat) → + Option (NTT.FittingDomain R requiredLen)) : + ModContext R := + reversal (NTTFast.FastMulLow.withFallback bestDomainForLength?) + +end ModContext + +end CPolynomial +end CompPoly diff --git a/CompPoly/Univariate/EuclideanAlgorithm.lean b/CompPoly/Univariate/EuclideanAlgorithm.lean index 4a89c633..2844c2a8 100644 --- a/CompPoly/Univariate/EuclideanAlgorithm.lean +++ b/CompPoly/Univariate/EuclideanAlgorithm.lean @@ -1,7 +1,7 @@ /- Copyright (c) 2026 CompPoly. All rights reserved. Released under Apache 2.0 license as described in the file LICENSE. -Authors: Juan Conejero +Authors: Juan Conejero, Valerii Huhnin -/ import CompPoly.Univariate.Basic import CompPoly.Univariate.DivisionCorrectness @@ -175,6 +175,24 @@ def normXgcd [Field R] [BEq R] [LawfulBEq R] let c := res.1.leadingCoeff⁻¹ (c • res.1, c • res.2.1, c • res.2.2) +/-- The normalized extended-gcd output satisfies the Bezout identity. -/ +theorem normXgcd_bezout [Field R] [BEq R] [LawfulBEq R] + (p q : CPolynomial R) (threshold : ℕ) : + Bezout p q (normXgcd p q threshold) := by + unfold normXgcd + have h := xgcd_bezout p q threshold + simp only [Bezout] at h ⊢ + rw [h] + apply toPolyLinearEquiv.injective + change (((p.xgcd q threshold).2.1 * p + (p.xgcd q threshold).2.2 * q).leadingCoeff⁻¹ • + ((p.xgcd q threshold).2.1 * p + (p.xgcd q threshold).2.2 * q)).toPoly = + (((p.xgcd q threshold).2.1 * p + (p.xgcd q threshold).2.2 * q).leadingCoeff⁻¹ • + (p.xgcd q threshold).2.1 * p + + ((p.xgcd q threshold).2.1 * p + (p.xgcd q threshold).2.2 * q).leadingCoeff⁻¹ • + (p.xgcd q threshold).2.2 * q).toPoly + simp only [toPoly_smul, toPoly_add, toPoly_mul, Polynomial.smul_eq_C_mul] + ring + /-- The gcd component of CompPoly's `normXgcd` is the normalization of Mathlib's `EuclideanDomain.gcd` -/ theorem normXgcd_fst_toPoly @@ -187,6 +205,135 @@ theorem normXgcd_fst_toPoly rw [Polynomial.smul_eq_C_mul, normalize_apply, Polynomial.coe_normUnit_of_ne_zero h, _root_.mul_comm] +private theorem Raw.toPoly_smul [Semiring R] [BEq R] [LawfulBEq R] + (c : R) (p : CPolynomial.Raw R) : + (c • p).toPoly = c • p.toPoly := by + ext i + rw [Polynomial.coeff_smul] + rw [CPolynomial.Raw.coeff_toPoly, CPolynomial.Raw.coeff_toPoly] + exact CPolynomial.Raw.smul_coeff c p i + +/-- Monic normalization of computable polynomials agrees with Mathlib normalization. -/ +theorem monicNormalize_toPoly_eq_normalize + [Field R] [BEq R] [LawfulBEq R] [DecidableEq R] + (p : CPolynomial R) : + (CPolynomial.monicNormalize p).toPoly = normalize p.toPoly := by + unfold CPolynomial.monicNormalize CPolynomial.Raw.monicNormalize + rw [ofArray_toPoly, CPolynomial.trim_eq] + by_cases hpraw : ((p.val : CPolynomial.Raw R) == 0) + · have hp : p = 0 := CPolynomial.ext (LawfulBEq.eq_of_beq hpraw) + rw [if_pos hpraw, hp] + rw [CPolynomial.Raw.toPoly_zero, CPolynomial.toPoly_zero, normalize_zero] + · have hp : p ≠ 0 := by + intro hp + exact hpraw (by subst p; rfl) + rw [if_neg hpraw, Raw.toPoly_smul] + have hlead : CPolynomial.Raw.leadingCoeff p.val = p.leadingCoeff := by + simp [CPolynomial.Raw.leadingCoeff, CPolynomial.leadingCoeff, CPolynomial.trim_eq] + rw [hlead] + have hpoly : p.toPoly ≠ 0 := (toPoly_eq_zero_iff p).not.mpr hp + rw [Polynomial.smul_eq_C_mul, normalize_apply, + Polynomial.coe_normUnit_of_ne_zero hpoly, _root_.mul_comm] + rw [CPolynomial.leadingCoeff_toPoly] + change p.toPoly * Polynomial.C p.toPoly.leadingCoeff⁻¹ = + p.toPoly * Polynomial.C p.toPoly.leadingCoeff⁻¹ + rfl + +private theorem gcdMonicWithFuel_toPoly_eq_normalize_gcd + [Field R] [BEq R] [LawfulBEq R] [DecidableEq R] + (fuel : ℕ) (p q : CPolynomial R) + (hfuel : q.toPoly.degree < fuel) : + (CPolynomial.gcdMonicWithFuel fuel p q).toPoly = + normalize (EuclideanDomain.gcd p.toPoly q.toPoly) := by + induction fuel generalizing p q with + | zero => + have hqpoly : q.toPoly = 0 := + Polynomial.degree_eq_bot.mp (Nat.WithBot.lt_zero_iff.mp hfuel) + have hq : q = 0 := (toPoly_eq_zero_iff q).mp hqpoly + subst q + rw [CPolynomial.toPoly_zero] + change (CPolynomial.monicNormalize p).toPoly = + normalize (EuclideanDomain.gcd p.toPoly (0 : Polynomial R)) + rw [monicNormalize_toPoly_eq_normalize, EuclideanDomain.gcd_zero_right] + | succ fuel ih => + by_cases hq : q = 0 + · subst q + simp [CPolynomial.gcdMonicWithFuel, CPolynomial.Raw.gcdMonicWithFuel, + CPolynomial.trim_eq, CPolynomial.toPoly_zero] + have hzero : (↑(0 : CPolynomial R) : CPolynomial.Raw R) = (#[] : CPolynomial.Raw R) := rfl + rw [if_pos hzero] + change (CPolynomial.monicNormalize p).toPoly = normalize p.toPoly + exact monicNormalize_toPoly_eq_normalize p + · have hqraw : ¬((q.val : CPolynomial.Raw R) == 0) := by + intro h + exact hq (CPolynomial.ext (LawfulBEq.eq_of_beq h)) + rw [CPolynomial.gcdMonicWithFuel, CPolynomial.Raw.gcdMonicWithFuel, + CPolynomial.trim_eq, CPolynomial.trim_eq, if_neg hqraw] + change (CPolynomial.gcdMonicWithFuel fuel q (p % q)).toPoly = + normalize (EuclideanDomain.gcd p.toPoly q.toPoly) + rw [ih] + · rw [show (p % q).toPoly = q.leadingCoeff⁻¹ • (p.toPoly % q.toPoly) by + exact mod_toPoly_eq_smul_mod p q] + have hunit : IsUnit (Polynomial.C q.leadingCoeff⁻¹ : Polynomial R) := by + exact Polynomial.isUnit_C.mpr + (isUnit_iff_ne_zero.mpr (inv_ne_zero (CPolynomial.leadingCoeff_ne_zero hq))) + have hassoc : + Associated (q.leadingCoeff⁻¹ • (p.toPoly % q.toPoly)) + (p.toPoly % q.toPoly) := by + simpa [Polynomial.smul_eq_C_mul] using + associated_unit_mul_left (p.toPoly % q.toPoly) + (Polynomial.C q.leadingCoeff⁻¹) hunit + refine normalize_eq_normalize_iff_associated.mpr + (associated_of_dvd_dvd ?_ ?_) + · refine EuclideanDomain.dvd_gcd ?_ (EuclideanDomain.gcd_dvd_left _ _) + exact (EuclideanDomain.dvd_mod_iff (EuclideanDomain.gcd_dvd_left + q.toPoly (q.leadingCoeff⁻¹ • (p.toPoly % q.toPoly)))).mp + ((EuclideanDomain.gcd_dvd_right q.toPoly + (q.leadingCoeff⁻¹ • (p.toPoly % q.toPoly))).trans hassoc.dvd) + · refine EuclideanDomain.dvd_gcd (EuclideanDomain.gcd_dvd_right _ _) ?_ + exact ((EuclideanDomain.dvd_mod_iff + (EuclideanDomain.gcd_dvd_right p.toPoly q.toPoly)).mpr + (EuclideanDomain.gcd_dvd_left p.toPoly q.toPoly)).trans hassoc.symm.dvd + · have hqpoly : q.toPoly ≠ 0 := (toPoly_eq_zero_iff q).not.mpr hq + have hmod : + (p % q).toPoly.degree ≤ (p.toPoly % q.toPoly).degree := by + rw [show (p % q).toPoly = q.leadingCoeff⁻¹ • (p.toPoly % q.toPoly) by + exact mod_toPoly_eq_smul_mod p q] + exact Polynomial.degree_smul_le _ _ + exact lt_of_le_of_lt hmod + (lt_of_lt_of_le (Polynomial.degree_mod_lt _ hqpoly) + (Order.le_of_lt_succ hfuel)) + +/-- The specialized monic gcd has the normalized Mathlib gcd as its `toPoly` +image. -/ +theorem gcdMonic_toPoly_eq_normalize_gcd + [Field R] [BEq R] [LawfulBEq R] [DecidableEq R] + (p q : CPolynomial R) : + (CPolynomial.gcdMonic p q).toPoly = + normalize (EuclideanDomain.gcd p.toPoly q.toPoly) := by + simpa [CPolynomial.gcdMonic, CPolynomial.Raw.gcdMonic] using + gcdMonicWithFuel_toPoly_eq_normalize_gcd + (p.val.size + q.val.size + 1) p q (by + have hqdeg : q.toPoly.degree < q.val.size := by + rw [← degree_toPoly] + exact mem_degreeLT_iff_size_le.mpr le_rfl + have hleNat : q.val.size ≤ p.val.size + q.val.size + 1 := by + omega + have hle : (q.val.size : WithBot ℕ) ≤ + (p.val.size + q.val.size + 1 : ℕ) := + WithBot.coe_le_coe.2 hleNat + exact lt_of_lt_of_le hqdeg hle) + +/-- The specialized monic gcd agrees with the gcd component of normalized +extended gcd. -/ +theorem gcdMonic_eq_normXgcd_fst + [Field R] [BEq R] [LawfulBEq R] [DecidableEq R] + (p q : CPolynomial R) : + CPolynomial.gcdMonic p q = (CPolynomial.normXgcd p q).1 := by + apply toPolyLinearEquiv.injective + change (CPolynomial.gcdMonic p q).toPoly = (CPolynomial.normXgcd p q).1.toPoly + rw [gcdMonic_toPoly_eq_normalize_gcd, normXgcd_fst_toPoly] + /-- The Bezout component of `normXgcd` under `toPoly` is Mathlib's `EuclideanDomain.xgcd` scaled by the inverse leading coefficient of the gcd -/ diff --git a/CompPoly/Univariate/ManyEval/Correctness.lean b/CompPoly/Univariate/ManyEval/Correctness.lean index 31c5a7b9..3317ba59 100644 --- a/CompPoly/Univariate/ManyEval/Correctness.lean +++ b/CompPoly/Univariate/ManyEval/Correctness.lean @@ -4,6 +4,7 @@ Released under Apache 2.0 license as described in the file LICENSE. Authors: Valerii Huhnin -/ import CompPoly.Univariate.ManyEval.Basic +import CompPoly.Data.Array.Lemmas /-! # Many-Polynomial Evaluation Correctness @@ -112,18 +113,6 @@ private lemma powerTable_getD_eq_pow [Semiring R] (x : R) {limit i : Nat} (by simp) (by intro i hi; omega) (by simp) i (by simpa using hi) simpa using h -private lemma list_foldl_congr {α β : Type*} {f g : α → β → α} - (xs : List β) (acc : α) (h : ∀ acc' x, x ∈ xs → f acc' x = g acc' x) : - xs.foldl f acc = xs.foldl g acc := by - induction xs generalizing acc with - | nil => simp - | cons x xs ih => - simp only [List.foldl_cons] - rw [h acc x (by simp)] - apply ih - intro acc' y hy - exact h acc' y (by simp [hy]) - private lemma evalWithPowersLoop_eq_foldl_range [Semiring R] (coeffs powers : Array R) (limit : Nat) (hcoeffs : limit ≤ coeffs.size) (hpowers : limit ≤ powers.size) (i : Nat) (acc : R) : @@ -160,7 +149,7 @@ private lemma foldl_zipIdx_eq_range'_getD_aux [Semiring R] rw [List.range'_succ, List.foldl_cons] simp only [Nat.sub_self, List.getD_cons_zero] rw [ih] - apply list_foldl_congr + apply List.foldl_congr_of_mem intro _acc' j hj have hsub : j - offset = (j - (offset + 1)) + 1 := by have hjmem := List.mem_range'.mp hj @@ -172,12 +161,6 @@ private lemma array_getD_toList {α : Type*} (a : Array α) (i : Nat) (d : α) : cases a simp -private lemma array_foldl_zipIdx_eq_list [Semiring R] (coeffs : Array R) (x acc : R) : - Array.foldl (fun acc ai ↦ acc + ai.1 * x ^ ai.2) acc coeffs.zipIdx 0 coeffs.size = - List.foldl (fun acc ai ↦ acc + ai.1 * x ^ ai.2) acc coeffs.toList.zipIdx := by - cases coeffs - simp - private lemma evalWithPowers_eq_eval [Semiring R] (coeffs powers : Array R) (x : R) (hsize : coeffs.size ≤ powers.size) (hpowers : ∀ i, i < coeffs.size → powers.getD i 0 = x ^ i) : @@ -185,9 +168,9 @@ private lemma evalWithPowers_eq_eval [Semiring R] (coeffs powers : Array R) (x : coeffs.zipIdx.foldl (fun acc ai ↦ acc + ai.1 * x ^ ai.2) 0 := by simp [ManyEval.evalWithPowers, Nat.min_eq_left hsize] rw [evalWithPowersLoop_eq_foldl_range] - rw [array_foldl_zipIdx_eq_list] + rw [Array.foldl_zipIdx_eq_foldl_toList_zipIdx_size] rw [foldl_zipIdx_eq_range'_getD_aux] - apply list_foldl_congr + apply List.foldl_congr_of_mem intro _acc' i hi rw [← array_getD_toList] rw [hpowers i] diff --git a/CompPoly/Univariate/Modular.lean b/CompPoly/Univariate/Modular.lean new file mode 100644 index 00000000..5af28b95 --- /dev/null +++ b/CompPoly/Univariate/Modular.lean @@ -0,0 +1,76 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Context + +/-! +# Modular Operations on Univariate Polynomials + +Reusable executable modular arithmetic for canonical `CPolynomial`s over public +`CPolynomial.MulContext` and `CPolynomial.ModContext` backends. +-/ + +namespace CompPoly + +namespace CPolynomial + +/-- Multiply modulo a polynomial, treating zero modulus as no reduction. -/ +def mulModWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : MulContext F) (D : ModContext F) + (modulus p q : CPolynomial F) : CPolynomial F := + let product := M.mul p q + if modulus == 0 then + product + else + D.modByMonic product (monicNormalize modulus) + +/-- Binary modular exponentiation accumulator. -/ +def powModBinaryAuxWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : MulContext F) (D : ModContext F) + (modulus : CPolynomial F) : + Nat → CPolynomial F → CPolynomial F → CPolynomial F + | 0, acc, _ => acc + | n + 1, acc, current => + let acc' := + if (n + 1) % 2 == 1 then + mulModWith M D modulus acc current + else + acc + let current' := mulModWith M D modulus current current + powModBinaryAuxWith M D modulus ((n + 1) / 2) acc' current' +termination_by exp _ _ => exp +decreasing_by + exact Nat.div_lt_self (Nat.succ_pos n) (by decide) + +/-- Modular exponentiation by repeated squaring. -/ +def powModWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : MulContext F) (D : ModContext F) + (modulus base : CPolynomial F) (exponent : Nat) : CPolynomial F := + let oneMod := + if modulus == 0 then + (1 : CPolynomial F) + else + D.modByMonic (1 : CPolynomial F) (monicNormalize modulus) + powModBinaryAuxWith M D modulus exponent oneMod base + +/-- `X mod modulus`, with zero modulus treated as no reduction. -/ +def xModWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (D : ModContext F) : CPolynomial F → CPolynomial F + | modulus => + if modulus == 0 then + (CPolynomial.X : CPolynomial F) + else + D.modByMonic (CPolynomial.X : CPolynomial F) (monicNormalize modulus) + +/-- `(X^q mod modulus) - (X mod modulus)`, without materializing `X^q - X`. -/ +def xPowSubXModWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : MulContext F) (D : ModContext F) + (q : Nat) (modulus : CPolynomial F) : CPolynomial F := + powModWith M D modulus (CPolynomial.X : CPolynomial F) q - xModWith D modulus + +end CPolynomial + +end CompPoly diff --git a/CompPoly/Univariate/NTT/KoalaBear.lean b/CompPoly/Univariate/NTT/KoalaBear.lean index a721a202..aeba091b 100644 --- a/CompPoly/Univariate/NTT/KoalaBear.lean +++ b/CompPoly/Univariate/NTT/KoalaBear.lean @@ -49,6 +49,55 @@ def domainOfLogN (logN : Nat) (hlogN : logN ≤ KoalaBear.twoAdicity) : KoalaBear.isPrimitiveRoot_twoAdicGenerator (bitsOfLogN logN hlogN) natCast_ne_zero := twoPowNatCast_ne_zero logN hlogN +/-- KoalaBear NTT domain lookup for dynamic multiplication contexts. -/ +def bestDomainForLength? (requiredLen : Nat) : + Option (FittingDomain KoalaBear.Field requiredLen) := + CPolynomial.NTT.bestDomainForLength? KoalaBear.twoAdicity + domainOfLogN (by intro _ _; rfl) requiredLen + +/-- Fast KoalaBear two-adic generators are primitive roots of the same orders. -/ +theorem fast_isPrimitiveRoot_twoAdicGenerator + (bits : Fin (KoalaBear.twoAdicity + 1)) : + IsPrimitiveRoot (KoalaBear.Fast.ofField KoalaBear.twoAdicGenerators[bits]) + (2 ^ (bits : Nat)) := by + have hbasic : IsPrimitiveRoot KoalaBear.twoAdicGenerators[bits] (2 ^ (bits : Nat)) := + KoalaBear.isPrimitiveRoot_twoAdicGenerator bits + have hfast : + IsPrimitiveRoot (KoalaBear.Fast.ringEquiv.symm KoalaBear.twoAdicGenerators[bits]) + (2 ^ (bits : Nat)) := by + exact hbasic.map_of_injective KoalaBear.Fast.ringEquiv.symm.injective + simpa using hfast + +/-- The fast KoalaBear NTT domain size is nonzero for supported two-adic sizes. -/ +theorem fast_twoPowNatCast_ne_zero + (logN : Nat) (hlogN : logN ≤ KoalaBear.twoAdicity) : + (((2 ^ logN : Nat) : KoalaBear.Fast.Field) ≠ 0) := by + intro hzero + exact twoPowNatCast_ne_zero logN hlogN (by + calc + (((2 ^ logN : Nat) : KoalaBear.Field)) = + KoalaBear.Fast.toField (((2 ^ logN : Nat) : KoalaBear.Fast.Field)) := by + rw [KoalaBear.Fast.toField_natCast] + _ = KoalaBear.Fast.toField 0 := congrArg KoalaBear.Fast.toField hzero + _ = 0 := KoalaBear.Fast.toField_zero) + +/-- Fast KoalaBear radix-2 NTT domain for a supported two-adic size. -/ +def fastDomainOfLogN (logN : Nat) (hlogN : logN ≤ KoalaBear.twoAdicity) : + Domain KoalaBear.Fast.Field where + logN := logN + omega := KoalaBear.Fast.ofField + KoalaBear.twoAdicGenerators[bitsOfLogN logN hlogN] + primitive := by + simpa [bitsOfLogN] using + fast_isPrimitiveRoot_twoAdicGenerator (bitsOfLogN logN hlogN) + natCast_ne_zero := fast_twoPowNatCast_ne_zero logN hlogN + +/-- Fast KoalaBear NTT domain lookup for dynamic multiplication contexts. -/ +def fastBestDomainForLength? (requiredLen : Nat) : + Option (FittingDomain KoalaBear.Fast.Field requiredLen) := + CPolynomial.NTT.bestDomainForLength? KoalaBear.twoAdicity + fastDomainOfLogN (by intro _ _; rfl) requiredLen + end KoalaBear end NTT end CPolynomial diff --git a/CompPoly/Univariate/Raw/Context.lean b/CompPoly/Univariate/Raw/Context.lean new file mode 100644 index 00000000..6f625739 --- /dev/null +++ b/CompPoly/Univariate/Raw/Context.lean @@ -0,0 +1,133 @@ +/- +Copyright (c) 2026 CompPoly. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ +import CompPoly.Univariate.DivisionCorrectness +import CompPoly.Univariate.NTT.FastMul +import CompPoly.Univariate.NTT.FastMulLow +import CompPoly.Univariate.NTTFast.Correctness.Pipeline +import CompPoly.Univariate.NTTFast.FastMulLow + +/-! +# Raw Univariate Algorithm Contexts + +Array-level execution dictionaries for reusable raw univariate polynomial +kernels. +-/ + +namespace CompPoly + +namespace CPolynomial + +namespace Raw + +variable {R : Type*} + +/-- Internal raw multiplication backend for array-level polynomial kernels. -/ +structure MulContext (R : Type*) [Semiring R] [BEq R] [LawfulBEq R] where + /-- Multiply two raw polynomials, returning the trimmed raw product. -/ + mul : CPolynomial.Raw R → CPolynomial.Raw R → CPolynomial.Raw R + /-- The backend agrees with raw polynomial multiplication. -/ + mul_eq_mul : ∀ p q, mul p q = p * q + +/-- Internal raw monic-remainder backend for array-level polynomial kernels. -/ +structure ModContext (R : Type*) [Field R] [BEq R] [LawfulBEq R] where + /-- Reduce the first raw polynomial modulo the second raw monic divisor. -/ + modByMonic : CPolynomial.Raw R → CPolynomial.Raw R → CPolynomial.Raw R + /-- The backend agrees with raw monic remainders for canonical raw inputs. -/ + modByMonic_eq_modByMonic : + ∀ p q, p.trim = p → q.trim = q → modByMonic p q = CPolynomial.Raw.modByMonic p q + +namespace MulContext + +/-- The default raw multiplication context, backed by raw polynomial multiplication. -/ +def naive [Semiring R] [BEq R] [LawfulBEq R] : MulContext R where + mul p q := p * q + mul_eq_mul _ _ := rfl + +/-- NTT-backed raw multiplication context with raw multiplication for unsupported lengths. -/ +def ntt [Field R] [BEq R] [LawfulBEq R] + (bestDomainForLength? : (requiredLen : Nat) → + Option (NTT.FittingDomain R requiredLen)) : + MulContext R where + mul p q := + let requiredLen := NTT.Domain.requiredLength p q + match bestDomainForLength? requiredLen with + | some ⟨D, _⟩ => (NTT.FastMul.Raw.fastMulImpl D p q).trim + | none => p * q + mul_eq_mul p q := by + let requiredLen := NTT.Domain.requiredLength p q + cases hdomain : bestDomainForLength? requiredLen with + | none => + simp [hdomain, requiredLen] + | some fitted => + rcases fitted with ⟨D, hfit⟩ + simp [hdomain, requiredLen, NTT.FastMul.Raw.fastMulImpl_trim_eq_mul D p q (by + simpa [NTT.Domain.fits] using hfit)] + +/-- NTTFast-backed raw multiplication context with raw multiplication for unsupported lengths. -/ +def nttFast [Field R] [BEq R] [LawfulBEq R] + (bestDomainForLength? : (requiredLen : Nat) → + Option (NTT.FittingDomain R requiredLen)) : + MulContext R where + mul p q := + let requiredLen := NTT.Domain.requiredLength p q + match bestDomainForLength? requiredLen with + | some ⟨D, _⟩ => (NTTFast.Raw.fastMulImpl D p q).trim + | none => p * q + mul_eq_mul p q := by + let requiredLen := NTT.Domain.requiredLength p q + cases hdomain : bestDomainForLength? requiredLen with + | none => + simp [hdomain, requiredLen] + | some fitted => + rcases fitted with ⟨D, hfit⟩ + simp [hdomain, requiredLen, NTTFast.Raw.fastMulImpl_trim_eq_mul D p q (by + simpa [NTT.Domain.fits] using hfit)] + +end MulContext + +namespace ModContext + +/-- The default raw monic-remainder context, backed by raw `modByMonic`. -/ +def naive [Field R] [BEq R] [LawfulBEq R] : ModContext R where + modByMonic p q := CPolynomial.Raw.modByMonic p q + modByMonic_eq_modByMonic _ _ _ _ := rfl + +/-- Raw remainder-only monic-remainder context. -/ +def remainderOnly [Field R] [BEq R] [LawfulBEq R] : ModContext R where + modByMonic p q := CPolynomial.Raw.modByMonicRemainderOnly p q + modByMonic_eq_modByMonic p q _ _ := CPolynomial.Raw.modByMonicRemainderOnly_eq_modByMonic p q + +/-- Raw reversal-based monic-remainder context parameterized by low-product multiplication. -/ +def reversal [Field R] [BEq R] [LawfulBEq R] + (M : Raw.MulLowContext R) : ModContext R where + modByMonic p q := CPolynomial.Raw.modByMonicByReversal M p q + modByMonic_eq_modByMonic p q hp hq := by + have h := congrArg Subtype.val + (CPolynomial.modByMonicByReversal_eq_modByMonic M + (CPolynomial.ofArray p) (CPolynomial.ofArray q)) + simpa [CPolynomial.ofArray, hp, hq] using h + +/-- Raw monic remainders by reversal, using an NTT low-product backend. -/ +def reversalNtt [Field R] [BEq R] [LawfulBEq R] + (bestDomainForLength? : (requiredLen : Nat) → + Option (NTT.FittingDomain R requiredLen)) : + ModContext R := + reversal (NTT.FastMulLow.withFallback bestDomainForLength?) + +/-- Raw monic remainders by reversal, using an NTTFast low-product backend. -/ +def reversalNttFast [Field R] [BEq R] [LawfulBEq R] + (bestDomainForLength? : (requiredLen : Nat) → + Option (NTT.FittingDomain R requiredLen)) : + ModContext R := + reversal (NTTFast.FastMulLow.withFallback bestDomainForLength?) + +end ModContext + +end Raw + +end CPolynomial + +end CompPoly diff --git a/CompPoly/Univariate/Raw/Division.lean b/CompPoly/Univariate/Raw/Division.lean index d8fcc46f..53bbb5e1 100644 --- a/CompPoly/Univariate/Raw/Division.lean +++ b/CompPoly/Univariate/Raw/Division.lean @@ -105,8 +105,8 @@ Remainder by a monic polynomial through reversal and truncated products. For canonical monic inputs this computes the quotient from the reversed divisor inverse modulo `X^k`, then subtracts only the low coefficients needed for the -remainder. Inputs outside that executable contract fall back to the -remainder-only implementation. +remainder. Inputs outside the fast-path guard use the simple monic-remainder +implementation. -/ @[inline, specialize] def modByMonicByReversal [Field R] [LawfulBEq R] (M : MulLowContext R) @@ -138,6 +138,30 @@ def mod [Field R] (p q : CPolynomial.Raw R) : CPolynomial.Raw R := instance [Field R] : Div (CPolynomial.Raw R) := ⟨div⟩ instance [Field R] : Mod (CPolynomial.Raw R) := ⟨mod⟩ +/-- Normalize a nonzero raw polynomial to monic form. The zero polynomial stays zero. -/ +def monicNormalize [Field R] (p : CPolynomial.Raw R) : CPolynomial.Raw R := + let p := p.trim + if p == 0 then + 0 + else + p.leadingCoeff⁻¹ • p + +/-- Raw Euclidean gcd with explicit fuel, normalized to a monic result. -/ +def gcdMonicWithFuel [Field R] : + Nat → CPolynomial.Raw R → CPolynomial.Raw R → CPolynomial.Raw R + | 0, p, _ => monicNormalize p + | fuel + 1, p, q => + let p := p.trim + let q := q.trim + if q == 0 then + monicNormalize p + else + gcdMonicWithFuel fuel q (p % q) + +/-- Raw monic Euclidean gcd. -/ +def gcdMonic [Field R] (p q : CPolynomial.Raw R) : CPolynomial.Raw R := + gcdMonicWithFuel (p.size + q.size + 1) p q + end Division end CPolynomial.Raw diff --git a/CompPoly/Univariate/Raw/Modular.lean b/CompPoly/Univariate/Raw/Modular.lean new file mode 100644 index 00000000..2c93b17d --- /dev/null +++ b/CompPoly/Univariate/Raw/Modular.lean @@ -0,0 +1,76 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Raw.Context + +/-! +# Raw Modular Operations on Univariate Polynomials + +Context-parametric modular multiplication and exponentiation over raw +univariate polynomials. +-/ + +namespace CompPoly + +namespace CPolynomial.Raw + +/-- Raw multiplication modulo a polynomial, treating zero modulus as no reduction. -/ +def mulModWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (modulus p q : CPolynomial.Raw F) : CPolynomial.Raw F := + let product := M.mul p q + if modulus.trim == 0 then + product + else + D.modByMonic product (monicNormalize modulus) + +/-- Raw binary modular exponentiation accumulator. -/ +def powModBinaryAuxWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (modulus : CPolynomial.Raw F) : + Nat → CPolynomial.Raw F → CPolynomial.Raw F → CPolynomial.Raw F + | 0, acc, _ => acc + | n + 1, acc, current => + let acc' := + if (n + 1) % 2 == 1 then + mulModWith M D modulus acc current + else + acc + let current' := mulModWith M D modulus current current + powModBinaryAuxWith M D modulus ((n + 1) / 2) acc' current' +termination_by exp _ _ => exp +decreasing_by + exact Nat.div_lt_self (Nat.succ_pos n) (by decide) + +/-- Raw modular exponentiation by repeated squaring. -/ +def powModWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (modulus base : CPolynomial.Raw F) (exponent : Nat) : CPolynomial.Raw F := + let oneMod := + if modulus.trim == 0 then + (1 : CPolynomial.Raw F) + else + D.modByMonic (1 : CPolynomial.Raw F) (monicNormalize modulus) + powModBinaryAuxWith M D modulus exponent oneMod base + +/-- Raw `X mod modulus`, with zero modulus treated as no reduction. -/ +def xModWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (D : CPolynomial.Raw.ModContext F) : CPolynomial.Raw F → CPolynomial.Raw F + | modulus => + if modulus.trim == 0 then + (CPolynomial.Raw.X : CPolynomial.Raw F) + else + D.modByMonic (CPolynomial.Raw.X : CPolynomial.Raw F) (monicNormalize modulus) + +/-- Raw `(X^q mod modulus) - (X mod modulus)`, without materializing `X^q - X`. -/ +def xPowSubXModWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (q : Nat) (modulus : CPolynomial.Raw F) : CPolynomial.Raw F := + powModWith M D modulus (CPolynomial.Raw.X : CPolynomial.Raw F) q - xModWith D modulus + +end CPolynomial.Raw + +end CompPoly diff --git a/CompPoly/Univariate/Roots.lean b/CompPoly/Univariate/Roots.lean new file mode 100644 index 00000000..761c2d0c --- /dev/null +++ b/CompPoly/Univariate/Roots.lean @@ -0,0 +1,15 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Roots.Correctness +import CompPoly.Univariate.Roots.Enumeration +import CompPoly.Univariate.Roots.SmoothSubgroup + +/-! +# Univariate Root Finding + +Public univariate root-finding surface. +-/ diff --git a/CompPoly/Univariate/Roots/Backend.lean b/CompPoly/Univariate/Roots/Backend.lean new file mode 100644 index 00000000..3aa98d1f --- /dev/null +++ b/CompPoly/Univariate/Roots/Backend.lean @@ -0,0 +1,56 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Roots.Extraction +import CompPoly.Univariate.Roots.RootProduct + +/-! +# Finite-Field Root Backend + +Executable field-root extraction over finite fields. The public operation +handles zero, constant, and linear cases explicitly, computes the finite-field +root product modulo the input polynomial, splits the product into linear +factors, then validates and deduplicates candidates against the original input. +-/ + +namespace CompPoly + +namespace CPolynomial + +namespace Roots + +namespace FiniteField + +/-- Executable roots of a univariate polynomial over a finite field. -/ +def rootsInFiniteFieldWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : FiniteFieldContext F) (splitter : LinearFactorProductSplitter F) + (p : CPolynomial F) : Array F := + if p == 0 then + #[] + else if p.val.size ≤ 1 then + #[] + else if p.val.size = 2 && !(p.coeff 1 == 0) then + CPolynomial.rootsFromLinearFactors p #[p] + else + let rootProduct := finiteFieldRootProductWith M D ctx p + let factors := splitter.splitLinearFactors ctx.q rootProduct + CPolynomial.rootsFromLinearFactors p factors + +/-- Executable roots using the default raw multiplication and monic-remainder backends. -/ +def rootsInFiniteField {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (ctx : FiniteFieldContext F) (splitter : LinearFactorProductSplitter F) + (p : CPolynomial F) : Array F := + rootsInFiniteFieldWith CPolynomial.Raw.MulContext.naive + CPolynomial.Raw.ModContext.naive ctx splitter p + +end FiniteField + +end Roots + +end CPolynomial + +end CompPoly diff --git a/CompPoly/Univariate/Roots/Context.lean b/CompPoly/Univariate/Roots/Context.lean new file mode 100644 index 00000000..5716c51d --- /dev/null +++ b/CompPoly/Univariate/Roots/Context.lean @@ -0,0 +1,77 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Basic +import Mathlib.FieldTheory.Finite.Basic + +/-! +# Finite-Field Root Contexts + +Explicit contexts for executable univariate root finding over finite fields. +The algorithms use the cardinality carried here; downstream proofs use the +finite-field and splitter contracts without unfolding concrete implementations. +-/ + +namespace CompPoly + +namespace CPolynomial + +namespace Roots + +namespace FiniteField + +/-- Finite-field facts needed by executable root extraction. -/ +structure FiniteFieldContext (F : Type*) [Field F] where + q : Nat + finite : Finite F + card_eq : Nat.card F = q + frobenius_fixed : ∀ a : F, a ^ q = a + +/-- A polynomial represented as a nonconstant linear factor. -/ +def IsLinearFactor {F : Type*} [Field F] (factor : CPolynomial F) : Prop := + factor.val.size ≤ 2 ∧ factor.coeff 1 ≠ 0 + +/-- A linear factor whose extracted root is `a`. -/ +def IsLinearRootFactorCandidate {F : Type*} [Field F] + (factor : CPolynomial F) (a : F) : Prop := + IsLinearFactor factor ∧ factor.coeff 0 + factor.coeff 1 * a = 0 + +/-- A deterministic splitter for squarefree products of linear factors. + +The executable function consumes the field cardinality and the current factor. +Completeness of the public root backend depends on splitter completeness for +every linear-factor product reached by the recursion. +-/ +structure LinearFactorProductSplitter (F : Type*) [Field F] [BEq F] [LawfulBEq F] where + splitLinearFactors : Nat → CPolynomial F → Array (CPolynomial F) + /-- + The precondition under which `complete` is claimed for the splitter input. + + Root backends should establish this predicate for the field-root product they + pass to the splitter. Executable splitters may remain defensive outside this + predicate, but completeness is only part of the contract under it. + -/ + validInput : Nat → CPolynomial F → Prop := fun _ _ ↦ True + sound : + ∀ q p factor, + factor ∈ (splitLinearFactors q p).toList → + IsLinearFactor factor + complete : + ∀ q p a, + validInput q p → + p ≠ 0 → + CPolynomial.eval a p = 0 → + ∃ factor, + factor ∈ (splitLinearFactors q p).toList ∧ + IsLinearRootFactorCandidate factor a + +end FiniteField + +end Roots + +end CPolynomial + +end CompPoly diff --git a/CompPoly/Univariate/Roots/Correctness.lean b/CompPoly/Univariate/Roots/Correctness.lean new file mode 100644 index 00000000..c6300d77 --- /dev/null +++ b/CompPoly/Univariate/Roots/Correctness.lean @@ -0,0 +1,1092 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Roots.Backend +import CompPoly.Univariate.Roots.Splitter +import CompPoly.Univariate.EuclideanAlgorithm +import CompPoly.Univariate.ToPoly.Impl + +/-! +# Finite-Field Root Correctness + +Theorem statements and certified context constructors for the executable +finite-field root backend. +-/ + +namespace CompPoly + +namespace CPolynomial + +private lemma mem_eraseDups_fold {α : Type*} [BEq α] [LawfulBEq α] + (a : α) : ∀ (xs : List α) (out : Array α), + a ∈ List.foldl (fun out x ↦ if x ∈ out then out else out.push x) out xs → + a ∈ out ∨ a ∈ xs := by + intro xs + induction xs with + | nil => + intro out h + exact Or.inl h + | cons x xs ih => + intro out h + simp only [List.foldl_cons] at h + by_cases hx : x ∈ out + · have h' := ih out (by simpa [hx] using h) + cases h' with + | inl hout => exact Or.inl hout + | inr hxs => exact Or.inr (by simp [hxs]) + · have h' := ih (out.push x) (by simpa [hx] using h) + cases h' with + | inl hout => + simp at hout + cases hout with + | inl hout => exact Or.inl hout + | inr hax => exact Or.inr (by simp [hax]) + | inr hxs => exact Or.inr (by simp [hxs]) + +private lemma mem_eraseDups {α : Type*} [BEq α] [LawfulBEq α] + {xs : Array α} {a : α} (h : a ∈ xs.eraseDups) : a ∈ xs := by + unfold Array.eraseDups at h + rcases xs with ⟨l⟩ + simp at h ⊢ + have hh := mem_eraseDups_fold a l #[] h + simpa using hh + +private lemma mem_eraseDups_fold_of_mem {α : Type*} [BEq α] [LawfulBEq α] + (a : α) : ∀ (xs : List α) (out : Array α), + a ∈ out ∨ a ∈ xs → + a ∈ List.foldl (fun out x ↦ if x ∈ out then out else out.push x) out xs := by + intro xs + induction xs with + | nil => + intro out h + simpa using h + | cons x xs ih => + intro out h + simp only [List.foldl_cons] + by_cases hx : x ∈ out + · simp [hx] + apply ih out + cases h with + | inl hout => exact Or.inl hout + | inr hmem => + simp at hmem + cases hmem with + | inl hax => exact Or.inl (by simpa [hax] using hx) + | inr hxs => exact Or.inr hxs + · simp [hx] + apply ih (out.push x) + cases h with + | inl hout => exact Or.inl (by simp [hout]) + | inr hmem => + simp at hmem + cases hmem with + | inl hax => exact Or.inl (by simp [hax]) + | inr hxs => exact Or.inr hxs + +private lemma mem_eraseDups_of_mem {α : Type*} [BEq α] [LawfulBEq α] + {xs : Array α} {a : α} (h : a ∈ xs) : a ∈ xs.eraseDups := by + unfold Array.eraseDups + rcases xs with ⟨l⟩ + simp at h ⊢ + exact mem_eraseDups_fold_of_mem a l #[] (Or.inr h) + +private theorem Raw.toPoly_smul {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (c : F) (p : CPolynomial.Raw F) : + (CPolynomial.Raw.smul c p).toPoly = c • p.toPoly := by + ext i + rw [Polynomial.coeff_smul] + rw [CPolynomial.Raw.coeff_toPoly, CPolynomial.Raw.coeff_toPoly] + exact CPolynomial.Raw.smul_coeff c p i + +private theorem Raw.eval_smul {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (c a : F) (p : CPolynomial.Raw F) : + (CPolynomial.Raw.smul c p).eval a = c * p.eval a := by + rw [← CPolynomial.Raw.eval_toPoly_eq_eval, Raw.toPoly_smul] + rw [Polynomial.smul_eq_C_mul, Polynomial.eval_mul, Polynomial.eval_C, + CPolynomial.Raw.eval_toPoly_eq_eval] + +private theorem Raw.eval_mul {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (a : F) (p q : CPolynomial.Raw F) : + (p * q).eval a = p.eval a * q.eval a := by + rw [← CPolynomial.Raw.eval_toPoly_eq_eval, CPolynomial.Raw.toPoly_mul, + Polynomial.eval_mul, CPolynomial.Raw.eval_toPoly_eq_eval, + CPolynomial.Raw.eval_toPoly_eq_eval] + +private theorem Raw.eval_C {F : Type*} [Field F] (a c : F) : + (CPolynomial.Raw.C c).eval a = c := by + rw [← CPolynomial.Raw.eval_toPoly_eq_eval] + simp [CPolynomial.Raw.toPoly_C] + +private theorem Raw.eval_C_smul {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (c a : F) (p : CPolynomial.Raw F) : + (CPolynomial.Raw.C c • p).eval a = c * p.eval a := by + rw [smul_eq_mul, Raw.eval_mul, Raw.eval_C] + +private theorem Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p : CPolynomial.Raw F} {a : F} (hp : p.eval a = 0) : + (CPolynomial.Raw.monicNormalize p).eval a = 0 := by + unfold CPolynomial.Raw.monicNormalize + by_cases hzero : p.trim = (#[] : CPolynomial.Raw F) + · simp [hzero, CPolynomial.Raw.eval, CPolynomial.Raw.eval₂] + · simp [hzero] + change CPolynomial.Raw.eval a (CPolynomial.Raw.smul p.trim.leadingCoeff⁻¹ p.trim) = 0 + rw [Raw.eval_smul, CPolynomial.Raw.eval_trim_eq_eval, hp] + simp + +private theorem Raw.eval_mod_eq_zero_of_left_right {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p q : CPolynomial.Raw F} {a : F} + (hp : p.eval a = 0) (hq : q.eval a = 0) : + (p % q).eval a = 0 := by + change (CPolynomial.Raw.mod p q).eval a = 0 + unfold CPolynomial.Raw.mod + have hqscaled : (CPolynomial.Raw.C (q.leadingCoeff)⁻¹ * q).eval a = 0 := by + rw [Raw.eval_mul, Raw.eval_C, hq] + simp + have hmod := CPolynomial.Raw.eval_modByMonic_eq_self_of_eval_eq_zero + (CPolynomial.Raw.C (q.leadingCoeff)⁻¹ • p) + (CPolynomial.Raw.C (q.leadingCoeff)⁻¹ * q) hqscaled + rw [hmod, Raw.eval_C_smul, hp] + simp + +private theorem Raw.eval_gcdMonicWithFuel_eq_zero_of_left_right {F : Type*} + [Field F] [BEq F] [LawfulBEq F] {a : F} : + ∀ fuel (p q : CPolynomial.Raw F), + p.eval a = 0 → q.eval a = 0 → + (CPolynomial.Raw.gcdMonicWithFuel fuel p q).eval a = 0 := by + intro fuel + induction fuel with + | zero => + intro p q hp _hq + unfold CPolynomial.Raw.gcdMonicWithFuel + exact Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero hp + | succ fuel ih => + intro p q hp hq + unfold CPolynomial.Raw.gcdMonicWithFuel + by_cases hqzero : q.trim = (#[] : CPolynomial.Raw F) + · simp [hqzero] + exact Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero (by + rw [CPolynomial.Raw.eval_trim_eq_eval] + exact hp) + · simp [hqzero] + apply ih + · rw [CPolynomial.Raw.eval_trim_eq_eval] + exact hq + · exact Raw.eval_mod_eq_zero_of_left_right + (by + rw [CPolynomial.Raw.eval_trim_eq_eval] + exact hp) + (by + rw [CPolynomial.Raw.eval_trim_eq_eval] + exact hq) + +/-- Monic normalization preserves roots of nonzero polynomials. -/ +theorem monicNormalize_root_iff {F : Type*} [Field F] [BEq F] [LawfulBEq F] + {p : CPolynomial F} {a : F} (hp : p ≠ 0) : + CPolynomial.eval a (monicNormalize p) = 0 ↔ CPolynomial.eval a p = 0 := by + unfold monicNormalize CPolynomial.Raw.monicNormalize + have hpempty : ¬ (p.val : CPolynomial.Raw F) = #[] := by + intro h + apply hp + apply CPolynomial.ext + simpa using h + have hsize : 0 < p.val.size := by + cases hs : p.val.size with + | zero => + have hval : p.val = (#[] : CPolynomial.Raw F) := + Array.eq_empty_of_size_eq_zero hs + exact (hpempty hval).elim + | succ _ => omega + have hgetLastD : + (p.val : CPolynomial.Raw F).getLastD 0 = + (p.val : CPolynomial.Raw F).getLast hsize := by + unfold Array.getLastD Array.getLast + simp [hsize] + have hlc : (p.val : CPolynomial.Raw F).getLastD 0 ≠ 0 := by + rw [hgetLastD] + exact p.property hsize + simp [CPolynomial.Raw.leadingCoeff, hpempty] + change + CPolynomial.Raw.eval a + (((Array.getLastD (p.val : CPolynomial.Raw F) 0)⁻¹ • p.val).trim) = + 0 ↔ + CPolynomial.Raw.eval a p.val = 0 + rw [CPolynomial.Raw.eval_trim_eq_eval] + change + CPolynomial.Raw.eval a + (CPolynomial.Raw.smul ((Array.getLastD (p.val : CPolynomial.Raw F) 0)⁻¹) + p.val) = + 0 ↔ + CPolynomial.Raw.eval a p.val = 0 + rw [CPolynomial.Raw.eval_smul] + constructor + · intro h + exact (mul_eq_zero.mp h).resolve_left (inv_ne_zero hlc) + · intro h + simp [h] + +/-- Monic normalization preserves roots, including the zero polynomial case. -/ +theorem monicNormalize_root_of_root {F : Type*} [Field F] [BEq F] [LawfulBEq F] + {p : CPolynomial F} {a : F} (hp : CPolynomial.eval a p = 0) : + CPolynomial.eval a (monicNormalize p) = 0 := by + unfold monicNormalize + change CPolynomial.Raw.eval a (CPolynomial.Raw.monicNormalize p.val).trim = 0 + rw [CPolynomial.Raw.eval_trim_eq_eval] + exact Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero hp + +/-- The normalized polynomial divides the original polynomial. -/ +theorem toPoly_monicNormalize_dvd_self {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (p : CPolynomial F) : + (CPolynomial.monicNormalize p).toPoly ∣ p.toPoly := by + letI : DecidableEq F := instDecidableEqOfLawfulBEq + rw [CPolynomial.monicNormalize_toPoly_eq_normalize] + exact (normalize_associated p.toPoly).dvd + +/-- The monic gcd contains every common root. -/ +theorem gcdMonic_root_of_left_right {F : Type*} [Field F] [BEq F] [LawfulBEq F] + {p q : CPolynomial F} {a : F} + (hp : CPolynomial.eval a p = 0) (hq : CPolynomial.eval a q = 0) : + CPolynomial.eval a (gcdMonic p q) = 0 := by + unfold gcdMonic + change CPolynomial.Raw.eval a (CPolynomial.Raw.gcdMonic p.val q.val).trim = 0 + rw [CPolynomial.Raw.eval_trim_eq_eval] + unfold CPolynomial.Raw.gcdMonic + exact Raw.eval_gcdMonicWithFuel_eq_zero_of_left_right _ p.val q.val hp hq + +/-- The monic gcd divides its left operand. -/ +theorem toPoly_gcdMonic_dvd_left {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (p q : CPolynomial F) : + (CPolynomial.gcdMonic p q).toPoly ∣ p.toPoly := by + letI : DecidableEq F := instDecidableEqOfLawfulBEq + rw [CPolynomial.gcdMonic_toPoly_eq_normalize_gcd] + exact (normalize_associated (EuclideanDomain.gcd p.toPoly q.toPoly)).dvd.trans + (EuclideanDomain.gcd_dvd_left p.toPoly q.toPoly) + +/-- The monic gcd divides its right operand. -/ +theorem toPoly_gcdMonic_dvd_right {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (p q : CPolynomial F) : + (CPolynomial.gcdMonic p q).toPoly ∣ q.toPoly := by + letI : DecidableEq F := instDecidableEqOfLawfulBEq + rw [CPolynomial.gcdMonic_toPoly_eq_normalize_gcd] + exact (normalize_associated (EuclideanDomain.gcd p.toPoly q.toPoly)).dvd.trans + (EuclideanDomain.gcd_dvd_right p.toPoly q.toPoly) + +/-- The monic gcd vanishes exactly at common roots. -/ +theorem gcdMonic_root_iff_left_right {F : Type*} + [Field F] [BEq F] [LawfulBEq F] [DecidableEq F] + {p q : CPolynomial F} {a : F} : + CPolynomial.eval a (gcdMonic p q) = 0 ↔ + CPolynomial.eval a p = 0 ∧ CPolynomial.eval a q = 0 := by + rw [CPolynomial.eval_toPoly, CPolynomial.gcdMonic_toPoly_eq_normalize_gcd] + rw [← Polynomial.IsRoot.def] + rw [normalize_apply] + by_cases h : EuclideanDomain.gcd p.toPoly q.toPoly = 0 + · rw [h] + rw [Polynomial.IsRoot.def] + rw [Polynomial.eval_mul, Polynomial.eval_zero] + ring_nf + have hleft : p.toPoly = 0 := (EuclideanDomain.gcd_eq_zero_iff.mp h).1 + have hright : q.toPoly = 0 := (EuclideanDomain.gcd_eq_zero_iff.mp h).2 + simp only [CPolynomial.eval_toPoly, hleft, hright, Polynomial.eval_zero, and_self] + · rw [Polynomial.coe_normUnit_of_ne_zero h] + rw [Polynomial.IsRoot.def] + rw [Polynomial.eval_mul, Polynomial.eval_C, mul_eq_zero] + have hlead : (EuclideanDomain.gcd p.toPoly q.toPoly).leadingCoeff ≠ 0 := + Polynomial.leadingCoeff_ne_zero.mpr h + have hinv : (EuclideanDomain.gcd p.toPoly q.toPoly).leadingCoeff⁻¹ ≠ 0 := + inv_ne_zero hlead + rw [or_iff_left hinv] + rw [← Polynomial.IsRoot.def] + have hiff := Polynomial.isRoot_gcd_iff_isRoot_left_right + (f := p.toPoly) (g := q.toPoly) (α := a) + simpa [Polynomial.IsRoot, CPolynomial.eval_toPoly] using hiff + +/-- A root of a dividend is a root of the exact quotient when the divisor does +not vanish at that point. -/ +theorem eval_div_eq_zero_of_dvd_of_root_of_ne_root {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p q : CPolynomial F} {a : F} + (hdiv : q.toPoly ∣ p.toPoly) + (hp : CPolynomial.eval a p = 0) + (hq : CPolynomial.eval a q ≠ 0) : + CPolynomial.eval a (p / q) = 0 := by + rcases hdiv with ⟨r, hr⟩ + have hqPoly : q.toPoly ≠ 0 := by + intro hqPoly + apply hq + rw [CPolynomial.eval_toPoly, hqPoly, Polynomial.eval_zero] + rw [CPolynomial.eval_toPoly] + change Polynomial.eval a (CPolynomial.div p q).toPoly = 0 + rw [CPolynomial.div_toPoly_eq_div] + have hdivPoly : p.toPoly / q.toPoly = r := by + exact (EuclideanDomain.eq_div_of_mul_eq_right hqPoly hr.symm).symm + rw [hdivPoly] + have hpEval : Polynomial.eval a p.toPoly = 0 := by + simpa [CPolynomial.eval_toPoly] using hp + rw [hr, Polynomial.eval_mul] at hpEval + exact (mul_eq_zero.mp hpEval).resolve_left (by + simpa [CPolynomial.eval_toPoly] using hq) + +/-- Monic normalization of an exact quotient preserves the quotient root from +`eval_div_eq_zero_of_dvd_of_root_of_ne_root`. -/ +theorem monicNormalize_div_root_of_dvd_of_root_of_ne_root {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p q : CPolynomial F} {a : F} + (hdiv : q.toPoly ∣ p.toPoly) + (hp : CPolynomial.eval a p = 0) + (hq : CPolynomial.eval a q ≠ 0) : + CPolynomial.eval a (CPolynomial.monicNormalize (p / q)) = 0 := + monicNormalize_root_of_root + (eval_div_eq_zero_of_dvd_of_root_of_ne_root hdiv hp hq) + +/-- The normalized extended gcd contains every common root. -/ +theorem normXgcd_root_of_left_right {F : Type*} [Field F] [BEq F] [LawfulBEq F] + {p q : CPolynomial F} {a : F} + (hp : CPolynomial.eval a p = 0) (hq : CPolynomial.eval a q = 0) : + CPolynomial.eval a (normXgcd p q).1 = 0 := by + have hbez := normXgcd_bezout p q 0 + simp only [Bezout] at hbez + rw [hbez] + rw [CPolynomial.eval_toPoly, CPolynomial.toPoly_add, CPolynomial.toPoly_mul, + CPolynomial.toPoly_mul, Polynomial.eval_add, Polynomial.eval_mul, Polynomial.eval_mul, + ← CPolynomial.eval_toPoly, ← CPolynomial.eval_toPoly, + ← CPolynomial.eval_toPoly, ← CPolynomial.eval_toPoly] + simp [hp, hq] + +/-- The specialized monic gcd and normalized extended gcd have the same root +predicate. -/ +theorem eval_gcdMonic_eq_zero_iff_normXgcd_fst_eq_zero + {F : Type*} [Field F] [BEq F] [LawfulBEq F] [DecidableEq F] + (a : F) (p q : CPolynomial F) : + CPolynomial.eval a (CPolynomial.gcdMonic p q) = 0 ↔ + CPolynomial.eval a ((CPolynomial.normXgcd p q).1) = 0 := by + rw [CPolynomial.gcdMonic_eq_normXgcd_fst] + +/-- The normalized extended-gcd first component contains every common root. -/ +theorem normXgcd_fst_root_of_left_right + {F : Type*} [Field F] [BEq F] [LawfulBEq F] [DecidableEq F] + {a : F} {p q : CPolynomial F} + (hp : CPolynomial.eval a p = 0) + (hq : CPolynomial.eval a q = 0) : + CPolynomial.eval a ((CPolynomial.normXgcd p q).1) = 0 := by + exact normXgcd_root_of_left_right hp hq + +/-- Roots extracted from a linear factor satisfy that factor. -/ +theorem linearRootOfFactor?_sound {F : Type*} [Field F] [BEq F] [LawfulBEq F] + {factor : CPolynomial F} {a : F} + (h : linearRootOfFactor? factor = some a) : + CPolynomial.eval a factor = 0 := by + rcases factor with ⟨⟨xs⟩, hcanon⟩ + cases xs with + | nil => simp [linearRootOfFactor?] at h + | cons x xs => + cases xs with + | nil => simp [linearRootOfFactor?] at h + | cons y xs => + cases xs with + | nil => + simp [linearRootOfFactor?, CPolynomial.eval] at h ⊢ + rcases h with ⟨hy, ha⟩ + rw [← ha] + field_simp [hy] + ring + | cons z xs => simp [linearRootOfFactor?] at h + +/-- Validation makes returned candidates sound for the original polynomial. -/ +theorem mem_validateRootCandidates_eval_eq_zero {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p : CPolynomial F} {candidates : Array F} {a : F} + (h : a ∈ (validateRootCandidates p candidates).toList) : + CPolynomial.eval a p = 0 := by + rw [validateRootCandidates] at h + simp at h + simpa [CPolynomial.eval_horner_eq_eval] using h.2 + +private theorem mem_rootsFromLinearFactors_eval_eq_zero {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p : CPolynomial F} {factors : Array (CPolynomial F)} {a : F} + (h : a ∈ (rootsFromLinearFactors p factors).toList) : + CPolynomial.eval a p = 0 := by + rw [rootsFromLinearFactors] at h + have h' : + a ∈ (validateRootCandidates p + (List.filterMap linearRootOfFactor? factors.toList).toArray).toList := by + have hm : + a ∈ validateRootCandidates p + (List.filterMap linearRootOfFactor? factors.toList).toArray := + mem_eraseDups (by simpa using h) + simpa using hm + exact mem_validateRootCandidates_eval_eq_zero h' + +namespace Roots + +namespace FiniteField + +private theorem linearRootOfFactor?_eq_some_of_candidate {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {factor : CPolynomial F} {a : F} + (h : IsLinearRootFactorCandidate factor a) : + CPolynomial.linearRootOfFactor? factor = some a := by + rw [CPolynomial.linearRootOfFactor?] + have hcond : factor.val.size ≤ 2 ∧ factor.coeff 1 ≠ 0 := h.1 + rw [if_pos] + · congr + change -(factor.coeff 0) / factor.coeff 1 = a + apply (div_eq_iff hcond.2).2 + rw [neg_eq_iff_add_eq_zero] + rw [_root_.mul_comm a (factor.coeff 1)] + exact h.2 + · simp [hcond] + intro hbad + exact hcond.2 (by simpa [CPolynomial.coeff, CPolynomial.Raw.coeff] using hbad) + +private theorem mem_rootsFromLinearFactors_of_candidate {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p : CPolynomial F} {factors : Array (CPolynomial F)} {factor : CPolynomial F} {a : F} + (hmem : factor ∈ factors.toList) (hcand : IsLinearRootFactorCandidate factor a) + (hroot : CPolynomial.eval a p = 0) : + a ∈ (CPolynomial.rootsFromLinearFactors p factors).toList := by + rw [CPolynomial.rootsFromLinearFactors] + have hvalid : a ∈ CPolynomial.validateRootCandidates p + (List.filterMap CPolynomial.linearRootOfFactor? factors.toList).toArray := by + rw [CPolynomial.validateRootCandidates] + simp [hroot, CPolynomial.eval_horner_eq_eval] + exact ⟨factor, by simpa using hmem, linearRootOfFactor?_eq_some_of_candidate hcand⟩ + have herase := mem_eraseDups_of_mem hvalid + simpa using herase + +private theorem raw_eval_mulMod_naive_eq_mul {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {modulus p q : CPolynomial.Raw F} {a : F} + (hmod : modulus.eval a = 0) : + (CPolynomial.Raw.mulModWith CPolynomial.Raw.MulContext.naive + CPolynomial.Raw.ModContext.naive modulus p q).eval a = + p.eval a * q.eval a := by + unfold CPolynomial.Raw.mulModWith CPolynomial.Raw.MulContext.naive + CPolynomial.Raw.ModContext.naive + by_cases hzero : modulus.trim = (#[] : CPolynomial.Raw F) + · simp [hzero, CPolynomial.Raw.eval_mul] + · simp [hzero] + have hroot : (CPolynomial.Raw.monicNormalize modulus).eval a = 0 := + CPolynomial.Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero hmod + have hmodBy := CPolynomial.Raw.eval_modByMonic_eq_self_of_eval_eq_zero + (p * q) (CPolynomial.Raw.monicNormalize modulus) hroot + rw [hmodBy, CPolynomial.Raw.eval_mul] + +private theorem raw_eval_powModBinaryAux_naive {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {modulus acc current : CPolynomial.Raw F} {a : F} + (hmod : modulus.eval a = 0) : + ∀ n, + (CPolynomial.Raw.powModBinaryAuxWith CPolynomial.Raw.MulContext.naive + CPolynomial.Raw.ModContext.naive modulus n acc current).eval a = + acc.eval a * current.eval a ^ n := by + intro n + induction n using Nat.strongRecOn generalizing acc current with + | ind n ih => + cases n with + | zero => + simp [CPolynomial.Raw.powModBinaryAuxWith] + | succ n => + rw [CPolynomial.Raw.powModBinaryAuxWith] + have ih' := ih ((n + 1) / 2) + (Nat.div_lt_self (Nat.succ_pos n) (by decide)) + (acc := + if (n + 1) % 2 == 1 then + CPolynomial.Raw.mulModWith CPolynomial.Raw.MulContext.naive + CPolynomial.Raw.ModContext.naive modulus acc current + else + acc) + (current := + CPolynomial.Raw.mulModWith CPolynomial.Raw.MulContext.naive + CPolynomial.Raw.ModContext.naive modulus current current) + rw [ih'] + rw [raw_eval_mulMod_naive_eq_mul hmod] + by_cases hodd : (n + 1) % 2 == 1 + · simp [hodd, raw_eval_mulMod_naive_eq_mul hmod] + have hoddNat : (n + 1) % 2 = 1 := by + simpa using hodd + have hpow : + (Raw.eval a current * Raw.eval a current) ^ ((n + 1) / 2) = + Raw.eval a current ^ n := by + rw [mul_pow, ← pow_add] + congr 1 + omega + simp [hpow, pow_succ] + ring + · simp [hodd] + have hoddNat : (n + 1) % 2 = 0 := by + simpa using hodd + have hpow : + (Raw.eval a current * Raw.eval a current) ^ ((n + 1) / 2) = + Raw.eval a current ^ (n + 1) := by + rw [mul_pow, ← pow_add] + congr 1 + omega + simp [hpow] + +private theorem raw_eval_one {F : Type*} [Field F] (a : F) : + (1 : CPolynomial.Raw F).eval a = 1 := by + rw [← CPolynomial.Raw.eval_toPoly_eq_eval] + simp + +private theorem raw_eval_X {F : Type*} [Field F] (a : F) : + (CPolynomial.Raw.X : CPolynomial.Raw F).eval a = a := by + rw [← CPolynomial.Raw.eval_toPoly_eq_eval] + simp [CPolynomial.Raw.toPoly_X] + +private theorem raw_eval_powMod_naive_X {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {modulus : CPolynomial.Raw F} {a : F} (hmod : modulus.eval a = 0) (q : Nat) : + (CPolynomial.Raw.powModWith CPolynomial.Raw.MulContext.naive + CPolynomial.Raw.ModContext.naive modulus CPolynomial.Raw.X q).eval a = + a ^ q := by + unfold CPolynomial.Raw.powModWith + by_cases hzero : modulus.trim = (#[] : CPolynomial.Raw F) + · simp [hzero] + rw [raw_eval_powModBinaryAux_naive hmod, raw_eval_one, raw_eval_X] + simp + · simp [hzero] + have hroot : (CPolynomial.Raw.monicNormalize modulus).eval a = 0 := + CPolynomial.Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero hmod + have hmodBy := CPolynomial.Raw.eval_modByMonic_eq_self_of_eval_eq_zero + (1 : CPolynomial.Raw F) (CPolynomial.Raw.monicNormalize modulus) hroot + rw [raw_eval_powModBinaryAux_naive hmod] + simp [CPolynomial.Raw.ModContext.naive, hmodBy, raw_eval_one, raw_eval_X] + +private theorem raw_eval_xMod_naive {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {modulus : CPolynomial.Raw F} {a : F} (hmod : modulus.eval a = 0) : + (CPolynomial.Raw.xModWith CPolynomial.Raw.ModContext.naive modulus).eval a = a := by + unfold CPolynomial.Raw.xModWith CPolynomial.Raw.ModContext.naive + by_cases hzero : modulus.trim = (#[] : CPolynomial.Raw F) + · simp [hzero, raw_eval_X] + · simp [hzero] + have hroot : (CPolynomial.Raw.monicNormalize modulus).eval a = 0 := + CPolynomial.Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero hmod + have hmodBy := CPolynomial.Raw.eval_modByMonic_eq_self_of_eval_eq_zero + (CPolynomial.Raw.X : CPolynomial.Raw F) (CPolynomial.Raw.monicNormalize modulus) hroot + rw [hmodBy, raw_eval_X] + +private theorem raw_eval_sub {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (a : F) (p q : CPolynomial.Raw F) : + (p - q).eval a = p.eval a - q.eval a := by + rw [← CPolynomial.Raw.eval_toPoly_eq_eval, CPolynomial.Raw.toPoly_sub, + Polynomial.eval_sub, CPolynomial.Raw.eval_toPoly_eq_eval, + CPolynomial.Raw.eval_toPoly_eq_eval] + +private theorem raw_eval_xPowSubXMod_naive {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (ctx : FiniteFieldContext F) {modulus : CPolynomial.Raw F} {a : F} + (hmod : modulus.eval a = 0) : + (CPolynomial.Raw.xPowSubXModWith CPolynomial.Raw.MulContext.naive + CPolynomial.Raw.ModContext.naive ctx.q modulus).eval a = 0 := by + unfold CPolynomial.Raw.xPowSubXModWith + rw [raw_eval_sub, raw_eval_powMod_naive_X hmod, raw_eval_xMod_naive hmod, + ctx.frobenius_fixed a] + simp + +private theorem raw_monicNormalize_trim {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (p : CPolynomial.Raw F) : + (CPolynomial.Raw.monicNormalize p).trim = CPolynomial.Raw.monicNormalize p := by + unfold CPolynomial.Raw.monicNormalize + let q := p.trim + by_cases hzero : q == (0 : CPolynomial.Raw F) + · have hzeroEq : q = 0 := by + simpa using hzero + have hpzeroEq : p.trim = (#[] : CPolynomial.Raw F) := by + simpa [q] using hzeroEq + simp [hpzeroEq, CPolynomial.Raw.Trim.canonical_empty] + · have hzeroNe : q ≠ (0 : CPolynomial.Raw F) := by + intro hq + exact hzero (by simp [hq]) + have hpzeroNe : ¬p.trim = (#[] : CPolynomial.Raw F) := by + intro hp + exact hzeroNe (by simpa [q] using hp) + simp [hpzeroNe] + change (CPolynomial.Raw.mk + (Array.map (fun r ↦ (p.trim).leadingCoeff⁻¹ * r) p.trim)).trim = + CPolynomial.Raw.mk (Array.map (fun r ↦ (p.trim).leadingCoeff⁻¹ * r) p.trim) + apply CPolynomial.Raw.Trim.non_zero_map (fun r ↦ q.leadingCoeff⁻¹ * r) + · intro r hr + apply mul_eq_zero.mp at hr + rcases hr with hinv | hr + · have hlead0 : q.leadingCoeff = 0 := by + exact inv_eq_zero.mp hinv + have hcanon : q.trim = q := by + simpa [q] using CPolynomial.Raw.Trim.trim_twice p + have hcrit := (CPolynomial.Raw.Trim.trim_eq_iff_size_eq_zero_or_getLastD_ne_zero + (p := q)).mp hcanon + rcases hcrit with hsize | hlast + · have hqempty : q = (#[] : CPolynomial.Raw F) := by + apply Array.eq_empty_of_size_eq_zero + exact hsize + exact False.elim (hzero (by simp [hqempty])) + · unfold CPolynomial.Raw.leadingCoeff at hlead0 + rw [hcanon] at hlead0 + exact (hlast hlead0).elim + · exact hr + · simpa [q] using CPolynomial.Raw.Trim.trim_twice p + +private theorem raw_eval_mulModWith_eq_mul {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + {modulus p q : CPolynomial.Raw F} {a : F} + (hmod : modulus.eval a = 0) : + (CPolynomial.Raw.mulModWith M D modulus p q).eval a = + p.eval a * q.eval a := by + unfold CPolynomial.Raw.mulModWith + by_cases hzero : modulus.trim = (#[] : CPolynomial.Raw F) + · simp [hzero, M.mul_eq_mul, CPolynomial.Raw.eval_mul] + · simp [hzero] + have hproductTrim : (M.mul p q).trim = M.mul p q := by + rw [M.mul_eq_mul] + exact CPolynomial.Raw.mul_is_trimmed p q + rw [D.modByMonic_eq_modByMonic _ _ hproductTrim (raw_monicNormalize_trim modulus)] + have hroot : (CPolynomial.Raw.monicNormalize modulus).eval a = 0 := + CPolynomial.Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero hmod + have hmodBy := CPolynomial.Raw.eval_modByMonic_eq_self_of_eval_eq_zero + (M.mul p q) (CPolynomial.Raw.monicNormalize modulus) hroot + rw [hmodBy, M.mul_eq_mul, CPolynomial.Raw.eval_mul] + +private theorem raw_eval_powModBinaryAuxWith {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + {modulus acc current : CPolynomial.Raw F} {a : F} + (hmod : modulus.eval a = 0) : + ∀ n, + (CPolynomial.Raw.powModBinaryAuxWith M D modulus n acc current).eval a = + acc.eval a * current.eval a ^ n := by + intro n + induction n using Nat.strongRecOn generalizing acc current with + | ind n ih => + cases n with + | zero => + simp [CPolynomial.Raw.powModBinaryAuxWith] + | succ n => + rw [CPolynomial.Raw.powModBinaryAuxWith] + have ih' := ih ((n + 1) / 2) + (Nat.div_lt_self (Nat.succ_pos n) (by decide)) + (acc := + if (n + 1) % 2 == 1 then + CPolynomial.Raw.mulModWith M D modulus acc current + else + acc) + (current := CPolynomial.Raw.mulModWith M D modulus current current) + rw [ih'] + rw [raw_eval_mulModWith_eq_mul M D hmod] + by_cases hodd : (n + 1) % 2 == 1 + · simp [hodd, raw_eval_mulModWith_eq_mul M D hmod] + have hoddNat : (n + 1) % 2 = 1 := by + simpa using hodd + have hpow : + (Raw.eval a current * Raw.eval a current) ^ ((n + 1) / 2) = + Raw.eval a current ^ n := by + rw [mul_pow, ← pow_add] + congr 1 + omega + simp [hpow, pow_succ] + ring + · simp [hodd] + have hoddNat : (n + 1) % 2 = 0 := by + simpa using hodd + have hpow : + (Raw.eval a current * Raw.eval a current) ^ ((n + 1) / 2) = + Raw.eval a current ^ (n + 1) := by + rw [mul_pow, ← pow_add] + congr 1 + omega + simp [hpow] + +private theorem raw_eval_powModWith_X {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + {modulus : CPolynomial.Raw F} {a : F} (hmod : modulus.eval a = 0) (q : Nat) : + (CPolynomial.Raw.powModWith M D modulus CPolynomial.Raw.X q).eval a = + a ^ q := by + unfold CPolynomial.Raw.powModWith + by_cases hzero : modulus.trim = (#[] : CPolynomial.Raw F) + · simp [hzero] + rw [raw_eval_powModBinaryAuxWith M D hmod, raw_eval_one, raw_eval_X] + simp + · simp [hzero] + have hroot : (CPolynomial.Raw.monicNormalize modulus).eval a = 0 := + CPolynomial.Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero hmod + have hmodBy := CPolynomial.Raw.eval_modByMonic_eq_self_of_eval_eq_zero + (1 : CPolynomial.Raw F) (CPolynomial.Raw.monicNormalize modulus) hroot + rw [raw_eval_powModBinaryAuxWith M D hmod] + have hOneTrim : (1 : CPolynomial.Raw F).trim = 1 := by + change CPolynomial.Raw.trim (#[] |>.push (1 : F)) = (#[] |>.push (1 : F)) + apply CPolynomial.Raw.Trim.push_trim + simp + rw [D.modByMonic_eq_modByMonic _ _ hOneTrim (raw_monicNormalize_trim modulus)] + rw [hmodBy, raw_eval_one, raw_eval_X] + simp + +/-- Raw modular exponentiation evaluates as ordinary exponentiation at roots of +the modulus. -/ +theorem raw_eval_powModWith_eq_pow {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + {modulus base : CPolynomial.Raw F} {a : F} + (hmod : modulus.eval a = 0) (q : Nat) : + (CPolynomial.Raw.powModWith M D modulus base q).eval a = + base.eval a ^ q := by + unfold CPolynomial.Raw.powModWith + by_cases hzero : modulus.trim = (#[] : CPolynomial.Raw F) + · simp [hzero] + rw [raw_eval_powModBinaryAuxWith M D hmod, raw_eval_one] + simp + · simp [hzero] + have hroot : (CPolynomial.Raw.monicNormalize modulus).eval a = 0 := + CPolynomial.Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero hmod + have hmodBy := CPolynomial.Raw.eval_modByMonic_eq_self_of_eval_eq_zero + (1 : CPolynomial.Raw F) (CPolynomial.Raw.monicNormalize modulus) hroot + rw [raw_eval_powModBinaryAuxWith M D hmod] + have hOneTrim : (1 : CPolynomial.Raw F).trim = 1 := by + change CPolynomial.Raw.trim (#[] |>.push (1 : F)) = (#[] |>.push (1 : F)) + apply CPolynomial.Raw.Trim.push_trim + simp + rw [D.modByMonic_eq_modByMonic _ _ hOneTrim (raw_monicNormalize_trim modulus)] + rw [hmodBy, raw_eval_one] + simp + +private theorem raw_eval_xModWith {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (D : CPolynomial.Raw.ModContext F) + {modulus : CPolynomial.Raw F} {a : F} (hmod : modulus.eval a = 0) : + (CPolynomial.Raw.xModWith D modulus).eval a = a := by + unfold CPolynomial.Raw.xModWith + by_cases hzero : modulus.trim = (#[] : CPolynomial.Raw F) + · simp [hzero, raw_eval_X] + · simp [hzero] + have hroot : (CPolynomial.Raw.monicNormalize modulus).eval a = 0 := + CPolynomial.Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero hmod + have hmodBy := CPolynomial.Raw.eval_modByMonic_eq_self_of_eval_eq_zero + (CPolynomial.Raw.X : CPolynomial.Raw F) (CPolynomial.Raw.monicNormalize modulus) hroot + have hXTrim : (CPolynomial.Raw.X : CPolynomial.Raw F).trim = CPolynomial.Raw.X := by + exact CPolynomial.Raw.X_canonical + rw [D.modByMonic_eq_modByMonic _ _ hXTrim (raw_monicNormalize_trim modulus)] + rw [hmodBy, raw_eval_X] + +private theorem raw_eval_xPowSubXModWith {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : FiniteFieldContext F) {modulus : CPolynomial.Raw F} {a : F} + (hmod : modulus.eval a = 0) : + (CPolynomial.Raw.xPowSubXModWith M D ctx.q modulus).eval a = 0 := by + unfold CPolynomial.Raw.xPowSubXModWith + rw [raw_eval_sub, raw_eval_powModWith_X M D hmod, raw_eval_xModWith D hmod, + ctx.frobenius_fixed a] + simp + +/-- Every root of `p` is a root of the finite-field root product. -/ +theorem finiteFieldRootProductWith_complete {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : FiniteFieldContext F) {p : CPolynomial F} {a : F} + (hp : p ≠ 0) (hroot : CPolynomial.eval a p = 0) : + CPolynomial.eval a (finiteFieldRootProductWith M D ctx p) = 0 := by + unfold finiteFieldRootProductWith + change + CPolynomial.Raw.eval a + (CPolynomial.Raw.Roots.FiniteField.finiteFieldRootProductWith M D + ctx p.val).trim = + 0 + rw [CPolynomial.Raw.eval_trim_eq_eval] + unfold CPolynomial.Raw.Roots.FiniteField.finiteFieldRootProductWith + by_cases hpempty : p.val = (#[] : CPolynomial.Raw F) + · have hp0 : p = 0 := by + apply CPolynomial.ext + simpa using hpempty + exact (hp hp0).elim + · simp [hpempty, CPolynomial.trim_eq] + have hmonicRoot : (CPolynomial.Raw.monicNormalize p.val).eval a = 0 := + CPolynomial.Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero hroot + unfold CPolynomial.Raw.gcdMonic + exact Raw.eval_gcdMonicWithFuel_eq_zero_of_left_right + _ (CPolynomial.Raw.monicNormalize p.val) + (CPolynomial.Raw.xPowSubXModWith M D ctx.q (CPolynomial.Raw.monicNormalize p.val)) + hmonicRoot (raw_eval_xPowSubXModWith M D ctx hmonicRoot) + +/-- Every root of `p` is a root of the finite-field root product. -/ +theorem finiteFieldRootProduct_complete {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (ctx : FiniteFieldContext F) {p : CPolynomial F} {a : F} + (hp : p ≠ 0) (hroot : CPolynomial.eval a p = 0) : + CPolynomial.eval a (finiteFieldRootProduct ctx p) = 0 := by + unfold finiteFieldRootProduct finiteFieldRootProductWith + change + CPolynomial.Raw.eval a + (CPolynomial.Raw.Roots.FiniteField.finiteFieldRootProductWith + CPolynomial.Raw.MulContext.naive CPolynomial.Raw.ModContext.naive + ctx p.val).trim = + 0 + rw [CPolynomial.Raw.eval_trim_eq_eval] + unfold CPolynomial.Raw.Roots.FiniteField.finiteFieldRootProductWith + by_cases hpempty : p.val = (#[] : CPolynomial.Raw F) + · have hp0 : p = 0 := by + apply CPolynomial.ext + simpa using hpempty + exact (hp hp0).elim + · simp [hpempty, CPolynomial.trim_eq] + have hmonicRoot : (CPolynomial.Raw.monicNormalize p.val).eval a = 0 := + CPolynomial.Raw.eval_monicNormalize_eq_zero_of_eval_eq_zero hroot + unfold CPolynomial.Raw.gcdMonic + exact Raw.eval_gcdMonicWithFuel_eq_zero_of_left_right + _ (CPolynomial.Raw.monicNormalize p.val) + (CPolynomial.Raw.xPowSubXModWith CPolynomial.Raw.MulContext.naive + CPolynomial.Raw.ModContext.naive ctx.q (CPolynomial.Raw.monicNormalize p.val)) + hmonicRoot (raw_eval_xPowSubXMod_naive ctx hmonicRoot) + +/-- Every validated root extracted from the root product is a root of `p`. -/ +theorem finiteFieldRootProduct_validated_sound {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (_ctx : FiniteFieldContext F) {p : CPolynomial F} + {candidates : Array F} {a : F} + (h : a ∈ (CPolynomial.validateRootCandidates p candidates).toList) : + CPolynomial.eval a p = 0 := by + exact CPolynomial.mem_validateRootCandidates_eval_eq_zero h + +/-- Returned finite-field roots are roots of the original polynomial. -/ +theorem rootsInFiniteFieldWith_sound {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : FiniteFieldContext F) (splitter : LinearFactorProductSplitter F) + {p : CPolynomial F} {a : F} + (h : a ∈ (rootsInFiniteFieldWith M D ctx splitter p).toList) : + CPolynomial.eval a p = 0 := by + rw [rootsInFiniteFieldWith] at h + split at h + · simp at h + · split at h + · simp at h + · split at h + · exact mem_rootsFromLinearFactors_eval_eq_zero h + · exact mem_rootsFromLinearFactors_eval_eq_zero h + +/-- Returned finite-field roots are roots of the original polynomial. -/ +theorem rootsInFiniteField_sound {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (ctx : FiniteFieldContext F) (splitter : LinearFactorProductSplitter F) + {p : CPolynomial F} {a : F} + (h : a ∈ (rootsInFiniteField ctx splitter p).toList) : + CPolynomial.eval a p = 0 := by + exact rootsInFiniteFieldWith_sound + (M := CPolynomial.Raw.MulContext.naive) (D := CPolynomial.Raw.ModContext.naive) + ctx splitter h + +private theorem eq_zero_of_size_le_one_root {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p : CPolynomial F} {a : F} + (hsize : p.val.size ≤ 1) (hroot : CPolynomial.eval a p = 0) : + p = 0 := by + rcases p with ⟨⟨xs⟩, hcanon⟩ + cases xs with + | nil => + rfl + | cons x xs => + cases xs with + | nil => + simp [CPolynomial.eval] at hroot + have hxne := hcanon (by simp) + simp [Array.getLast] at hxne + exact (hxne hroot).elim + | cons _ _ => + simp at hsize + +private theorem linear_candidate_self_of_size_two_root {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p : CPolynomial F} {a : F} + (hsize : p.val.size = 2) (hcoeff : p.coeff 1 ≠ 0) + (hroot : CPolynomial.eval a p = 0) : + IsLinearRootFactorCandidate p a := by + rcases p with ⟨⟨xs⟩, hcanon⟩ + cases xs with + | nil => + simp at hsize + | cons x xs => + cases xs with + | nil => + simp at hsize + | cons y xs => + cases xs with + | nil => + simp [IsLinearRootFactorCandidate, IsLinearFactor, CPolynomial.eval, + CPolynomial.coeff, CPolynomial.Raw.coeff] at hroot hcoeff ⊢ + exact ⟨hcoeff, hroot⟩ + | cons _ _ => + simp at hsize + +/-- Every explicit `X - a` factor is represented as a nonconstant linear factor. -/ +theorem linearFactor_isLinearFactor {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (a : F) : IsLinearFactor (CPolynomial.linearFactor a) := by + constructor + · rw [CPolynomial.linearFactor] + calc + ((CPolynomial.C (-a) : CPolynomial F) + CPolynomial.X).val.size + = (CPolynomial.Raw.addRaw (CPolynomial.C (-a)).val CPolynomial.X.val).trim.size := by + rfl + _ ≤ (CPolynomial.Raw.addRaw (CPolynomial.C (-a)).val CPolynomial.X.val).size := + CPolynomial.Raw.Trim.size_le_size _ + _ = max (CPolynomial.C (-a)).val.size CPolynomial.X.val.size := + CPolynomial.Raw.add_size + _ ≤ 2 := by + have hC : (CPolynomial.C (-a) : CPolynomial F).val.size ≤ 1 := by + unfold CPolynomial.C CPolynomial.Raw.C + exact CPolynomial.Raw.Trim.size_le_size _ + have hX : (CPolynomial.X : CPolynomial F).val.size ≤ 2 := by + rfl + exact max_le (le_trans hC (by omega)) hX + · rw [CPolynomial.linearFactor, CPolynomial.coeff_add, CPolynomial.coeff_C] + simp [CPolynomial.X, CPolynomial.Raw.X, CPolynomial.coeff, CPolynomial.Raw.coeff] + +/-- The represented-linear recognizer accepts every explicit `X - a` factor. -/ +theorem linearFactor_isRepresentedLinearFactor {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (a : F) : isRepresentedLinearFactor (CPolynomial.linearFactor a) = true := by + unfold isRepresentedLinearFactor + have hlin := linearFactor_isLinearFactor a + have hcoeff : ¬(CPolynomial.linearFactor a).val[1]?.getD 0 = 0 := by + simpa [CPolynomial.coeff, CPolynomial.Raw.coeff] using hlin.2 + simp [hlin.1, hcoeff, CPolynomial.coeff, CPolynomial.Raw.coeff] + +/-- The explicit `X - a` factor represents the root `a`. -/ +theorem linearFactor_isRootFactorCandidate {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (a : F) : IsLinearRootFactorCandidate (CPolynomial.linearFactor a) a := by + constructor + · exact linearFactor_isLinearFactor a + · have h0 : (CPolynomial.linearFactor a).coeff 0 = -a := by + rw [CPolynomial.linearFactor, CPolynomial.coeff_add, CPolynomial.coeff_C] + simp [CPolynomial.X, CPolynomial.Raw.X, CPolynomial.coeff, CPolynomial.Raw.coeff] + have h1 : (CPolynomial.linearFactor a).coeff 1 = 1 := by + rw [CPolynomial.linearFactor, CPolynomial.coeff_add, CPolynomial.coeff_C] + simp [CPolynomial.X, CPolynomial.Raw.X, CPolynomial.coeff, CPolynomial.Raw.coeff] + rw [h0, h1] + simp + +/-- If the represented-linear recognizer accepts `p`, then any root of `p` is its root candidate. -/ +theorem representedLinearFactor_candidate_of_root {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p : CPolynomial F} {a : F} + (hlin : isRepresentedLinearFactor p = true) + (hroot : CPolynomial.eval a p = 0) : + IsLinearRootFactorCandidate p a := by + rcases p with ⟨⟨xs⟩, hcanon⟩ + cases xs with + | nil => + simp [isRepresentedLinearFactor, CPolynomial.coeff, CPolynomial.Raw.coeff] at hlin + | cons x xs => + cases xs with + | nil => + simp [isRepresentedLinearFactor, CPolynomial.coeff, CPolynomial.Raw.coeff] at hlin + | cons y xs => + cases xs with + | nil => + simp [isRepresentedLinearFactor, IsLinearRootFactorCandidate, IsLinearFactor, + CPolynomial.eval, CPolynomial.coeff, CPolynomial.Raw.coeff] at hroot hlin ⊢ + exact ⟨hlin, hroot⟩ + | cons _ _ => + simp [isRepresentedLinearFactor] at hlin + +/-- Evaluation of a modular `X^n` at a root of the modulus agrees with `a^n`. -/ +theorem eval_xPowModWith_eq_pow {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + {p : CPolynomial F} {a : F} (hroot : CPolynomial.eval a p = 0) (n : Nat) : + CPolynomial.eval a (xPowModWith M D p n) = a ^ n := by + unfold xPowModWith + change CPolynomial.Raw.eval a + (CPolynomial.Raw.powModWith M D p.val CPolynomial.Raw.X n).trim = a ^ n + rw [CPolynomial.Raw.eval_trim_eq_eval] + exact raw_eval_powModWith_X M D hroot n + +/-- Every root of a nonzero polynomial is returned by the finite-field backend. -/ +theorem rootsInFiniteFieldWith_complete {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : FiniteFieldContext F) (splitter : LinearFactorProductSplitter F) + (splitterValid : + ∀ {p : CPolynomial F}, p ≠ 0 → + splitter.validInput ctx.q (finiteFieldRootProductWith M D ctx p)) + {p : CPolynomial F} {a : F} + (hp : p ≠ 0) (hroot : CPolynomial.eval a p = 0) : + a ∈ (rootsInFiniteFieldWith M D ctx splitter p).toList := by + rw [rootsInFiniteFieldWith] + split + · rename_i hzero + have hp0 : p = 0 := by + simpa using hzero + exact (hp hp0).elim + · split + · rename_i _ hsmall + exact (hp (eq_zero_of_size_le_one_root hsmall hroot)).elim + · split + · rename_i _ _ hlinear + simp at hlinear + have hcoeff : p.coeff 1 ≠ 0 := by + simpa [CPolynomial.coeff, CPolynomial.Raw.coeff] using hlinear.2 + exact mem_rootsFromLinearFactors_of_candidate (p := p) (factors := #[p]) + (factor := p) (by simp) + (linear_candidate_self_of_size_two_root hlinear.1 hcoeff hroot) hroot + · have hprodRoot : CPolynomial.eval a (finiteFieldRootProductWith M D ctx p) = 0 := + finiteFieldRootProductWith_complete M D ctx hp hroot + have hprodNe : finiteFieldRootProductWith M D ctx p ≠ 0 := + finiteFieldRootProductWith_ne_zero_of_ne_zero M D ctx hp + have hprodValid : + splitter.validInput ctx.q (finiteFieldRootProductWith M D ctx p) := by + exact splitterValid hp + rcases splitter.complete ctx.q (finiteFieldRootProductWith M D ctx p) a + hprodValid hprodNe hprodRoot with + ⟨factor, hmem, hcand⟩ + exact mem_rootsFromLinearFactors_of_candidate (p := p) + (factors := splitter.splitLinearFactors ctx.q (finiteFieldRootProductWith M D ctx p)) + (factor := factor) hmem hcand hroot + +/-- Every root of a nonzero polynomial is returned by the finite-field backend. -/ +theorem rootsInFiniteField_complete {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (ctx : FiniteFieldContext F) (splitter : LinearFactorProductSplitter F) + (splitterValid : + ∀ {p : CPolynomial F}, p ≠ 0 → + splitter.validInput ctx.q (finiteFieldRootProduct ctx p)) + {p : CPolynomial F} {a : F} + (hp : p ≠ 0) (hroot : CPolynomial.eval a p = 0) : + a ∈ (rootsInFiniteField ctx splitter p).toList := by + exact rootsInFiniteFieldWith_complete + (M := CPolynomial.Raw.MulContext.naive) (D := CPolynomial.Raw.ModContext.naive) + ctx splitter (by + intro p hp + exact splitterValid hp) hp hroot + +/-- The complete executable finite-field root pipeline is sound and complete +for nonzero inputs under the finite-field and splitter contracts. -/ +theorem rootsInFiniteField_spec {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (ctx : FiniteFieldContext F) (splitter : LinearFactorProductSplitter F) + (splitterValid : + ∀ {p : CPolynomial F}, p ≠ 0 → + splitter.validInput ctx.q (finiteFieldRootProduct ctx p)) + {p : CPolynomial F} {a : F} (hp : p ≠ 0) : + a ∈ (rootsInFiniteField ctx splitter p).toList ↔ + CPolynomial.eval a p = 0 := by + constructor + · intro h + exact rootsInFiniteField_sound ctx splitter h + · intro h + exact rootsInFiniteField_complete ctx splitter splitterValid hp h + +end FiniteField + +end Roots + +end CPolynomial + +end CompPoly diff --git a/CompPoly/Univariate/Roots/Enumeration.lean b/CompPoly/Univariate/Roots/Enumeration.lean new file mode 100644 index 00000000..15490940 --- /dev/null +++ b/CompPoly/Univariate/Roots/Enumeration.lean @@ -0,0 +1,124 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Roots.Correctness + +/-! +# Exhaustive Finite-Field Root Enumeration + +Reusable lazy field-enumeration contexts for finite-field root search. The +context stores an indexing function rather than an array of all elements; array +inputs are adapted through `fieldEnumerationOfArray` for tests and small +callers. +-/ + +namespace CompPoly + +namespace CPolynomial + +namespace Roots + +namespace FiniteField + +/-- A lazy complete enumeration of a finite field. -/ +structure FieldEnumeration (F : Type*) where + size : Nat + elem : Fin size → F + complete : ∀ a : F, ∃ i : Fin size, elem i = a + +/-- An array contains every field element. Duplicate entries are allowed. -/ +def ContainsAllFieldElements {F : Type*} (elements : Array F) : Prop := + ∀ a : F, a ∈ elements.toList + +/-- Adapt an explicit element array to a lazy enumeration context. -/ +def fieldEnumerationOfArray {F : Type*} (elements : Array F) + (hElements : ContainsAllFieldElements elements) : + FieldEnumeration F where + size := elements.size + elem i := elements[i] + complete := by + intro a + rcases List.mem_iff_get.mp (hElements a) with ⟨i, hi⟩ + refine ⟨⟨i.val, by simpa only [Array.length_toList] using i.isLt⟩, ?_⟩ + simpa using hi + +/-- Roots by exhaustive evaluation over a lazy field enumeration. -/ +def rootsInFieldByEnumeration {F : Type*} [Semiring F] [BEq F] + (enumeration : FieldEnumeration F) (p : CPolynomial F) : Array F := + (Array.ofFn enumeration.elem).filter fun a ↦ CPolynomial.evalHorner a p == 0 + +/-- Exhaustive enumeration only returns actual roots. -/ +theorem rootsInFieldByEnumeration_sound {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {enumeration : FieldEnumeration F} {p : CPolynomial F} {a : F} + (h : a ∈ (rootsInFieldByEnumeration enumeration p).toList) : + CPolynomial.eval a p = 0 := by + rw [rootsInFieldByEnumeration] at h + simp at h + simpa [CPolynomial.eval_horner_eq_eval] using h.2 + +/-- Complete enumeration finds every root. -/ +theorem rootsInFieldByEnumeration_complete {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (enumeration : FieldEnumeration F) {p : CPolynomial F} {a : F} + (hroot : CPolynomial.eval a p = 0) : + a ∈ (rootsInFieldByEnumeration enumeration p).toList := by + rw [rootsInFieldByEnumeration] + rcases enumeration.complete a with ⟨i, hi⟩ + simpa [CPolynomial.eval_horner_eq_eval, hroot] using ⟨i, hi⟩ + +/-- Linear factors for every enumerated root of `p`. -/ +def enumeratedLinearFactors {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (enumeration : FieldEnumeration F) (p : CPolynomial F) : + Array (CPolynomial F) := + (rootsInFieldByEnumeration enumeration p).map CPolynomial.linearFactor + +/-- Every factor emitted by exhaustive enumeration is represented linear. -/ +theorem enumeratedLinearFactors_sound {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {enumeration : FieldEnumeration F} {p factor : CPolynomial F} + (h : factor ∈ (enumeratedLinearFactors enumeration p).toList) : + IsLinearFactor factor := by + rw [enumeratedLinearFactors] at h + simp at h + rcases h with ⟨a, _hmem, rfl⟩ + exact linearFactor_isLinearFactor a + +/-- Exhaustive enumeration emits the linear factor for every root. -/ +theorem enumeratedLinearFactors_complete {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (enumeration : FieldEnumeration F) {p : CPolynomial F} {a : F} + (hroot : CPolynomial.eval a p = 0) : + ∃ factor, + factor ∈ (enumeratedLinearFactors enumeration p).toList ∧ + IsLinearRootFactorCandidate factor a := by + refine ⟨CPolynomial.linearFactor a, ?_, linearFactor_isRootFactorCandidate a⟩ + rw [enumeratedLinearFactors] + simpa using + (List.mem_map.mpr + ⟨a, rootsInFieldByEnumeration_complete enumeration hroot, rfl⟩) + +/-- Exhaustive enumeration packaged as a linear-factor product splitter. -/ +def enumeratingLinearFactorProductSplitter {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (enumeration : FieldEnumeration F) : + LinearFactorProductSplitter F where + splitLinearFactors := fun _ p ↦ enumeratedLinearFactors enumeration p + validInput := fun _ _ ↦ True + sound := by + intro _q p factor h + exact enumeratedLinearFactors_sound h + complete := by + intro _q p a _hvalid _hp hroot + exact enumeratedLinearFactors_complete enumeration hroot + +end FiniteField + +end Roots + +end CPolynomial + +end CompPoly diff --git a/CompPoly/Univariate/Roots/Extraction.lean b/CompPoly/Univariate/Roots/Extraction.lean new file mode 100644 index 00000000..1dc3a51e --- /dev/null +++ b/CompPoly/Univariate/Roots/Extraction.lean @@ -0,0 +1,41 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Basic + +/-! +# Root Extraction Helpers + +Candidate extraction, validation, and deduplication helpers for univariate root +finding. +-/ + +namespace CompPoly + +namespace CPolynomial + +/-- Extract the root of a represented nonconstant linear factor. -/ +def linearRootOfFactor? {F : Type*} [Field F] [BEq F] + (factor : CPolynomial F) : Option F := + if factor.val.size ≤ 2 && !(factor.coeff 1 == 0) then + some (-(factor.coeff 0) / factor.coeff 1) + else + none + +/-- Keep only candidates that are actual roots of the original polynomial. -/ +def validateRootCandidates {F : Type*} [Field F] [BEq F] + (p : CPolynomial F) (candidates : Array F) : Array F := + candidates.filter fun a ↦ CPolynomial.evalHorner a p == 0 + +/-- Extract, validate, and deduplicate roots from a list of linear factors. -/ +def rootsFromLinearFactors {F : Type*} [Field F] [BEq F] + (p : CPolynomial F) (factors : Array (CPolynomial F)) : Array F := + let candidates := (factors.toList.filterMap linearRootOfFactor?).toArray + (validateRootCandidates p candidates).eraseDups + +end CPolynomial + +end CompPoly diff --git a/CompPoly/Univariate/Roots/RootProduct.lean b/CompPoly/Univariate/Roots/RootProduct.lean new file mode 100644 index 00000000..794034be --- /dev/null +++ b/CompPoly/Univariate/Roots/RootProduct.lean @@ -0,0 +1,643 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Raw.Modular +import CompPoly.Univariate.EuclideanAlgorithm +import CompPoly.Univariate.Roots.Context + +/-! +# Finite-Field Root Products + +Executable construction of `gcd(p, X^q - X)` as +`gcd(p, (X^q mod p) - (X mod p))`, so large finite fields never materialize the +dense polynomial `X^q - X`. +-/ + +namespace CompPoly + +namespace CPolynomial + +namespace Raw + +namespace Roots + +namespace FiniteField + +/-- Raw squarefree product of the linear factors of `p` whose roots lie in the field. -/ +def finiteFieldRootProductWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : CPolynomial.Roots.FiniteField.FiniteFieldContext F) + (p : CPolynomial.Raw F) : CPolynomial.Raw F := + if p.trim == 0 then + 0 + else + let pMonic := CPolynomial.Raw.monicNormalize p + CPolynomial.Raw.gcdMonic pMonic (CPolynomial.Raw.xPowSubXModWith M D ctx.q pMonic) + +/-- Raw squarefree product using the default raw backends. -/ +def finiteFieldRootProduct {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (ctx : CPolynomial.Roots.FiniteField.FiniteFieldContext F) + (p : CPolynomial.Raw F) : CPolynomial.Raw F := + finiteFieldRootProductWith CPolynomial.Raw.MulContext.naive + CPolynomial.Raw.ModContext.naive ctx p + +end FiniteField + +end Roots + +end Raw + +namespace Roots + +namespace FiniteField + +/-- The squarefree product of the linear factors of `p` whose roots lie in the field. -/ +def finiteFieldRootProductWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : FiniteFieldContext F) (p : CPolynomial F) : CPolynomial F := + CPolynomial.ofArray (CPolynomial.Raw.Roots.FiniteField.finiteFieldRootProductWith M D ctx p.val) + +/-- The squarefree product of the linear factors of `p` using the default raw backends. -/ +def finiteFieldRootProduct {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (ctx : FiniteFieldContext F) (p : CPolynomial F) : CPolynomial F := + finiteFieldRootProductWith CPolynomial.Raw.MulContext.naive + CPolynomial.Raw.ModContext.naive ctx p + +private theorem raw_monicNormalize_ne_zero_of_trim_ne_zero {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p : CPolynomial.Raw F} (hp : p.trim ≠ (0 : CPolynomial.Raw F)) : + CPolynomial.Raw.monicNormalize p ≠ 0 := by + unfold CPolynomial.Raw.monicNormalize + by_cases hzero : p.trim == (0 : CPolynomial.Raw F) + · have hzeroEq : p.trim = 0 := by + simpa using hzero + exact (hp hzeroEq).elim + · have hpzeroNe : ¬p.trim = (#[] : CPolynomial.Raw F) := by + intro hpzero + exact hzero (by simp [hpzero]) + rw [if_neg hzero] + intro hsmul + have hsize_smul : + (CPolynomial.Raw.smul (p.trim.leadingCoeff)⁻¹ p.trim).size = p.trim.size := by + simp [CPolynomial.Raw.smul] + change CPolynomial.Raw.smul (p.trim.leadingCoeff)⁻¹ p.trim = 0 at hsmul + rw [hsmul] at hsize_smul + have hsize : p.trim.size = 0 := by + simpa using hsize_smul.symm + have hpempty : p.trim = (#[] : CPolynomial.Raw F) := by + apply Array.eq_empty_of_size_eq_zero + exact hsize + exact hp (by simpa using hpempty) + +private theorem raw_monicNormalize_trim {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (p : CPolynomial.Raw F) : + (CPolynomial.Raw.monicNormalize p).trim = CPolynomial.Raw.monicNormalize p := by + unfold CPolynomial.Raw.monicNormalize + let q := p.trim + by_cases hzero : q == (0 : CPolynomial.Raw F) + · have hzeroEq : q = 0 := by + simpa using hzero + have hpzeroEq : p.trim = (#[] : CPolynomial.Raw F) := by + simpa [q] using hzeroEq + simp [hpzeroEq, CPolynomial.Raw.Trim.canonical_empty] + · have hzeroNe : q ≠ (0 : CPolynomial.Raw F) := by + intro hq + exact hzero (by simp [hq]) + have hpzeroNe : ¬p.trim = (#[] : CPolynomial.Raw F) := by + intro hp + exact hzeroNe (by simpa [q] using hp) + simp [hpzeroNe] + change (CPolynomial.Raw.mk + (Array.map (fun r ↦ (p.trim).leadingCoeff⁻¹ * r) p.trim)).trim = + CPolynomial.Raw.mk (Array.map (fun r ↦ (p.trim).leadingCoeff⁻¹ * r) p.trim) + apply CPolynomial.Raw.Trim.non_zero_map (fun r ↦ q.leadingCoeff⁻¹ * r) + · intro r hr + apply mul_eq_zero.mp at hr + rcases hr with hinv | hr + · have hlead0 : q.leadingCoeff = 0 := by + exact inv_eq_zero.mp hinv + have hcanon : q.trim = q := by + simpa [q] using CPolynomial.Raw.Trim.trim_twice p + have hcrit := (CPolynomial.Raw.Trim.trim_eq_iff_size_eq_zero_or_getLastD_ne_zero + (p := q)).mp hcanon + rcases hcrit with hsize | hlast + · have hqempty : q = (#[] : CPolynomial.Raw F) := by + apply Array.eq_empty_of_size_eq_zero + exact hsize + exact False.elim (hzero (by simp [hqempty])) + · unfold CPolynomial.Raw.leadingCoeff at hlead0 + rw [hcanon] at hlead0 + exact (hlast hlead0).elim + · exact hr + · simpa [q] using CPolynomial.Raw.Trim.trim_twice p + +private theorem raw_xPowSubXModWith_trim {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (q : Nat) (modulus : CPolynomial.Raw F) : + (CPolynomial.Raw.xPowSubXModWith M D q modulus).trim = + CPolynomial.Raw.xPowSubXModWith M D q modulus := by + unfold CPolynomial.Raw.xPowSubXModWith + exact CPolynomial.Raw.Trim.trim_twice _ + +/-- The finite-field root product agrees with the normalized Mathlib gcd of the +monic input and its represented Frobenius witness. -/ +theorem finiteFieldRootProductWith_toPoly_eq_normalize_gcd {F : Type*} + [Field F] [BEq F] [LawfulBEq F] [DecidableEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : FiniteFieldContext F) {p : CPolynomial F} + (hp : p ≠ 0) : + let pMonic := CPolynomial.monicNormalize p + (finiteFieldRootProductWith M D ctx p).toPoly = + normalize (EuclideanDomain.gcd pMonic.toPoly + (CPolynomial.ofArray + (CPolynomial.Raw.xPowSubXModWith M D ctx.q pMonic.val)).toPoly) := by + dsimp + unfold finiteFieldRootProductWith + unfold CPolynomial.Raw.Roots.FiniteField.finiteFieldRootProductWith + have hpraw : p.val.trim ≠ (0 : CPolynomial.Raw F) := by + intro hpraw0 + apply hp + apply CPolynomial.ext + rw [CPolynomial.trim_eq] at hpraw0 + simpa using hpraw0 + have hpzero : ¬(p.val.trim == (0 : CPolynomial.Raw F)) := by + intro hzero + exact hpraw (LawfulBEq.eq_of_beq hzero) + rw [if_neg hpzero] + have hpMonicVal : + (CPolynomial.monicNormalize p).val = CPolynomial.Raw.monicNormalize p.val := by + unfold CPolynomial.monicNormalize CPolynomial.ofArray + change (CPolynomial.Raw.monicNormalize p.val).trim = + CPolynomial.Raw.monicNormalize p.val + exact raw_monicNormalize_trim p.val + rw [← hpMonicVal] + let witnessRaw := + CPolynomial.Raw.xPowSubXModWith M D ctx.q (CPolynomial.monicNormalize p).val + have hwitnessTrim : witnessRaw.trim = witnessRaw := by + dsimp [witnessRaw] + exact raw_xPowSubXModWith_trim M D ctx.q (CPolynomial.monicNormalize p).val + simpa [CPolynomial.gcdMonic, CPolynomial.ofArray, witnessRaw, hwitnessTrim] using + CPolynomial.gcdMonic_toPoly_eq_normalize_gcd + (CPolynomial.monicNormalize p) (CPolynomial.ofArray witnessRaw) + +private theorem raw_monicNormalize_trim_arg {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (p : CPolynomial.Raw F) : + CPolynomial.Raw.monicNormalize p.trim = CPolynomial.Raw.monicNormalize p := by + unfold CPolynomial.Raw.monicNormalize + rw [CPolynomial.Raw.Trim.trim_twice] + +private theorem raw_monicNormalize_toPoly_eq_normalize {F : Type*} + [Field F] [BEq F] [LawfulBEq F] [DecidableEq F] + (p : CPolynomial.Raw F) : + (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize p)).toPoly = + normalize (CPolynomial.ofArray p).toPoly := by + have h := CPolynomial.monicNormalize_toPoly_eq_normalize (CPolynomial.ofArray p) + unfold CPolynomial.monicNormalize at h + unfold CPolynomial.ofArray at h + rw [raw_monicNormalize_trim_arg] at h + exact h + +private theorem raw_monicNormalize_toPoly_monic {F : Type*} + [Field F] [BEq F] [LawfulBEq F] [DecidableEq F] + {p : CPolynomial.Raw F} (hp : p.trim ≠ 0) : + (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize p)).toPoly.Monic := by + rw [raw_monicNormalize_toPoly_eq_normalize] + have hpC : CPolynomial.ofArray p ≠ 0 := by + intro h + have hval := congrArg Subtype.val h + unfold CPolynomial.ofArray at hval + exact hp (by simpa using hval) + have hpPoly : (CPolynomial.ofArray p).toPoly ≠ 0 := + (CPolynomial.toPoly_eq_zero_iff (CPolynomial.ofArray p)).not.mpr hpC + exact Polynomial.monic_normalize hpPoly + +private theorem raw_modContext_toPoly_eq_modByMonic {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (D : CPolynomial.Raw.ModContext F) {p q : CPolynomial.Raw F} + (hpTrim : p.trim = p) (hqTrim : q.trim = q) + (hqMonic : (CPolynomial.ofArray q).monic) : + (CPolynomial.ofArray (D.modByMonic p q)).toPoly = + (CPolynomial.ofArray p).toPoly %ₘ (CPolynomial.ofArray q).toPoly := by + rw [D.modByMonic_eq_modByMonic p q hpTrim hqTrim] + have hpval : (CPolynomial.ofArray p).val = p := by + unfold CPolynomial.ofArray + exact hpTrim + have hqval : (CPolynomial.ofArray q).val = q := by + unfold CPolynomial.ofArray + exact hqTrim + have h := CPolynomial.modByMonic_toPoly_eq_modByMonic + (CPolynomial.ofArray p) (CPolynomial.ofArray q) hqMonic + simpa [CPolynomial.modByMonic, hpval, hqval, CPolynomial.ofArray_toPoly] using h + +private theorem raw_mulModWith_toPoly_eq_modByMonic {F : Type*} + [Field F] [BEq F] [LawfulBEq F] [DecidableEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + {modulus p q : CPolynomial.Raw F} (hmod : modulus.trim ≠ 0) : + (CPolynomial.ofArray (CPolynomial.Raw.mulModWith M D modulus p q)).toPoly = + ((CPolynomial.ofArray p).toPoly * (CPolynomial.ofArray q).toPoly) %ₘ + (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).toPoly := by + unfold CPolynomial.Raw.mulModWith + have hzero : ¬modulus.trim == (0 : CPolynomial.Raw F) := by + intro h + exact hmod (LawfulBEq.eq_of_beq h) + rw [if_neg hzero] + have hproductTrim : (M.mul p q).trim = M.mul p q := by + rw [M.mul_eq_mul] + exact CPolynomial.Raw.mul_is_trimmed p q + have hqMonic : (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).monic := + (CPolynomial.monic_toPoly_iff _).mpr (raw_monicNormalize_toPoly_monic hmod) + rw [raw_modContext_toPoly_eq_modByMonic D hproductTrim + (raw_monicNormalize_trim modulus) hqMonic] + rw [M.mul_eq_mul, CPolynomial.ofArray_toPoly, CPolynomial.ofArray_toPoly, + CPolynomial.Raw.toPoly_mul] + rw [CPolynomial.ofArray_toPoly, CPolynomial.ofArray_toPoly] + +private theorem polynomial_modByMonic_idem {F : Type*} [Field F] + {p q : Polynomial F} (hq : q.Monic) : + (p %ₘ q) %ₘ q = p %ₘ q := by + refine Polynomial.modByMonic_eq_of_dvd_sub hq ?_ + refine ⟨-(p /ₘ q), ?_⟩ + calc + p %ₘ q - p = p %ₘ q - (p %ₘ q + q * (p /ₘ q)) := by + rw [Polynomial.modByMonic_add_div p] + _ = q * -(p /ₘ q) := by + ring + +private theorem polynomial_mul_modByMonic_congr {F : Type*} [Field F] + {a b c d m : Polynomial F} + (ha : a %ₘ m = b %ₘ m) (hc : c %ₘ m = d %ₘ m) : + (a * c) %ₘ m = (b * d) %ₘ m := by + calc + (a * c) %ₘ m = (a %ₘ m * (c %ₘ m)) %ₘ m := by + rw [Polynomial.mul_modByMonic] + _ = (b %ₘ m * (d %ₘ m)) %ₘ m := by + rw [ha, hc] + _ = (b * d) %ₘ m := by + exact (Polynomial.mul_modByMonic b d m).symm + +private theorem polynomial_pow_modByMonic_congr {F : Type*} [Field F] + {a b m : Polynomial F} (h : a %ₘ m = b %ₘ m) : + ∀ n, a ^ n %ₘ m = b ^ n %ₘ m := by + intro n + induction n with + | zero => simp + | succ n ih => + rw [pow_succ, pow_succ] + exact polynomial_mul_modByMonic_congr ih h + +private theorem raw_powModBinaryAuxWith_toPoly_modByMonic {F : Type*} + [Field F] [BEq F] [LawfulBEq F] [DecidableEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + {modulus acc current : CPolynomial.Raw F} (hmod : modulus.trim ≠ 0) : + ∀ n, + (CPolynomial.ofArray + (CPolynomial.Raw.powModBinaryAuxWith M D modulus n acc current)).toPoly %ₘ + (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).toPoly = + ((CPolynomial.ofArray acc).toPoly * (CPolynomial.ofArray current).toPoly ^ n) %ₘ + (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).toPoly := by + intro n + induction n using Nat.strongRecOn generalizing acc current with + | ind n ih => + cases n with + | zero => + rw [CPolynomial.Raw.powModBinaryAuxWith, pow_zero] + congr 1 + ring + | succ n => + rw [CPolynomial.Raw.powModBinaryAuxWith] + let m := (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).toPoly + let acc' := + if (n + 1) % 2 == 1 then + CPolynomial.Raw.mulModWith M D modulus acc current + else + acc + let current' := CPolynomial.Raw.mulModWith M D modulus current current + have ih' := ih ((n + 1) / 2) + (Nat.div_lt_self (Nat.succ_pos n) (by decide)) + (acc := acc') (current := current') + change + (CPolynomial.ofArray + (CPolynomial.Raw.powModBinaryAuxWith M D modulus + ((n + 1) / 2) acc' current')).toPoly %ₘ + m = + ((CPolynomial.ofArray acc).toPoly * + (CPolynomial.ofArray current).toPoly ^ (n + 1)) %ₘ m + rw [ih'] + have hcurrent : + (CPolynomial.ofArray current').toPoly %ₘ m = + ((CPolynomial.ofArray current).toPoly * + (CPolynomial.ofArray current).toPoly) %ₘ m := by + dsimp [current', m] + rw [raw_mulModWith_toPoly_eq_modByMonic M D hmod] + exact polynomial_modByMonic_idem (raw_monicNormalize_toPoly_monic hmod) + by_cases hodd : (n + 1) % 2 == 1 + · have hacc : + (CPolynomial.ofArray acc').toPoly %ₘ m = + ((CPolynomial.ofArray acc).toPoly * + (CPolynomial.ofArray current).toPoly) %ₘ m := by + dsimp [acc', m] + rw [if_pos hodd] + rw [raw_mulModWith_toPoly_eq_modByMonic M D hmod] + exact polynomial_modByMonic_idem (raw_monicNormalize_toPoly_monic hmod) + have hpow : + ((CPolynomial.ofArray current).toPoly * + (CPolynomial.ofArray current).toPoly) ^ ((n + 1) / 2) = + (CPolynomial.ofArray current).toPoly ^ n := by + rw [mul_pow, ← pow_add] + congr 1 + have hoddNat : (n + 1) % 2 = 1 := by + simpa using hodd + omega + have hcurrentPow : + (CPolynomial.ofArray current').toPoly ^ ((n + 1) / 2) %ₘ m = + (CPolynomial.ofArray current).toPoly ^ n %ₘ m := by + have h := polynomial_pow_modByMonic_congr hcurrent ((n + 1) / 2) + rwa [hpow] at h + calc + ((CPolynomial.ofArray acc').toPoly * + (CPolynomial.ofArray current').toPoly ^ ((n + 1) / 2)) %ₘ m = + (((CPolynomial.ofArray acc).toPoly * + (CPolynomial.ofArray current).toPoly) * + (CPolynomial.ofArray current).toPoly ^ n) %ₘ m := by + exact polynomial_mul_modByMonic_congr hacc hcurrentPow + _ = ((CPolynomial.ofArray acc).toPoly * + (CPolynomial.ofArray current).toPoly ^ (n + 1)) %ₘ m := by + rw [pow_succ] + ring_nf + · have hacc : + (CPolynomial.ofArray acc').toPoly %ₘ m = + (CPolynomial.ofArray acc).toPoly %ₘ m := by + dsimp [acc'] + rw [if_neg hodd] + have hpow : + ((CPolynomial.ofArray current).toPoly * + (CPolynomial.ofArray current).toPoly) ^ ((n + 1) / 2) = + (CPolynomial.ofArray current).toPoly ^ (n + 1) := by + rw [mul_pow, ← pow_add] + congr 1 + have hoddNat : (n + 1) % 2 = 0 := by + simpa using hodd + omega + have hcurrentPow : + (CPolynomial.ofArray current').toPoly ^ ((n + 1) / 2) %ₘ m = + (CPolynomial.ofArray current).toPoly ^ (n + 1) %ₘ m := by + have h := polynomial_pow_modByMonic_congr hcurrent ((n + 1) / 2) + rwa [hpow] at h + exact polynomial_mul_modByMonic_congr hacc hcurrentPow + +private theorem raw_powModWith_X_toPoly_modByMonic {F : Type*} + [Field F] [BEq F] [LawfulBEq F] [DecidableEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + {modulus : CPolynomial.Raw F} (hmod : modulus.trim ≠ 0) (q : Nat) : + (CPolynomial.ofArray + (CPolynomial.Raw.powModWith M D modulus CPolynomial.Raw.X q)).toPoly %ₘ + (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).toPoly = + ((Polynomial.X : Polynomial F) ^ q) %ₘ + (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).toPoly := by + unfold CPolynomial.Raw.powModWith + have hzero : ¬modulus.trim == (0 : CPolynomial.Raw F) := by + intro h + exact hmod (LawfulBEq.eq_of_beq h) + rw [if_neg hzero] + let m := (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).toPoly + let oneMod := D.modByMonic (1 : CPolynomial.Raw F) (CPolynomial.Raw.monicNormalize modulus) + change + (CPolynomial.ofArray + (CPolynomial.Raw.powModBinaryAuxWith M D modulus q oneMod CPolynomial.Raw.X)).toPoly %ₘ + m = + ((Polynomial.X : Polynomial F) ^ q) %ₘ m + rw [raw_powModBinaryAuxWith_toPoly_modByMonic M D hmod] + have hOneTrim : (1 : CPolynomial.Raw F).trim = 1 := by + change CPolynomial.Raw.trim (#[] |>.push (1 : F)) = (#[] |>.push (1 : F)) + apply CPolynomial.Raw.Trim.push_trim + simp + have hone : + (CPolynomial.ofArray oneMod).toPoly %ₘ m = + (1 : Polynomial F) %ₘ m := by + dsimp [oneMod, m] + rw [raw_modContext_toPoly_eq_modByMonic D hOneTrim + (raw_monicNormalize_trim modulus) + ((CPolynomial.monic_toPoly_iff _).mpr (raw_monicNormalize_toPoly_monic hmod))] + rw [CPolynomial.ofArray_toPoly, CPolynomial.Raw.toPoly_one] + exact polynomial_modByMonic_idem (raw_monicNormalize_toPoly_monic hmod) + have hX : + (CPolynomial.ofArray CPolynomial.Raw.X).toPoly ^ q %ₘ m = + ((Polynomial.X : Polynomial F) ^ q) %ₘ m := by + rw [CPolynomial.ofArray_toPoly, CPolynomial.Raw.toPoly_X] + have hmul := polynomial_mul_modByMonic_congr hone hX + simpa [one_mul] using hmul + +private theorem raw_xModWith_toPoly_modByMonic {F : Type*} + [Field F] [BEq F] [LawfulBEq F] [DecidableEq F] + (D : CPolynomial.Raw.ModContext F) + {modulus : CPolynomial.Raw F} (hmod : modulus.trim ≠ 0) : + (CPolynomial.ofArray (CPolynomial.Raw.xModWith D modulus)).toPoly %ₘ + (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).toPoly = + (Polynomial.X : Polynomial F) %ₘ + (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).toPoly := by + unfold CPolynomial.Raw.xModWith + have hzero : ¬modulus.trim == (0 : CPolynomial.Raw F) := by + intro h + exact hmod (LawfulBEq.eq_of_beq h) + rw [if_neg hzero] + let m := (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).toPoly + have hXTrim : (CPolynomial.Raw.X : CPolynomial.Raw F).trim = CPolynomial.Raw.X := by + exact CPolynomial.Raw.X_canonical + rw [raw_modContext_toPoly_eq_modByMonic D hXTrim + (raw_monicNormalize_trim modulus) + ((CPolynomial.monic_toPoly_iff _).mpr (raw_monicNormalize_toPoly_monic hmod))] + rw [CPolynomial.ofArray_toPoly, CPolynomial.Raw.toPoly_X] + exact polynomial_modByMonic_idem (raw_monicNormalize_toPoly_monic hmod) + +private theorem raw_xPowSubXModWith_toPoly_modByMonic {F : Type*} + [Field F] [BEq F] [LawfulBEq F] [DecidableEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + {modulus : CPolynomial.Raw F} (hmod : modulus.trim ≠ 0) (q : Nat) : + (CPolynomial.ofArray + (CPolynomial.Raw.xPowSubXModWith M D q modulus)).toPoly %ₘ + (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).toPoly = + (((Polynomial.X : Polynomial F) ^ q - Polynomial.X) %ₘ + (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize modulus)).toPoly) := by + unfold CPolynomial.Raw.xPowSubXModWith + rw [CPolynomial.ofArray_toPoly, CPolynomial.Raw.toPoly_sub, Polynomial.sub_modByMonic] + rw [← CPolynomial.ofArray_toPoly + (CPolynomial.Raw.powModWith M D modulus CPolynomial.Raw.X q), + ← CPolynomial.ofArray_toPoly (CPolynomial.Raw.xModWith D modulus)] + rw [raw_powModWith_X_toPoly_modByMonic M D hmod q, + raw_xModWith_toPoly_modByMonic D hmod] + rw [Polynomial.sub_modByMonic] + +/-- The finite-field root product divides the finite-field Frobenius polynomial. + +This is the modular-gcd bridge: the executable root product uses the congruent +modular witness `(X^q mod p) - (X mod p)` instead of materializing `X^q - X`. +-/ +theorem finiteFieldRootProductWith_dvd_frobenius_of_context {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : FiniteFieldContext F) {p : CPolynomial F} (_hp : p ≠ 0) : + (finiteFieldRootProductWith M D ctx p).toPoly ∣ + ((Polynomial.X : Polynomial F) ^ ctx.q - Polynomial.X) := by + letI : DecidableEq F := instDecidableEqOfLawfulBEq + rw [finiteFieldRootProductWith_toPoly_eq_normalize_gcd M D ctx _hp] + apply dvd_trans (normalize_associated _).dvd + let pMonic := CPolynomial.monicNormalize p + let witness := + (CPolynomial.ofArray + (CPolynomial.Raw.xPowSubXModWith M D ctx.q pMonic.val)).toPoly + have hgLeft : EuclideanDomain.gcd pMonic.toPoly witness ∣ pMonic.toPoly := + EuclideanDomain.gcd_dvd_left _ _ + have hgRight : EuclideanDomain.gcd pMonic.toPoly witness ∣ witness := + EuclideanDomain.gcd_dvd_right _ _ + have hpPoly : p.toPoly ≠ 0 := + (CPolynomial.toPoly_eq_zero_iff p).not.mpr _hp + have hpMonicPoly : pMonic.toPoly.Monic := by + dsimp [pMonic] + rw [CPolynomial.monicNormalize_toPoly_eq_normalize] + exact Polynomial.monic_normalize hpPoly + have hpMonicNe : pMonic ≠ 0 := by + intro hzero + dsimp [pMonic] at hzero + unfold CPolynomial.monicNormalize at hzero + have hval := congrArg Subtype.val hzero + simp [CPolynomial.ofArray] at hval + have hpraw : p.val.trim ≠ (0 : CPolynomial.Raw F) := by + intro hpraw0 + apply _hp + apply CPolynomial.ext + rw [CPolynomial.trim_eq] at hpraw0 + simpa using hpraw0 + rw [raw_monicNormalize_trim] at hval + exact raw_monicNormalize_ne_zero_of_trim_ne_zero hpraw hval + have hmod : pMonic.val.trim ≠ 0 := by + intro htrim + apply hpMonicNe + apply CPolynomial.ext + rw [CPolynomial.trim_eq] at htrim + simpa using htrim + have hnormModulus : + (CPolynomial.ofArray (CPolynomial.Raw.monicNormalize pMonic.val)).toPoly = + pMonic.toPoly := by + rw [raw_monicNormalize_toPoly_eq_normalize] + rw [CPolynomial.ofArray_toPoly] + exact hpMonicPoly.normalize_eq_self + have hwitnessMod : + witness %ₘ pMonic.toPoly = + (((Polynomial.X : Polynomial F) ^ ctx.q - Polynomial.X) %ₘ pMonic.toPoly) := by + dsimp [witness] + have h := raw_xPowSubXModWith_toPoly_modByMonic M D hmod ctx.q + rwa [hnormModulus] at h + have hpMonic : + pMonic.toPoly ∣ (Polynomial.X : Polynomial F) ^ ctx.q - Polynomial.X - witness := by + rw [← Polynomial.modByMonic_eq_zero_iff_dvd hpMonicPoly] + rw [Polynomial.sub_modByMonic] + rw [← hwitnessMod] + simp + simpa [sub_eq_add_neg, add_assoc, add_left_comm, add_comm] using + dvd_add (hgLeft.trans hpMonic) hgRight + +private theorem raw_gcdMonicWithFuel_trim_ne_zero_of_left {F : Type*} + [Field F] [BEq F] [LawfulBEq F] : + ∀ fuel (p q : CPolynomial.Raw F), p.trim ≠ (0 : CPolynomial.Raw F) → + (CPolynomial.Raw.gcdMonicWithFuel fuel p q).trim ≠ 0 := by + intro fuel + induction fuel with + | zero => + intro p q hp + unfold CPolynomial.Raw.gcdMonicWithFuel + rw [raw_monicNormalize_trim] + exact raw_monicNormalize_ne_zero_of_trim_ne_zero hp + | succ fuel ih => + intro p q hp + rw [CPolynomial.Raw.gcdMonicWithFuel] + by_cases hqzero : q.trim == (0 : CPolynomial.Raw F) + · rw [if_pos hqzero] + rw [raw_monicNormalize_trim] + apply raw_monicNormalize_ne_zero_of_trim_ne_zero + simpa [CPolynomial.Raw.Trim.trim_twice] using hp + · have hqzeroNe : ¬q.trim = (0 : CPolynomial.Raw F) := by + intro hqtrim0 + exact hqzero (by simp [hqtrim0]) + rw [if_neg hqzero] + apply ih + intro hqtrim0 + rw [CPolynomial.Raw.Trim.trim_twice q] at hqtrim0 + exact hqzeroNe hqtrim0 + +/-- Monic normalization preserves nonzeroness. -/ +theorem monicNormalize_ne_zero_of_ne_zero {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p : CPolynomial F} (hp : p ≠ 0) : + CPolynomial.monicNormalize p ≠ 0 := by + intro hzero + unfold CPolynomial.monicNormalize at hzero + have hval := congrArg Subtype.val hzero + simp [CPolynomial.ofArray] at hval + have hpraw : p.val.trim ≠ (0 : CPolynomial.Raw F) := by + intro hpraw0 + apply hp + apply CPolynomial.ext + rw [CPolynomial.trim_eq] at hpraw0 + simpa using hpraw0 + rw [raw_monicNormalize_trim] at hval + exact raw_monicNormalize_ne_zero_of_trim_ne_zero hpraw hval + +/-- The monic gcd of a nonzero left operand is nonzero. -/ +theorem gcdMonic_ne_zero_of_left {F : Type*} [Field F] [BEq F] [LawfulBEq F] + {p q : CPolynomial F} (hp : p ≠ 0) : + CPolynomial.gcdMonic p q ≠ 0 := by + intro hzero + unfold CPolynomial.gcdMonic at hzero + have hval := congrArg Subtype.val hzero + simp [CPolynomial.ofArray] at hval + have hpraw : p.val.trim ≠ (0 : CPolynomial.Raw F) := by + intro hpraw0 + apply hp + apply CPolynomial.ext + rw [CPolynomial.trim_eq] at hpraw0 + simpa using hpraw0 + unfold CPolynomial.Raw.gcdMonic at hval + exact raw_gcdMonicWithFuel_trim_ne_zero_of_left _ p.val q.val hpraw hval + +/-- The finite-field root product of a nonzero polynomial is nonzero. -/ +theorem finiteFieldRootProductWith_ne_zero_of_ne_zero {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : FiniteFieldContext F) {p : CPolynomial F} + (hp : p ≠ 0) : + finiteFieldRootProductWith M D ctx p ≠ 0 := by + intro hzeroProduct + unfold finiteFieldRootProductWith at hzeroProduct + have hval := congrArg Subtype.val hzeroProduct + simp [CPolynomial.ofArray] at hval + have hpraw : p.val.trim ≠ (0 : CPolynomial.Raw F) := by + intro hpraw0 + apply hp + apply CPolynomial.ext + rw [CPolynomial.trim_eq] at hpraw0 + simpa using hpraw0 + unfold CPolynomial.Raw.Roots.FiniteField.finiteFieldRootProductWith at hval + by_cases hpempty : p.val = (#[] : CPolynomial.Raw F) + · apply hp + apply CPolynomial.ext + simpa using hpempty + · simp [hpempty, CPolynomial.trim_eq] at hval + unfold CPolynomial.Raw.gcdMonic at hval + exact raw_gcdMonicWithFuel_trim_ne_zero_of_left + _ (CPolynomial.Raw.monicNormalize p.val) + (CPolynomial.Raw.xPowSubXModWith M D ctx.q (CPolynomial.Raw.monicNormalize p.val)) + (by + rw [raw_monicNormalize_trim] + exact raw_monicNormalize_ne_zero_of_trim_ne_zero hpraw) hval + +end FiniteField + +end Roots + +end CPolynomial + +end CompPoly diff --git a/CompPoly/Univariate/Roots/SmoothSubgroup.lean b/CompPoly/Univariate/Roots/SmoothSubgroup.lean new file mode 100644 index 00000000..8e88e42a --- /dev/null +++ b/CompPoly/Univariate/Roots/SmoothSubgroup.lean @@ -0,0 +1,14 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Roots.SmoothSubgroup.Basic +import CompPoly.Univariate.Roots.SmoothSubgroup.Correctness + +/-! +# Smooth Multiplicative-Subgroup Root Splitting + +Public module for the smooth cyclic subgroup splitter. +-/ diff --git a/CompPoly/Univariate/Roots/SmoothSubgroup/Basic.lean b/CompPoly/Univariate/Roots/SmoothSubgroup/Basic.lean new file mode 100644 index 00000000..8b4bb4cb --- /dev/null +++ b/CompPoly/Univariate/Roots/SmoothSubgroup/Basic.lean @@ -0,0 +1,225 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.BatchEval.Context +import CompPoly.Univariate.Roots.Splitter +import Mathlib.GroupTheory.OrderOfElement + +/-! +# Smooth Multiplicative-Subgroup Linear-Factor Splitting + +Executable splitter for field-root products over finite fields whose +multiplicative group admits a smooth cyclic refinement schedule, following the +subgroup-refinement root-finding method of [MOV92]. + +The generic root pipeline consumes a `LinearFactorProductSplitter`. This module +provides a contract-bearing smooth context plus an adapter from that context to +the shared splitter interface. + +## References + +* [Menezes, A. J., van Oorschot, P. C., and Vanstone, S. A., *Subgroup + Refinement Algorithms for Root Finding in GF(q)*][MOV92] +-/ + +namespace CompPoly + +namespace CPolynomial + +namespace Roots + +namespace FiniteField + +/-- +Splitter-input predicate for a smooth cyclic splitter. + +This predicate records the mathematical input contract for the splitter. The +field-root pipeline proves it for `gcd(p, X^q - X)` values before using the +splitter completeness theorem. +-/ +def smoothSplitterInput {F : Type*} [Field F] + (q : Nat) (_generator : F) (_schedule : Array Nat) (p : CPolynomial F) : Prop := + p ≠ 0 ∧ ∀ a : F, CPolynomial.eval a p = 0 → a ^ q = a + +/-- Package the executable smooth splitter with its field and schedule facts. -/ +structure SmoothCyclicRootContext (F : Type*) [Field F] [BEq F] [LawfulBEq F] where + q : Nat + generator : F + schedule : Array Nat + leafEvaluator : BatchEvalContext F + validInput : CPolynomial F → Prop + card_eq : Nat.card F = q + generator_order : orderOf generator = q - 1 + schedule_complete : + schedule.toList.foldl (fun order ell ↦ order / ell) (q - 1) = 1 + splitLinearFactorsWith : + CPolynomial.Raw.MulContext F → + CPolynomial.Raw.ModContext F → + CPolynomial F → Array (CPolynomial F) + sound : + ∀ M D p factor, + factor ∈ (splitLinearFactorsWith M D p).toList → + IsLinearFactor factor + complete : + ∀ M D p a, + validInput p → + p ≠ 0 → + CPolynomial.eval a p = 0 → + ∃ factor, + factor ∈ (splitLinearFactorsWith M D p).toList ∧ + IsLinearRootFactorCandidate factor a + +/-- Ordered elements of the coset `alpha * ` with the supplied order. -/ +def smoothCosetPoints {F : Type*} [Field F] + (alpha gamma : F) (order : Nat) : Array F := + Array.ofFn fun j : Fin order ↦ alpha * gamma ^ j.val + +/-- Emit linear factors for the points whose evaluated value is zero. -/ +def linearFactorsFromLeafValues {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (points values : Array F) : Array (CPolynomial F) := + points.zipIdx.foldl + (fun factors ⟨point, idx⟩ ↦ + if values.getD idx 1 == 0 then + factors.push (CPolynomial.linearFactor point) + else + factors) + #[] + +/-- Leaf extraction by evaluating all elements of the current coset. -/ +def smoothLeafLinearFactors {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (E : BatchEvalContext F) (alpha gamma : F) (order : Nat) + (p : CPolynomial F) : Array (CPolynomial F) := + let points := smoothCosetPoints alpha gamma order + linearFactorsFromLeafValues points (E.evalBatchWith p points) + +/-- Schedule-driven nonzero-root refinement inside one multiplicative coset. -/ +def smoothCosetLinearFactorsWithSchedule {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (E : BatchEvalContext F) : + List Nat → Nat → F → F → CPolynomial F → Array (CPolynomial F) + | [], order, alpha, gamma, p => + let p := CPolynomial.monicNormalize p + if p == 0 || p == 1 then + #[] + else if isRepresentedLinearFactor p then + #[p] + else + smoothLeafLinearFactors E alpha gamma order p + | ell :: rest, order, alpha, gamma, p => + let p := CPolynomial.monicNormalize p + if p == 0 || p == 1 then + #[] + else if isRepresentedLinearFactor p then + #[p] + else if ell ≤ 1 then + smoothCosetLinearFactorsWithSchedule M D E rest order alpha gamma p + else if ell = order then + smoothLeafLinearFactors E alpha gamma order p + else + let childOrder := order / ell + let tau := gamma ^ childOrder + let xPow := xPowModWith M D p childOrder + (List.range ell).foldl + (fun factors j ↦ + let beta := alpha ^ childOrder * tau ^ j + let witness := xPow - CPolynomial.C beta + let child := CPolynomial.monicNormalize (CPolynomial.gcdMonic p witness) + if child == 0 || child == 1 then + factors + else + factors ++ + smoothCosetLinearFactorsWithSchedule M D E rest childOrder + (alpha * gamma ^ j) (gamma ^ ell) child) + #[] + +/-- Nonzero root extraction by smooth multiplicative-subgroup refinement. -/ +def smoothNonzeroLinearFactorsWith {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (E : BatchEvalContext F) (q : Nat) (generator : F) (schedule : Array Nat) + (p : CPolynomial F) : Array (CPolynomial F) := + smoothCosetLinearFactorsWithSchedule M D E schedule.toList (q - 1) + (1 : F) generator p + +/-- Smooth linear-factor splitting algorithm, including separate handling of root `0`. -/ +def smoothLinearFactorsAlgorithmWith {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (E : BatchEvalContext F) (q : Nat) (generator : F) (schedule : Array Nat) + (p : CPolynomial F) : Array (CPolynomial F) := + let p := CPolynomial.monicNormalize p + if p == 0 || p == 1 then + #[] + else + let zeroSplit : Array (CPolynomial F) × CPolynomial F := + if p.coeff 0 == 0 then + (#[CPolynomial.linearFactor (0 : F)], CPolynomial.monicNormalize (CPolynomial.divX p)) + else + (#[], p) + zeroSplit.1 ++ + smoothNonzeroLinearFactorsWith M D E q generator schedule zeroSplit.2 + +/-- Build a smooth cyclic root context from executable constants and proof fields. -/ +def smoothCyclicRootContextOf {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (q : Nat) (generator : F) (schedule : Array Nat) + (leafEvaluator : BatchEvalContext F) + (validInput : CPolynomial F → Prop) + (card_eq : Nat.card F = q) + (generator_order : orderOf generator = q - 1) + (schedule_complete : + schedule.toList.foldl (fun order ell ↦ order / ell) (q - 1) = 1) + (sound : + ∀ M D p factor, + factor ∈ + (smoothLinearFactorsAlgorithmWith M D leafEvaluator q generator schedule p).toList → + IsLinearFactor factor) + (complete : + ∀ M D p a, + validInput p → + p ≠ 0 → + CPolynomial.eval a p = 0 → + ∃ factor, + factor ∈ + (smoothLinearFactorsAlgorithmWith M D leafEvaluator q generator schedule + p).toList ∧ + IsLinearRootFactorCandidate factor a) : + SmoothCyclicRootContext F where + q := q + generator := generator + schedule := schedule + leafEvaluator := leafEvaluator + validInput := validInput + card_eq := card_eq + generator_order := generator_order + schedule_complete := schedule_complete + splitLinearFactorsWith := fun M D ↦ + smoothLinearFactorsAlgorithmWith M D leafEvaluator q generator schedule + sound := sound + complete := complete + +/-- Adapt a smooth cyclic context to the generic splitter interface. -/ +def smoothLinearFactorProductSplitterWith {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : SmoothCyclicRootContext F) : + LinearFactorProductSplitter F where + splitLinearFactors := fun _qArg p ↦ ctx.splitLinearFactorsWith M D p + validInput := fun _qArg p ↦ ctx.validInput p + sound := by + intro _qArg p factor h + exact ctx.sound M D p factor h + complete := by + intro _qArg p a hvalid hp hroot + exact ctx.complete M D p a hvalid hp hroot + +end FiniteField + +end Roots + +end CPolynomial + +end CompPoly diff --git a/CompPoly/Univariate/Roots/SmoothSubgroup/Correctness.lean b/CompPoly/Univariate/Roots/SmoothSubgroup/Correctness.lean new file mode 100644 index 00000000..33834bed --- /dev/null +++ b/CompPoly/Univariate/Roots/SmoothSubgroup/Correctness.lean @@ -0,0 +1,774 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Roots.Correctness +import CompPoly.Univariate.Roots.RootProduct +import CompPoly.Univariate.Roots.SmoothSubgroup.Basic +import Mathlib.Algebra.Group.Subgroup.Finite + +/-! +# Smooth Multiplicative-Subgroup Splitter Correctness + +Correctness theorems for the smooth cyclic splitter, including executable +contracts and adapter theorems. +-/ + +namespace CompPoly + +namespace CPolynomial + +namespace Roots + +namespace FiniteField + +/-- Field-root products satisfy the generic smooth-splitter input predicate. -/ +theorem finiteFieldRootProductWith_smoothSplitterInput {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : FiniteFieldContext F) (generator : F) (schedule : Array Nat) + {p : CPolynomial F} (hp : p ≠ 0) : + smoothSplitterInput ctx.q generator schedule + (finiteFieldRootProductWith M D ctx p) := by + constructor + · exact CPolynomial.Roots.FiniteField.finiteFieldRootProductWith_ne_zero_of_ne_zero + M D ctx hp + · intro a _ha + exact ctx.frobenius_fixed a + +/-- Soundness theorem for a smooth context adapted to the splitter interface. -/ +theorem smoothLinearFactorProductSplitterWith_sound {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : SmoothCyclicRootContext F) {q : Nat} {p factor : CPolynomial F} + (h : factor ∈ + ((smoothLinearFactorProductSplitterWith M D ctx).splitLinearFactors q p).toList) : + IsLinearFactor factor := by + exact (smoothLinearFactorProductSplitterWith M D ctx).sound q p factor h + +/-- Completeness theorem for a smooth context adapted to the splitter interface. -/ +theorem smoothLinearFactorProductSplitterWith_complete {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (ctx : SmoothCyclicRootContext F) {q : Nat} {p : CPolynomial F} {a : F} + (hvalid : (smoothLinearFactorProductSplitterWith M D ctx).validInput q p) + (hp : p ≠ 0) (hroot : CPolynomial.eval a p = 0) : + ∃ factor, + factor ∈ + ((smoothLinearFactorProductSplitterWith M D ctx).splitLinearFactors q p).toList ∧ + IsLinearRootFactorCandidate factor a := by + exact (smoothLinearFactorProductSplitterWith M D ctx).complete q p a hvalid hp hroot + +/-- Zero-root extraction is sound for the emitted `X` factor. -/ +theorem smooth_zero_root_factor_sound {F : Type*} [Field F] [BEq F] [LawfulBEq F] : + IsLinearRootFactorCandidate (CPolynomial.linearFactor (0 : F)) (0 : F) := by + have hC0 : CPolynomial.C (0 : F) = 0 := by + apply (CPolynomial.eq_iff_coeff).2 + intro i + simp + rw [CPolynomial.linearFactor, neg_zero, hC0, CPolynomial.zero_add] + simp [IsLinearRootFactorCandidate, IsLinearFactor, CPolynomial.X, CPolynomial.Raw.X, + CPolynomial.coeff, CPolynomial.Raw.coeff] + +/-- Leaf extraction emits only represented nonconstant linear factors. -/ +private theorem linearFactorsFromLeafValues_sound {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (points values : Array F) {factor : CPolynomial F} + (h : factor ∈ (linearFactorsFromLeafValues points values).toList) : + IsLinearFactor factor := by + unfold linearFactorsFromLeafValues at h + rw [Array.foldl_zipIdx_eq_foldl_toList_zipIdx] at h + have haux : ∀ (xs : List (F × Nat)) (acc : Array (CPolynomial F)), + (∀ factor, factor ∈ acc.toList → IsLinearFactor factor) → + factor ∈ (xs.foldl + (fun factors x ↦ + if values.getD x.2 1 == 0 then + factors.push (CPolynomial.linearFactor x.1) + else + factors) acc).toList → + IsLinearFactor factor := by + intro xs + induction xs with + | nil => + intro acc hacc hmem + exact hacc factor hmem + | cons x xs ih => + intro acc hacc hmem + simp only [List.foldl_cons] at hmem + refine ih (if values.getD x.2 1 == 0 then + acc.push (CPolynomial.linearFactor x.1) + else + acc) ?_ hmem + intro fac hfac + by_cases hcond : (values.getD x.2 1 == 0) = true + · rw [if_pos hcond] at hfac + by_cases hlast : fac = CPolynomial.linearFactor x.1 + · rw [hlast] + exact linearFactor_isLinearFactor x.1 + · apply hacc + simpa [hlast] using hfac + · rw [if_neg hcond] at hfac + exact hacc fac hfac + exact haux points.toList.zipIdx #[] (by simp) h + +/-- A point whose batch value is zero contributes its explicit linear factor. -/ +private theorem mem_linearFactorsFromLeafValues_of_get_eq_zero {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (points values : Array F) {idx : Nat} (hidx : idx < points.size) + (hval : values.getD idx 1 = 0) : + CPolynomial.linearFactor points[idx] ∈ + (linearFactorsFromLeafValues points values).toList := by + unfold linearFactorsFromLeafValues + rw [Array.foldl_zipIdx_eq_foldl_toList_zipIdx] + have htarget : (points[idx], idx) ∈ points.toList.zipIdx := by + have hex : ∃ x ∈ points.toList.zipIdx, x = (points[idx], idx) := by + rw [List.exists_mem_zipIdx'] + refine ⟨idx, ?_, ?_⟩ + · simpa using hidx + · simp [Array.getElem_toList] + rcases hex with ⟨x, hxmem, hxeq⟩ + simpa [hxeq] using hxmem + have haux : ∀ (xs : List (F × Nat)) (acc : Array (CPolynomial F)), + CPolynomial.linearFactor points[idx] ∈ acc.toList ∨ (points[idx], idx) ∈ xs → + CPolynomial.linearFactor points[idx] ∈ (xs.foldl + (fun factors x ↦ + if values.getD x.2 1 == 0 then + factors.push (CPolynomial.linearFactor x.1) + else + factors) acc).toList := by + intro xs + induction xs with + | nil => + intro acc h + simpa using h + | cons x xs ih => + intro acc h + simp only [List.foldl_cons] + apply ih + rcases h with hacc | hx + · left + by_cases hcond : (values.getD x.2 1 == 0) = true + · rw [if_pos hcond] + simp [hacc] + · rw [if_neg hcond] + exact hacc + · simp at hx + rcases hx with hx | hx + · subst x + left + have hbeq : (values.getD idx 1 == 0) = true := by + simp [hval] + rw [if_pos hbeq] + simp + · right + exact hx + exact haux points.toList.zipIdx #[] (Or.inr htarget) + +/-- Smooth leaf extraction emits only represented nonconstant linear factors. -/ +private theorem smoothLeafLinearFactors_sound {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (E : BatchEvalContext F) (alpha gamma : F) (order : Nat) (p : CPolynomial F) + {factor : CPolynomial F} + (h : factor ∈ (smoothLeafLinearFactors E alpha gamma order p).toList) : + IsLinearFactor factor := by + unfold smoothLeafLinearFactors at h + exact linearFactorsFromLeafValues_sound _ _ h + +/-- Smooth leaf extraction is complete for roots in the enumerated coset. -/ +theorem smoothLeafLinearFactors_complete {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (E : BatchEvalContext F) {alpha gamma a : F} {order : Nat} {p : CPolynomial F} + (hcoset : ∃ k : Nat, k < order ∧ a = alpha * gamma ^ k) + (hroot : CPolynomial.eval a p = 0) : + ∃ factor, + factor ∈ (smoothLeafLinearFactors E alpha gamma order p).toList ∧ + IsLinearRootFactorCandidate factor a := by + rcases hcoset with ⟨k, hk, ha⟩ + refine ⟨CPolynomial.linearFactor a, ?_, linearFactor_isRootFactorCandidate a⟩ + unfold smoothLeafLinearFactors + let points := smoothCosetPoints alpha gamma order + have hsize : k < points.size := by + simp [points, smoothCosetPoints] + exact hk + have hpoint : points[k] = a := by + simp [points, smoothCosetPoints, ha] + have hval : (E.evalBatchWith p points).getD k 1 = 0 := by + rw [E.correct] + unfold CPolynomial.evalBatch + rw [Array.getD_map_of_lt points (fun x ↦ CPolynomial.eval x p) 1 hsize] + rw [hpoint] + exact hroot + have hmem := mem_linearFactorsFromLeafValues_of_get_eq_zero points (E.evalBatchWith p points) + hsize hval + simpa [hpoint] using hmem + +/-- Schedule-driven smooth coset recursion emits only represented nonconstant linear factors. -/ +theorem smoothCosetLinearFactorsWithSchedule_sound {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (E : BatchEvalContext F) : + ∀ (schedule : List Nat) (order : Nat) (alpha gamma : F) (p factor : CPolynomial F), + factor ∈ + (smoothCosetLinearFactorsWithSchedule M D E schedule order alpha gamma p).toList → + IsLinearFactor factor := by + intro schedule + induction schedule with + | nil => + intro order alpha gamma p factor h + unfold smoothCosetLinearFactorsWithSchedule at h + by_cases hzero : (CPolynomial.monicNormalize p == 0 || + CPolynomial.monicNormalize p == 1) = true + · rw [if_pos hzero] at h + simp at h + · rw [if_neg hzero] at h + by_cases hlin : isRepresentedLinearFactor (CPolynomial.monicNormalize p) = true + · rw [if_pos hlin] at h + simp at h + rcases h with rfl + exact isRepresentedLinearFactor_sound hlin + · rw [if_neg hlin] at h + exact smoothLeafLinearFactors_sound E alpha gamma order (CPolynomial.monicNormalize p) h + | cons ell rest ih => + intro order alpha gamma p factor h + unfold smoothCosetLinearFactorsWithSchedule at h + by_cases hzero : (CPolynomial.monicNormalize p == 0 || + CPolynomial.monicNormalize p == 1) = true + · rw [if_pos hzero] at h + simp at h + · rw [if_neg hzero] at h + by_cases hlin : isRepresentedLinearFactor (CPolynomial.monicNormalize p) = true + · rw [if_pos hlin] at h + simp at h + rcases h with rfl + exact isRepresentedLinearFactor_sound hlin + · rw [if_neg hlin] at h + by_cases hellsmall : ell ≤ 1 + · rw [if_pos hellsmall] at h + exact ih order alpha gamma (CPolynomial.monicNormalize p) factor h + · rw [if_neg hellsmall] at h + by_cases helleq : ell = order + · rw [if_pos helleq] at h + exact smoothLeafLinearFactors_sound E alpha gamma order + (CPolynomial.monicNormalize p) h + · rw [if_neg helleq] at h + let childOrder := order / ell + let tau := gamma ^ childOrder + let xPow := xPowModWith M D (CPolynomial.monicNormalize p) childOrder + have hfold : ∀ (js : List Nat) (acc : Array (CPolynomial F)), + (∀ factor, factor ∈ acc.toList → IsLinearFactor factor) → + factor ∈ (js.foldl + (fun factors j ↦ + let beta := alpha ^ childOrder * tau ^ j + let witness := xPow - CPolynomial.C beta + let child := CPolynomial.monicNormalize + (CPolynomial.gcdMonic (CPolynomial.monicNormalize p) witness) + if child == 0 || child == 1 then + factors + else + factors ++ smoothCosetLinearFactorsWithSchedule M D E rest childOrder + (alpha * gamma ^ j) (gamma ^ ell) child) acc).toList → + IsLinearFactor factor := by + intro js + induction js with + | nil => + intro acc hacc hmem + exact hacc factor hmem + | cons j js ihjs => + intro acc hacc hmem + simp only [List.foldl_cons] at hmem + let beta := alpha ^ childOrder * tau ^ j + let witness := xPow - CPolynomial.C beta + let child := CPolynomial.monicNormalize + (CPolynomial.gcdMonic (CPolynomial.monicNormalize p) witness) + refine ihjs + (if child == 0 || child == 1 then + acc + else + acc ++ smoothCosetLinearFactorsWithSchedule M D E rest childOrder + (alpha * gamma ^ j) (gamma ^ ell) child) ?_ hmem + intro fac hfac + by_cases hskip : (child == 0 || child == 1) = true + · rw [if_pos hskip] at hfac + exact hacc fac hfac + · rw [if_neg hskip] at hfac + simp at hfac + rcases hfac with haccmem | hrec + · exact hacc fac (by simpa using haccmem) + · exact ih childOrder (alpha * gamma ^ j) (gamma ^ ell) child fac + (by simpa using hrec) + exact hfold (List.range ell) #[] (by simp) h + +/-- The top-level smooth splitter emits only represented nonconstant linear factors. -/ +theorem smoothLinearFactorsAlgorithmWith_sound {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (E : BatchEvalContext F) (q : Nat) (generator : F) (schedule : Array Nat) + {p factor : CPolynomial F} + (h : factor ∈ (smoothLinearFactorsAlgorithmWith M D E q generator schedule p).toList) : + IsLinearFactor factor := by + unfold smoothLinearFactorsAlgorithmWith at h + by_cases hzero : (CPolynomial.monicNormalize p == 0 || + CPolynomial.monicNormalize p == 1) = true + · rw [if_pos hzero] at h + simp at h + · rw [if_neg hzero] at h + by_cases hconst : (((CPolynomial.monicNormalize p).coeff 0 == 0) = true) + · rw [if_pos hconst] at h + unfold smoothNonzeroLinearFactorsWith at h + simp at h + rcases h with hzeroFactor | hnonzero + · rcases hzeroFactor with rfl + exact linearFactor_isLinearFactor 0 + · exact smoothCosetLinearFactorsWithSchedule_sound M D E schedule.toList (q - 1) + 1 generator (CPolynomial.monicNormalize (CPolynomial.divX (CPolynomial.monicNormalize p))) + factor (by simpa using hnonzero) + · rw [if_neg hconst] at h + unfold smoothNonzeroLinearFactorsWith at h + simp at h + exact smoothCosetLinearFactorsWithSchedule_sound M D E schedule.toList (q - 1) + 1 generator (CPolynomial.monicNormalize p) factor (by simpa using h) + +/-- A smooth coset split maps the residue class `k % ell` to the child-coset equation. -/ +theorem smooth_coset_split_root_partition_mod {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {alpha gamma : F} {order ell k : Nat} + (hgamma_pow : gamma ^ order = 1) + (hell_dvd : ell ∣ order) : + (alpha * gamma ^ k) ^ (order / ell) = + alpha ^ (order / ell) * (gamma ^ (order / ell)) ^ (k % ell) := by + let n := order / ell + have hellmul : ell * n = order := by + dsimp [n] + rw [Nat.mul_comm] + exact Nat.div_mul_cancel hell_dvd + have hgamma_reduce : gamma ^ (k * n) = gamma ^ ((k % ell) * n) := by + have hkdecomp : k = ell * (k / ell) + k % ell := (Nat.div_add_mod k ell).symm + calc + gamma ^ (k * n) = gamma ^ ((ell * (k / ell) + k % ell) * n) := by + rw [← hkdecomp] + _ = gamma ^ ((k / ell) * order + (k % ell) * n) := by + congr 1 + have hterm : (ell * (k / ell)) * n = (k / ell) * order := by + calc + (ell * (k / ell)) * n = (k / ell) * (ell * n) := by ac_rfl + _ = (k / ell) * order := by rw [hellmul] + rw [Nat.add_mul, hterm] + _ = gamma ^ ((k % ell) * n) := by + rw [pow_add] + rw [Nat.mul_comm (k / ell) order] + rw [pow_mul, hgamma_pow] + simp + calc + (alpha * gamma ^ k) ^ n = alpha ^ n * (gamma ^ k) ^ n := by rw [mul_pow] + _ = alpha ^ n * gamma ^ (k * n) := by rw [pow_mul] + _ = alpha ^ n * gamma ^ ((k % ell) * n) := by rw [hgamma_reduce] + _ = alpha ^ n * gamma ^ (n * (k % ell)) := by rw [Nat.mul_comm (k % ell) n] + _ = alpha ^ n * (gamma ^ n) ^ (k % ell) := by rw [pow_mul] + +/-- A smooth coset split partitions roots according to the child-coset equation. -/ +theorem smooth_coset_split_root_partition {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {alpha gamma x : F} {order ell : Nat} + (hgamma_order : orderOf gamma = order) + (hell_dvd : ell ∣ order) (hell_pos : 0 < ell) + (hx : ∃ k : Nat, k < order ∧ x = alpha * gamma ^ k) : + ∃ j : Nat, j < ell ∧ + x ^ (order / ell) = alpha ^ (order / ell) * (gamma ^ (order / ell)) ^ j := by + rcases hx with ⟨k, _hklt, rfl⟩ + refine ⟨k % ell, Nat.mod_lt k hell_pos, ?_⟩ + apply smooth_coset_split_root_partition_mod + · rw [← hgamma_order] + exact pow_orderOf_eq_one gamma + · exact hell_dvd + +/-- A finite-field generator of order `#F - 1` enumerates all nonzero elements. -/ +theorem exists_generator_pow_of_order_eq_card_sub_one {F : Type*} + [Field F] [Finite F] [DecidableEq F] + {q : Nat} {generator a : F} + (hcard : Nat.card F = q) + (hgen : orderOf generator = q - 1) (ha : a ≠ 0) : + ∃ k : Nat, k < q - 1 ∧ a = generator ^ k := by + have hq_gt_one : 1 < q := by + rw [← hcard] + exact Finite.one_lt_card + have hgen_ne_zero : generator ≠ 0 := by + intro hzero + have hpow : generator ^ (q - 1) = 1 := by + rw [← hgen] + exact pow_orderOf_eq_one generator + rw [hzero] at hpow + have hqsub : 0 < q - 1 := Nat.sub_pos_of_lt hq_gt_one + rcases Nat.exists_eq_add_of_lt hqsub with ⟨m, hm⟩ + rw [hm] at hpow + simp at hpow + let gu : Fˣ := Units.mk0 generator hgen_ne_zero + let au : Fˣ := Units.mk0 a ha + have hgu_order : orderOf (gu : Fˣ) = Nat.card Fˣ := by + rw [← orderOf_units (y := gu)] + simp [gu, hgen, Nat.card_units, hcard] + have htop : Subgroup.zpowers gu = ⊤ := by + rw [← Subgroup.card_eq_iff_eq_top (Subgroup.zpowers gu)] + rw [Nat.card_zpowers, hgu_order] + have hau_z : au ∈ Subgroup.zpowers gu := by + rw [htop] + exact Subgroup.mem_top au + have hau_range := (mem_zpowers_iff_mem_range_orderOf (x := gu) (y := au)).mp hau_z + rcases Finset.mem_image.mp hau_range with ⟨k, hk, hkpow⟩ + refine ⟨k, ?_, ?_⟩ + · have hklt : k < orderOf (gu : Fˣ) := by + simpa using Finset.mem_range.mp hk + rwa [hgu_order, Nat.card_units, hcard] at hklt + · have hval := congrArg (fun u : Fˣ ↦ (u : F)) hkpow + simpa [gu, au] using hval.symm + +/-- Roots of `p` are contained in the explicit coset `alpha * ` of order `order`. -/ +def SmoothCosetInvariant {F : Type*} [Field F] + (alpha gamma : F) (order : Nat) (p : CPolynomial F) : Prop := + alpha ≠ 0 ∧ + 0 < order ∧ + orderOf gamma = order ∧ + ∀ x : F, CPolynomial.eval x p = 0 → + ∃ k : Nat, k < order ∧ x = alpha * gamma ^ k + +/-- Schedule recursion preserves the smooth coset invariant. -/ +theorem smooth_schedule_recursion_preserves_coset_invariant {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + {order ell j : Nat} {alpha gamma : F} {p : CPolynomial F} + (hchild_alpha_ne_zero : alpha * gamma ^ j ≠ 0) + (hchild_order_pos : 0 < order / ell) + (hchild_generator_order : orderOf (gamma ^ ell) = order / ell) + (hchild_roots : + ∀ x : F, + CPolynomial.eval x + (CPolynomial.monicNormalize + (CPolynomial.gcdMonic p + (xPowModWith M D p (order / ell) - + CPolynomial.C + (alpha ^ (order / ell) * (gamma ^ (order / ell)) ^ j)))) = 0 → + ∃ k : Nat, k < order / ell ∧ + x = (alpha * gamma ^ j) * (gamma ^ ell) ^ k) : + SmoothCosetInvariant + (alpha * gamma ^ j) + (gamma ^ ell) + (order / ell) + (CPolynomial.monicNormalize + (CPolynomial.gcdMonic p + (xPowModWith M D p (order / ell) - + CPolynomial.C (alpha ^ (order / ell) * (gamma ^ (order / ell)) ^ j)))) := by + exact ⟨hchild_alpha_ne_zero, hchild_order_pos, hchild_generator_order, hchild_roots⟩ + +/-- Schedule factors divide the current coset order along executable smooth recursion. -/ +def SmoothScheduleDivides : List Nat → Nat → Prop + | [], _ => True + | ell :: rest, order => + if ell ≤ 1 then + SmoothScheduleDivides rest order + else if ell = order then + True + else + ell ∣ order ∧ SmoothScheduleDivides rest (order / ell) + +instance instDecidableSmoothScheduleDivides (schedule : List Nat) (order : Nat) : + Decidable (SmoothScheduleDivides schedule order) := by + induction schedule generalizing order with + | nil => + unfold SmoothScheduleDivides + infer_instance + | cons ell rest ih => + unfold SmoothScheduleDivides + by_cases hellsmall : ell ≤ 1 + · rw [if_pos hellsmall] + exact ih order + · rw [if_neg hellsmall] + by_cases helleq : ell = order + · rw [if_pos helleq] + infer_instance + · rw [if_neg helleq] + infer_instance + +/-- The declared smooth schedule reaches singleton cosets. -/ +theorem smooth_schedule_reaches_singleton {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (ctx : SmoothCyclicRootContext F) : + ctx.schedule.toList.foldl (fun order ell ↦ order / ell) (ctx.q - 1) = 1 := + ctx.schedule_complete + +/-- Path completeness for the schedule-driven smooth coset recursion. -/ +theorem smoothCosetLinearFactorsWithSchedule_complete {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (E : BatchEvalContext F) : + ∀ (schedule : List Nat) (order : Nat) (alpha gamma : F) (p : CPolynomial F) (a : F), + SmoothScheduleDivides schedule order → + gamma ^ order = 1 → + p ≠ 0 → + CPolynomial.eval a p = 0 → + (∃ k : Nat, k < order ∧ a = alpha * gamma ^ k) → + ∃ factor, + factor ∈ + (smoothCosetLinearFactorsWithSchedule M D E schedule order alpha gamma + p).toList ∧ + IsLinearRootFactorCandidate factor a := by + intro schedule + induction schedule with + | nil => + intro order alpha gamma p a _hsched _hgamma hp hroot hcoset + unfold smoothCosetLinearFactorsWithSchedule + let p' := CPolynomial.monicNormalize p + have hp' : p' ≠ 0 := monicNormalize_ne_zero_of_ne_zero hp + have hroot' : CPolynomial.eval a p' = 0 := (monicNormalize_root_iff hp).2 hroot + by_cases hzero : (p' == 0 || p' == 1) = true + · have hcases : p' = 0 ∨ p' = 1 := by + simpa [p'] using hzero + rcases hcases with h0 | h1 + · exact (hp' h0).elim + · rw [h1] at hroot' + rw [eval_one a] at hroot' + exact (one_ne_zero hroot').elim + · rw [if_neg hzero] + by_cases hlin : isRepresentedLinearFactor p' = true + · rw [if_pos hlin] + refine ⟨p', ?_, representedLinearFactor_candidate_of_root hlin hroot'⟩ + simp [p'] + · rw [if_neg hlin] + exact smoothLeafLinearFactors_complete E hcoset hroot' + | cons ell rest ih => + intro order alpha gamma p a hsched hgamma hp hroot hcoset + unfold smoothCosetLinearFactorsWithSchedule + let p' := CPolynomial.monicNormalize p + have hp' : p' ≠ 0 := monicNormalize_ne_zero_of_ne_zero hp + have hroot' : CPolynomial.eval a p' = 0 := (monicNormalize_root_iff hp).2 hroot + by_cases hzero : (p' == 0 || p' == 1) = true + · have hcases : p' = 0 ∨ p' = 1 := by + simpa [p'] using hzero + rcases hcases with h0 | h1 + · exact (hp' h0).elim + · rw [h1] at hroot' + rw [eval_one a] at hroot' + exact (one_ne_zero hroot').elim + · rw [if_neg hzero] + by_cases hlin : isRepresentedLinearFactor p' = true + · rw [if_pos hlin] + refine ⟨p', ?_, representedLinearFactor_candidate_of_root hlin hroot'⟩ + simp [p'] + · rw [if_neg hlin] + by_cases hellsmall : ell ≤ 1 + · rw [if_pos hellsmall] + have hschedRest : SmoothScheduleDivides rest order := by + simpa [SmoothScheduleDivides, hellsmall] using hsched + exact ih order alpha gamma p' a hschedRest hgamma hp' hroot' hcoset + · rw [if_neg hellsmall] + by_cases helleq : ell = order + · rw [if_pos helleq] + exact smoothLeafLinearFactors_complete E hcoset hroot' + · rw [if_neg helleq] + have hschedRest : ell ∣ order ∧ SmoothScheduleDivides rest (order / ell) := by + simpa [SmoothScheduleDivides, hellsmall, helleq] using hsched + rcases hschedRest with ⟨hell_dvd, hschedChild⟩ + rcases hcoset with ⟨k, hk, haeq⟩ + have hell_gt_one : 1 < ell := by omega + have hell_pos : 0 < ell := by omega + let childOrder := order / ell + let tau := gamma ^ childOrder + let xPow := xPowModWith M D p' childOrder + let j := k % ell + have hjlt : j < ell := Nat.mod_lt k hell_pos + have hjmem : j ∈ List.range ell := by + simpa [j] using hjlt + have hellmul : ell * childOrder = order := by + dsimp [childOrder] + rw [Nat.mul_comm] + exact Nat.div_mul_cancel hell_dvd + have hchildGamma : (gamma ^ ell) ^ childOrder = 1 := by + rw [← pow_mul, hellmul, hgamma] + have hchildCoset : ∃ t : Nat, t < childOrder ∧ + a = (alpha * gamma ^ j) * (gamma ^ ell) ^ t := by + refine ⟨k / ell, ?_, ?_⟩ + · dsimp [childOrder] + rw [Nat.div_lt_iff_lt_mul hell_pos] + rw [Nat.mul_comm] + exact lt_of_lt_of_eq hk hellmul.symm + · rw [haeq] + dsimp [j] + have hkdecomp : k % ell + ell * (k / ell) = k := + Nat.mod_add_div k ell + calc + alpha * gamma ^ k = + alpha * gamma ^ (k % ell + ell * (k / ell)) := by + rw [hkdecomp] + _ = alpha * (gamma ^ (k % ell) * gamma ^ (ell * (k / ell))) := by + rw [pow_add] + _ = alpha * gamma ^ (k % ell) * (gamma ^ ell) ^ (k / ell) := by + rw [pow_mul] + ring + let beta := alpha ^ childOrder * tau ^ j + let witness := xPow - CPolynomial.C beta + let child := CPolynomial.monicNormalize (CPolynomial.gcdMonic p' witness) + have hpoweq : a ^ childOrder = beta := by + rw [haeq] + dsimp [beta, tau, j, childOrder] + exact smooth_coset_split_root_partition_mod hgamma hell_dvd + have hwitnessRoot : CPolynomial.eval a witness = 0 := by + dsimp [witness, xPow, beta] + rw [eval_sub, eval_xPowModWith_eq_pow M D hroot' childOrder, eval_C, hpoweq] + ring + have hgcdRoot : CPolynomial.eval a (CPolynomial.gcdMonic p' witness) = 0 := + gcdMonic_root_of_left_right hroot' hwitnessRoot + have hgcdNe : CPolynomial.gcdMonic p' witness ≠ 0 := + gcdMonic_ne_zero_of_left hp' + have hchildNe : child ≠ 0 := monicNormalize_ne_zero_of_ne_zero hgcdNe + have hchildRoot : CPolynomial.eval a child = 0 := + (monicNormalize_root_iff hgcdNe).2 hgcdRoot + have hchildNotOne : child ≠ 1 := by + intro h1 + rw [h1] at hchildRoot + rw [eval_one a] at hchildRoot + exact (one_ne_zero hchildRoot).elim + have hskip : ¬ (child == 0 || child == 1) = true := by + intro hskipTrue + have hcases : child = 0 ∨ child = 1 := by + simpa [child] using hskipTrue + rcases hcases with h0 | h1 + · exact hchildNe h0 + · exact hchildNotOne h1 + have hrec := ih childOrder (alpha * gamma ^ j) (gamma ^ ell) child a + hschedChild hchildGamma hchildNe hchildRoot hchildCoset + have hfold : ∀ (js : List Nat) (acc : Array (CPolynomial F)), + (∃ factor, factor ∈ acc.toList ∧ IsLinearRootFactorCandidate factor a) ∨ + j ∈ js → + ∃ factor, + factor ∈ (js.foldl + (fun factors y ↦ + let beta := alpha ^ childOrder * tau ^ y + let witness := xPow - CPolynomial.C beta + let child := CPolynomial.monicNormalize (CPolynomial.gcdMonic p' witness) + if child == 0 || child == 1 then + factors + else + factors ++ smoothCosetLinearFactorsWithSchedule M D E rest childOrder + (alpha * gamma ^ y) (gamma ^ ell) child) acc).toList ∧ + IsLinearRootFactorCandidate factor a := by + intro js + induction js with + | nil => + intro acc h + rcases h with hacc | hjnil + · exact hacc + · simp at hjnil + | cons y ys ihys => + intro acc h + simp only [List.foldl_cons] + apply ihys + rcases h with hacc | hy + · left + rcases hacc with ⟨factor, hmem, hcand⟩ + refine ⟨factor, ?_, hcand⟩ + by_cases hskipY : (CPolynomial.monicNormalize + (CPolynomial.gcdMonic p' + (xPow - CPolynomial.C (alpha ^ childOrder * tau ^ y))) == 0 || + CPolynomial.monicNormalize + (CPolynomial.gcdMonic p' + (xPow - CPolynomial.C (alpha ^ childOrder * tau ^ y))) == 1) = true + · rw [if_pos hskipY] + exact hmem + · rw [if_neg hskipY] + let childY := CPolynomial.monicNormalize + (CPolynomial.gcdMonic p' + (xPow - CPolynomial.C (alpha ^ childOrder * tau ^ y))) + let tail := smoothCosetLinearFactorsWithSchedule M D E rest childOrder + (alpha * gamma ^ y) (gamma ^ ell) childY + simpa [childY, tail] using Array.mem_append_left tail + (by simpa using hmem) + · simp at hy + rcases hy with hyj | hyrest + · subst y + left + rcases hrec with ⟨factor, hmem, hcand⟩ + refine ⟨factor, ?_, hcand⟩ + rw [if_neg hskip] + simpa [child, witness, beta] using + Array.mem_append_right acc (by simpa using hmem) + · right + exact hyrest + exact hfold (List.range ell) #[] (Or.inr hjmem) + +/-- Completeness of the top-level smooth linear-factor splitter. -/ +theorem smoothLinearFactorsAlgorithmWith_complete {F : Type*} + [Field F] [Finite F] [DecidableEq F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (E : BatchEvalContext F) (q : Nat) (generator : F) (schedule : Array Nat) + (hcard : Nat.card F = q) + (hgenerator : orderOf generator = q - 1) + (hschedule : SmoothScheduleDivides schedule.toList (q - 1)) + {p : CPolynomial F} {a : F} + (hp : p ≠ 0) (hroot : CPolynomial.eval a p = 0) : + ∃ factor, + factor ∈ (smoothLinearFactorsAlgorithmWith M D E q generator schedule p).toList ∧ + IsLinearRootFactorCandidate factor a := by + unfold smoothLinearFactorsAlgorithmWith + let p' := CPolynomial.monicNormalize p + have hp' : p' ≠ 0 := monicNormalize_ne_zero_of_ne_zero hp + have hroot' : CPolynomial.eval a p' = 0 := (monicNormalize_root_iff hp).2 hroot + by_cases hzero : (p' == 0 || p' == 1) = true + · have hcases : p' = 0 ∨ p' = 1 := by + simpa [p'] using hzero + rcases hcases with h0 | h1 + · exact (hp' h0).elim + · rw [h1] at hroot' + rw [eval_one a] at hroot' + exact (one_ne_zero hroot').elim + · rw [if_neg hzero] + by_cases ha0 : a = 0 + · subst a + have hcoeff : p'.coeff 0 = 0 := by + simpa [eval_zero_eq_coeff_zero] using hroot' + have hconst : ((p'.coeff 0 == 0) = true) := by + simp [hcoeff] + rw [if_pos hconst] + refine ⟨CPolynomial.linearFactor (0 : F), ?_, linearFactor_isRootFactorCandidate 0⟩ + unfold smoothNonzeroLinearFactorsWith + simp + · have hgenPow : generator ^ (q - 1) = 1 := by + rw [← hgenerator] + exact pow_orderOf_eq_one generator + have hcoset : ∃ k : Nat, k < q - 1 ∧ a = (1 : F) * generator ^ k := by + rcases exists_generator_pow_of_order_eq_card_sub_one hcard hgenerator ha0 with + ⟨k, hk, ha⟩ + exact ⟨k, hk, by simp [ha]⟩ + by_cases hconst : ((p'.coeff 0 == 0) = true) + · rw [if_pos hconst] + have hcoeff : p'.coeff 0 = 0 := by simpa using hconst + have hdivRoot : CPolynomial.eval a (CPolynomial.divX p') = 0 := + eval_divX_eq_zero_of_ne_zero_root ha0 hcoeff hroot' + have hdivNe : CPolynomial.divX p' ≠ 0 := + divX_ne_zero_of_ne_zero_coeff_zero hp' hcoeff + have hchildRoot : + CPolynomial.eval a (CPolynomial.monicNormalize (CPolynomial.divX p')) = 0 := + (monicNormalize_root_iff hdivNe).2 hdivRoot + have hchildNe : CPolynomial.monicNormalize (CPolynomial.divX p') ≠ 0 := + monicNormalize_ne_zero_of_ne_zero hdivNe + rcases smoothCosetLinearFactorsWithSchedule_complete M D E schedule.toList (q - 1) + (1 : F) generator (CPolynomial.monicNormalize (CPolynomial.divX p')) a + hschedule hgenPow hchildNe hchildRoot hcoset with ⟨factor, hmem, hcand⟩ + refine ⟨factor, ?_, hcand⟩ + unfold smoothNonzeroLinearFactorsWith + simpa [p'] using Array.mem_append_right (#[CPolynomial.linearFactor (0 : F)]) + (by simpa using hmem) + · rw [if_neg hconst] + rcases smoothCosetLinearFactorsWithSchedule_complete M D E schedule.toList (q - 1) + (1 : F) generator p' a hschedule hgenPow hp' hroot' hcoset with + ⟨factor, hmem, hcand⟩ + refine ⟨factor, ?_, hcand⟩ + unfold smoothNonzeroLinearFactorsWith + simpa using hmem + +end FiniteField + +end Roots + +end CPolynomial + +end CompPoly diff --git a/CompPoly/Univariate/Roots/Splitter.lean b/CompPoly/Univariate/Roots/Splitter.lean new file mode 100644 index 00000000..97c99b05 --- /dev/null +++ b/CompPoly/Univariate/Roots/Splitter.lean @@ -0,0 +1,57 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Raw.Modular +import CompPoly.Univariate.Roots.Context + +/-! +# Linear-Factor Splitter Helpers + +Shared executable helpers for finite-field linear-factor splitters. +-/ + +namespace CompPoly + +namespace CPolynomial + +namespace Roots + +namespace FiniteField + +/-- Boolean recognizer for represented nonconstant linear factors. -/ +def isRepresentedLinearFactor {F : Type*} [Field F] [BEq F] + (p : CPolynomial F) : Bool := + decide (p.val.size ≤ 2) && !(p.coeff 1 == 0) + +/-- The represented-linear Boolean recognizer is sound. -/ +theorem isRepresentedLinearFactor_sound {F : Type*} + [Field F] [BEq F] [LawfulBEq F] + {p : CPolynomial F} (h : isRepresentedLinearFactor p = true) : + IsLinearFactor p := by + unfold isRepresentedLinearFactor at h + simp at h + simpa [IsLinearFactor, CPolynomial.coeff, CPolynomial.Raw.coeff] using h + +/-- Return `p` as a singleton array if it is represented as a nonconstant linear factor. -/ +def representedLinearFactorArray {F : Type*} [Field F] [BEq F] + (p : CPolynomial F) : Array (CPolynomial F) := + if isRepresentedLinearFactor p then #[p] else #[] + +/-- `X^exponent mod modulus`, lifted back to canonical polynomials. -/ +def xPowModWith {F : Type*} [Field F] [BEq F] [LawfulBEq F] + (M : CPolynomial.Raw.MulContext F) (D : CPolynomial.Raw.ModContext F) + (modulus : CPolynomial F) (exponent : Nat) : CPolynomial F := + CPolynomial.ofArray + (CPolynomial.Raw.powModWith M D modulus.val + (CPolynomial.Raw.X : CPolynomial.Raw F) exponent) + +end FiniteField + +end Roots + +end CPolynomial + +end CompPoly diff --git a/CompPoly/Univariate/ToPoly/Core.lean b/CompPoly/Univariate/ToPoly/Core.lean index df15dcf4..c0ba9f71 100644 --- a/CompPoly/Univariate/ToPoly/Core.lean +++ b/CompPoly/Univariate/ToPoly/Core.lean @@ -206,6 +206,12 @@ theorem trim_toImpl [LawfulBEq R] (p : R[X]) : p.toImpl.trim = p.toImpl := by end Raw +/-- `ofArray` preserves the raw polynomial's `toPoly` image. -/ +theorem ofArray_toPoly [LawfulBEq R] (p : CPolynomial.Raw R) : + (CPolynomial.ofArray p).toPoly = p.toPoly := by + unfold CPolynomial.ofArray + exact Raw.toPoly_trim + /-- On canonical polynomials, `toImpl` is a left-inverse of `toPoly`. This shows `toPoly` is a bijection from `CPolynomial R` to `Polynomial R`. -/ diff --git a/CompPoly/Univariate/ToPoly/Impl.lean b/CompPoly/Univariate/ToPoly/Impl.lean index 97aabc1d..c1e8ada6 100644 --- a/CompPoly/Univariate/ToPoly/Impl.lean +++ b/CompPoly/Univariate/ToPoly/Impl.lean @@ -51,6 +51,11 @@ theorem eval_toPoly [BEq R] [LawfulBEq R] (x : R) (p : CPolynomial R) : · rw [ Raw.eval_toPoly_eq_eval ]; rfl · convert Raw.eval_toPoly_eq_eval x p.val +/-- Evaluation of a constant computable polynomial. -/ +theorem eval_C [BEq R] [LawfulBEq R] (a c : R) : + CPolynomial.eval a (CPolynomial.C c) = c := by + rw [CPolynomial.eval_toPoly, CPolynomial.C_toPoly, Polynomial.eval_C] + /-- Raw.eval₂ is correct wrt the Mathlib spec. -/ theorem Raw.eval₂_toPoly {S : Type*} [Semiring S] (f : R →+* S) (x : S) (p : CPolynomial.Raw R) : @@ -78,6 +83,13 @@ theorem coeff_toPoly [BEq R] [LawfulBEq R] (p : CPolynomial R) (i : ℕ) : unfold toPoly coeff simp [Raw.coeff_toPoly] +/-- Evaluation at zero returns the constant coefficient. -/ +theorem eval_zero_eq_coeff_zero [BEq R] [LawfulBEq R] + (p : CPolynomial R) : CPolynomial.eval 0 p = p.coeff 0 := by + rw [CPolynomial.eval_toPoly] + rw [← Polynomial.coeff_zero_eq_eval_zero p.toPoly] + exact (CPolynomial.coeff_toPoly p 0).symm + /-- CPolynomial.divX is correct wrt the Mathlib spec. -/ theorem divX_toPoly [BEq R] [LawfulBEq R] (p : CPolynomial R) : (divX p).toPoly = p.toPoly.divX := by @@ -262,6 +274,52 @@ theorem eval_ext Polynomial.eq_zero_of_natDegree_lt_card_of_eval_eq_zero' r.toPoly T heval_zero (lt_of_le_of_lt hrPolyDeg hTcard) +/-- Evaluation preserves subtraction. -/ +theorem eval_sub [Ring R] [BEq R] [LawfulBEq R] + (a : R) (p q : CPolynomial R) : + CPolynomial.eval a (p - q) = CPolynomial.eval a p - CPolynomial.eval a q := by + rw [CPolynomial.eval_toPoly, CPolynomial.toPoly_sub, Polynomial.eval_sub, + ← CPolynomial.eval_toPoly, ← CPolynomial.eval_toPoly] + +/-- Evaluation of the constant one computable polynomial. -/ +theorem eval_one [Semiring R] [BEq R] [LawfulBEq R] [Nontrivial R] + (a : R) : CPolynomial.eval a (1 : CPolynomial R) = 1 := by + rw [CPolynomial.eval_toPoly, CPolynomial.toPoly_one, Polynomial.eval_one] + +/-- Dividing by `X` preserves a nonzero root when the constant coefficient vanishes. -/ +theorem eval_divX_eq_zero_of_ne_zero_root [Field R] [BEq R] [LawfulBEq R] + {p : CPolynomial R} {a : R} (ha : a ≠ 0) + (hcoeff : p.coeff 0 = 0) (hroot : CPolynomial.eval a p = 0) : + CPolynomial.eval a (CPolynomial.divX p) = 0 := by + have hdecomp := CPolynomial.X_mul_divX_add (p := p) + have hroot' : + CPolynomial.eval a (CPolynomial.X * CPolynomial.divX p + CPolynomial.C (p.coeff 0)) = + 0 := by + rw [← hdecomp] + exact hroot + rw [CPolynomial.eval_toPoly, CPolynomial.toPoly_add, CPolynomial.toPoly_mul, + CPolynomial.X_toPoly, CPolynomial.C_toPoly, Polynomial.eval_add, Polynomial.eval_mul, + Polynomial.eval_X, Polynomial.eval_C, ← CPolynomial.eval_toPoly] at hroot' + rw [hcoeff] at hroot' + simp at hroot' + rcases hroot' with hzero | hdivRoot + · exact (ha hzero).elim + · exact hdivRoot + +/-- If a nonzero polynomial has zero constant coefficient, its quotient by `X` is nonzero. -/ +theorem divX_ne_zero_of_ne_zero_coeff_zero [Field R] [BEq R] [LawfulBEq R] + {p : CPolynomial R} (hp : p ≠ 0) (hcoeff : p.coeff 0 = 0) : + CPolynomial.divX p ≠ 0 := by + intro hdiv + apply hp + have hC0 : CPolynomial.C (0 : R) = 0 := by + apply (CPolynomial.eq_iff_coeff).2 + intro i + simp + rw [CPolynomial.X_mul_divX_add (p := p), hdiv, hcoeff] + rw [hC0] + simp + namespace Raw theorem eval_sub_C_mul_X_pow_trim_eq_self_of_eval_eq_zero diff --git a/lake-manifest.json b/lake-manifest.json index 18f01401..b547213f 100644 --- a/lake-manifest.json +++ b/lake-manifest.json @@ -5,17 +5,17 @@ "type": "git", "subDir": null, "scope": "leanprover-community", - "rev": "5450b53e5ddc75d46418fabb605edbf36bd0beb6", + "rev": "c5ea00351c28e24afc9f0f84379aa41082b1188f", "name": "mathlib", "manifestFile": "lake-manifest.json", - "inputRev": "v4.30.0-rc2", + "inputRev": "v4.30.0", "inherited": false, "configFile": "lakefile.lean"}, {"url": "https://github.com/leanprover-community/plausible", "type": "git", "subDir": null, "scope": "leanprover-community", - "rev": "86210d4ad1b08b086d0bd638637a75246523dbb8", + "rev": "a456461b368b71d2accd95234832cd9c174b5437", "name": "plausible", "manifestFile": "lake-manifest.json", "inputRev": "main", @@ -35,7 +35,7 @@ "type": "git", "subDir": null, "scope": "leanprover-community", - "rev": "cdab3938ccabbdb044be6896e251b5814bec932e", + "rev": "515cf9d0c00ece5e661f6de4326a53dedc1e8ea1", "name": "importGraph", "manifestFile": "lake-manifest.json", "inputRev": "main", @@ -45,50 +45,50 @@ "type": "git", "subDir": null, "scope": "leanprover-community", - "rev": "2db6054a44326f8c0230ee0570e2ddb894816511", + "rev": "a84b3e2475d5c5ab979567b1ad8aea21b764bcf8", "name": "proofwidgets", "manifestFile": "lake-manifest.json", - "inputRev": "v0.0.98", + "inputRev": "v0.0.99", "inherited": true, "configFile": "lakefile.lean"}, {"url": "https://github.com/leanprover-community/aesop", "type": "git", "subDir": null, "scope": "leanprover-community", - "rev": "f0c6e183ea26531e82773feb4b73ab6595ca17a5", + "rev": "558915ae105bfd8074e22d597613d1961822adc2", "name": "aesop", "manifestFile": "lake-manifest.json", - "inputRev": "v4.30.0-rc2", + "inputRev": "v4.30.0", "inherited": true, "configFile": "lakefile.toml"}, {"url": "https://github.com/leanprover-community/quote4", "type": "git", "subDir": null, "scope": "leanprover-community", - "rev": "1cc7e819b9b9bc1e87c9edcccb62e0269e00a809", + "rev": "a6e6c34c4ef182f83b219a3a5a385f51f44bdc4c", "name": "Qq", "manifestFile": "lake-manifest.json", - "inputRev": "v4.30.0-rc2", + "inputRev": "v4.30.0", "inherited": true, "configFile": "lakefile.toml"}, {"url": "https://github.com/leanprover-community/batteries", "type": "git", "subDir": null, "scope": "leanprover-community", - "rev": "5c57f3857ba81924a88b2cdf4f062e34ec04ff11", + "rev": "32dc18cde3684679f3c003de608743b57498c56f", "name": "batteries", "manifestFile": "lake-manifest.json", - "inputRev": "v4.30.0-rc2", + "inputRev": "main", "inherited": true, "configFile": "lakefile.toml"}, {"url": "https://github.com/leanprover/lean4-cli", "type": "git", "subDir": null, "scope": "leanprover", - "rev": "13567aed1ac4f12aea9484178e07e51f8c9f7658", + "rev": "6b907cf12b2e445ccb7c24bc208ef04a1f39e84c", "name": "Cli", "manifestFile": "lake-manifest.json", - "inputRev": "v4.30.0-rc2", + "inputRev": "v4.30.0", "inherited": true, "configFile": "lakefile.toml"}], "name": "CompPoly", diff --git a/lakefile.lean b/lakefile.lean index 7e3cefc4..16d958cd 100644 --- a/lakefile.lean +++ b/lakefile.lean @@ -6,7 +6,7 @@ package CompPoly where version := v!"0.1.0" testDriver := "CompPolyTests" -require "leanprover-community" / mathlib @ git "v4.30.0-rc2" +require "leanprover-community" / mathlib @ git "v4.30.0" @[default_target] lean_lib CompPoly diff --git a/lean-toolchain b/lean-toolchain index 635bb953..af9e5d33 100644 --- a/lean-toolchain +++ b/lean-toolchain @@ -1 +1 @@ -leanprover/lean4:v4.30.0-rc2 \ No newline at end of file +leanprover/lean4:v4.30.0 diff --git a/tests/CompPolyTests.lean b/tests/CompPolyTests.lean index 9d4818d5..c8cb52dc 100644 --- a/tests/CompPolyTests.lean +++ b/tests/CompPolyTests.lean @@ -28,4 +28,6 @@ import CompPolyTests.Univariate.NTT.FastMul import CompPolyTests.Univariate.NTT.Forward import CompPolyTests.Univariate.NTT.Inverse import CompPolyTests.Univariate.Raw +import CompPolyTests.Univariate.Roots.Enumeration +import CompPolyTests.Univariate.Roots.FiniteField import CompPolyTests.Univariate.ToPoly diff --git a/tests/CompPolyTests/Univariate/Roots/Enumeration.lean b/tests/CompPolyTests/Univariate/Roots/Enumeration.lean new file mode 100644 index 00000000..72718c7c --- /dev/null +++ b/tests/CompPolyTests/Univariate/Roots/Enumeration.lean @@ -0,0 +1,84 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Roots.Enumeration +import Mathlib.Algebra.Field.ZMod + +/-! +# Exhaustive Enumeration Root Tests + +Executable coverage for lazy finite-field enumeration adapters over small +prime fields. +-/ + +namespace CompPolyTests + +open CompPoly +open CompPoly.CPolynomial.Roots.FiniteField + +namespace Univariate.Roots.Enumeration + +abbrev F3 := ZMod 3 +abbrev F5 := ZMod 5 + +instance : Fact (Nat.Prime 3) := + ⟨by decide⟩ + +instance : Fact (Nat.Prime 5) := + ⟨by decide⟩ + +private def f3Elements : Array F3 := + #[0, 1, 2] + +private theorem f3Elements_complete : ContainsAllFieldElements f3Elements := by + unfold ContainsAllFieldElements + intro a + fin_cases a <;> decide + +private def f3Enumeration : FieldEnumeration F3 := + fieldEnumerationOfArray f3Elements f3Elements_complete + +private def f3RootsPolynomial : CPolynomial F3 := + CPolynomial.linearFactor (0 : F3) * CPolynomial.linearFactor (2 : F3) + +private def f3Roots : Array F3 := + rootsInFieldByEnumeration f3Enumeration f3RootsPolynomial + +#guard f3Roots == #[(0 : F3), (2 : F3)] + +private def f5Elements : Array F5 := + #[0, 1, 2, 3, 4] + +private theorem f5Elements_complete : ContainsAllFieldElements f5Elements := by + unfold ContainsAllFieldElements + intro a + fin_cases a <;> decide + +private def f5Enumeration : FieldEnumeration F5 := + fieldEnumerationOfArray f5Elements f5Elements_complete + +private def f5RootsPolynomial : CPolynomial F5 := + CPolynomial.linearFactor (1 : F5) * + CPolynomial.linearFactor (4 : F5) * + CPolynomial.linearFactor (4 : F5) + +private def f5Roots : Array F5 := + rootsInFieldByEnumeration f5Enumeration f5RootsPolynomial + +#guard f5Roots == #[(1 : F5), (4 : F5)] + +private def f5EnumerationSplitter : LinearFactorProductSplitter F5 := + enumeratingLinearFactorProductSplitter f5Enumeration + +private def f5EnumerationFactors : Array (CPolynomial F5) := + f5EnumerationSplitter.splitLinearFactors 5 f5RootsPolynomial + +#guard (CPolynomial.rootsFromLinearFactors f5RootsPolynomial f5EnumerationFactors).contains (1 : F5) +#guard (CPolynomial.rootsFromLinearFactors f5RootsPolynomial f5EnumerationFactors).contains (4 : F5) + +end Univariate.Roots.Enumeration + +end CompPolyTests diff --git a/tests/CompPolyTests/Univariate/Roots/FiniteField.lean b/tests/CompPolyTests/Univariate/Roots/FiniteField.lean new file mode 100644 index 00000000..b00c6fb7 --- /dev/null +++ b/tests/CompPolyTests/Univariate/Roots/FiniteField.lean @@ -0,0 +1,87 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin +-/ + +import CompPoly.Univariate.Roots.Correctness +import Mathlib.Algebra.Field.ZMod + +/-! +# Finite-Field Univariate Root Tests + +Focused executable coverage for the generic finite-field root backend. +-/ + +namespace CompPolyTests + +open CompPoly +open CompPoly.CPolynomial.Roots.FiniteField + +namespace Univariate.Roots.FiniteField + +abbrev F3 := ZMod 3 +abbrev F5 := ZMod 5 + +instance : Fact (Nat.Prime 3) := + ⟨by decide⟩ + +instance : Fact (Nat.Prime 5) := + ⟨by decide⟩ + +private def f3Ctx : FiniteFieldContext F3 where + q := 3 + finite := by infer_instance + card_eq := by + simp [F3, Nat.card_eq_fintype_card, ZMod.card] + frobenius_fixed := by decide + +private def f5Ctx : FiniteFieldContext F5 where + q := 5 + finite := by infer_instance + card_eq := by + simp [F5, Nat.card_eq_fintype_card, ZMod.card] + frobenius_fixed := by decide + +private def f3Repeated : CPolynomial F3 := + CPolynomial.linearFactor (1 : F3) * + CPolynomial.linearFactor (2 : F3) * + CPolynomial.linearFactor (2 : F3) + +private def f3SquarefreeProduct : CPolynomial F3 := + CPolynomial.linearFactor (1 : F3) * + CPolynomial.linearFactor (2 : F3) + +private def f3RootProduct : CPolynomial F3 := + finiteFieldRootProduct f3Ctx f3Repeated + +#guard CPolynomial.evalHorner (1 : F3) f3RootProduct == 0 +#guard CPolynomial.evalHorner (2 : F3) f3RootProduct == 0 +#guard CPolynomial.evalHorner (0 : F3) f3RootProduct != 0 +#guard f3RootProduct == f3SquarefreeProduct + +private def f3NoRoot : CPolynomial F3 := + CPolynomial.ofArray #[(1 : F3), 0, 1] + +#guard finiteFieldRootProduct f3Ctx f3NoRoot == 1 + +private def f5MultipleRoots : CPolynomial F5 := + CPolynomial.linearFactor (1 : F5) * + CPolynomial.linearFactor (1 : F5) * + CPolynomial.linearFactor (3 : F5) + +private def f5LinearProduct : CPolynomial F5 := + CPolynomial.linearFactor (1 : F5) * + CPolynomial.linearFactor (3 : F5) + +private def f5RootProduct : CPolynomial F5 := + finiteFieldRootProduct f5Ctx f5MultipleRoots + +#guard CPolynomial.evalHorner (1 : F5) f5RootProduct == 0 +#guard CPolynomial.evalHorner (3 : F5) f5RootProduct == 0 +#guard CPolynomial.evalHorner (2 : F5) f5RootProduct != 0 +#guard f5RootProduct == f5LinearProduct + +end Univariate.Roots.FiniteField + +end CompPolyTests