diff --git a/CompPoly.lean b/CompPoly.lean index f5801449..2894a260 100644 --- a/CompPoly.lean +++ b/CompPoly.lean @@ -1,6 +1,10 @@ import CompPoly.Bivariate.Basic import CompPoly.Bivariate.CMvEquiv import CompPoly.Bivariate.ToPoly +import CompPoly.CodingTheory.GuruswamiSudan +import CompPoly.CodingTheory.GuruswamiSudan.Counterexamples +import CompPoly.CodingTheory.GuruswamiSudan.Correctness +import CompPoly.CodingTheory.GuruswamiSudanFast import CompPoly.Data.Array.Lemmas import CompPoly.Data.Classes.DCast import CompPoly.Data.ExtTreeMap.DTreeMap diff --git a/CompPoly/CodingTheory/GuruswamiSudan.lean b/CompPoly/CodingTheory/GuruswamiSudan.lean new file mode 100644 index 00000000..9eeff0aa --- /dev/null +++ b/CompPoly/CodingTheory/GuruswamiSudan.lean @@ -0,0 +1,102 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: OpenAI Codex +-/ + +import CompPoly.CodingTheory.GuruswamiSudan.Generic +import CompPoly.Fields.KoalaBear + +/-! +# Executable Guruswami-Sudan Decoder + +Canonical `KoalaBear.Field` instantiation of the generic executable +Guruswami-Sudan decoder. +-/ + +namespace CompPoly +namespace CodingTheory +namespace GuruswamiSudan + +abbrev F := KoalaBear.Field +abbrev UniPoly := Generic.UniPoly F +abbrev Bivariate := Generic.Bivariate F +abbrev ReedSolomonCode := Generic.ReedSolomonCode F +abbrev DecodeResult := Generic.DecodeResult F + +def maxCandidatesPerDepth : Nat := Generic.maxCandidatesPerDepth +def maxTotalRoots : Nat := Generic.maxTotalRoots + +namespace NatUtil + +export Generic.NatUtil (floorSqrt ceilSqrt ceilSqrtDiv binomial) + +end NatUtil + +namespace UniPoly + +export Generic.UniPoly + (zero one isZero trim coeff degree leadingCoeff ofNatArray neg add sub scale mul + mulXPow evaluate pow monomial longDivRem? coeffsEq) + +end UniPoly + +namespace Bivariate + +export Generic.Bivariate + (zero trimCoeffs ofCoeffs isZero yDegree maxXDegree coeff evaluate weightedDegree + evaluateYPolynomial fromMonomials) + +end Bivariate + +namespace ReedSolomonCode + +def withDomain (domain : Array F) (k : Nat) : ReedSolomonCode := + Generic.ReedSolomonCode.withDomain domain k + +def consecutiveDomain (n : Nat) : Array F := + Generic.ReedSolomonCode.consecutiveDomain (F := F) n + +def withConsecutiveDomain (n k : Nat) : ReedSolomonCode := + Generic.ReedSolomonCode.withConsecutiveDomain (F := F) n k + +def pow2Domain (logN : Nat) : Array F := + let omega := KoalaBear.twoAdicGenerators.toArray.getD logN (1 : F) + Generic.ReedSolomonCode.pow2Domain (F := F) logN omega + +def withRootsOfUnityDomain (logN k : Nat) : ReedSolomonCode := + let omega := KoalaBear.twoAdicGenerators.toArray.getD logN (1 : F) + Generic.ReedSolomonCode.withRootsOfUnityDomain (F := F) logN k omega + +def encodePolynomial (code : ReedSolomonCode) (poly : UniPoly) : Array F := + Generic.ReedSolomonCode.encodePolynomial code poly + +def encode (code : ReedSolomonCode) (message : Array F) : Array F := + Generic.ReedSolomonCode.encode code message + +end ReedSolomonCode + +export Generic + (gsDecodingRadius johnsonListBound? countMonomials chooseParameters + monomialsBelowWeightedDegree findKernelVector interpolateWithMultiplicity + extractSmallValue? lagrangeInterpolateAtZeroWithPoints? + lagrangeInterpolatePolynomial? substituteAndDivide findRootsLinearY + appendUniquePoly appendUniqueField enumerateSmallPolys + trySmallIntegerPolynomials tryDirectRoots tryInterpolatedCandidates + findUnivariateRootsWithHints findUnivariateRootsWithHintsAndDomain + rrSearchWithDomain findPolynomialRootsWithDomain agreement + gsListDecodeWithMultiplicity gsListDecode introduceErrors + introduceErrorsAtPositions) + +def polyToNats (p : UniPoly) : Array Nat := + Generic.polyToNats (F := F) ZMod.val p + +def valuesToNats (xs : Array F) : Array Nat := + Generic.valuesToNats (F := F) ZMod.val xs + +def candidatesToNats (xs : Array UniPoly) : Array (Array Nat) := + Generic.candidatesToNats (F := F) ZMod.val xs + +end GuruswamiSudan +end CodingTheory +end CompPoly diff --git a/CompPoly/CodingTheory/GuruswamiSudan/Correctness.lean b/CompPoly/CodingTheory/GuruswamiSudan/Correctness.lean new file mode 100644 index 00000000..27f2ed1f --- /dev/null +++ b/CompPoly/CodingTheory/GuruswamiSudan/Correctness.lean @@ -0,0 +1,197 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: OpenAI Codex +-/ + +import CompPoly.CodingTheory.GuruswamiSudan.Generic +import CompPoly.Univariate.ToPoly.Core + +/-! +# Executable Correctness Milestones for Guruswami-Sudan + +This file records proof targets that are true for the current executable +Guruswami-Sudan implementation. The root finder in `Generic.lean` intentionally +matches the Lambdaworks heuristic, so full completeness is stated conditionally +on that root finder returning every polynomial root of the interpolation +polynomial. +-/ + +namespace CompPoly +namespace CodingTheory +namespace GuruswamiSudan +namespace Generic + +variable {F : Type} [Field F] + +namespace UniPoly + +/-- Mathlib polynomial denoted by the executable GS coefficient array. -/ +noncomputable def toPolynomial (p : UniPoly F) : Polynomial F := + CPolynomial.Raw.toPoly (p : CPolynomial.Raw F) + +/-- The executable Horner evaluator agrees with Mathlib polynomial evaluation. -/ +theorem evaluate_eq_toPolynomial_eval (p : UniPoly F) (x : F) : + UniPoly.evaluate p x = (toPolynomial p).eval x := by + calc + UniPoly.evaluate p x = + CPolynomial.Raw.eval₂Horner (RingHom.id F) x (p : CPolynomial.Raw F) := rfl + _ = CPolynomial.Raw.eval x (p : CPolynomial.Raw F) := + CPolynomial.Raw.eval₂Horner_eq_eval₂ (RingHom.id F) x (p : CPolynomial.Raw F) + _ = (toPolynomial p).eval x := + (CPolynomial.Raw.eval_toPoly_eq_eval x (p : CPolynomial.Raw F)).symm + +end UniPoly + +variable [BEq F] [Inhabited F] + +/-- A returned candidate satisfies the executable Reed-Solomon list-decoding +acceptance predicate: degree below `k` and enough agreement with the received +word. -/ +def IsAcceptedCandidate + (code : ReedSolomonCode F) (received : Array F) (f : UniPoly F) : Prop := + UniPoly.degree f < code.k ∧ + agreement received code.domain f ≥ code.n - gsDecodingRadius code.n code.k + +/-- Dot product used to state the executable interpolation linear system. -/ +def dotProduct (row vector : Array F) : F := Id.run do + let mut acc : F := 0 + for i in [0:min row.size vector.size] do + acc := acc + row[i]! * vector[i]! + pure acc + +/-- A vector is in the right kernel of every row of an executable matrix. -/ +def VectorInKernel (matrix : Array (Array F)) (vector : Array F) : Prop := + ∀ row ∈ matrix, dotProduct row vector = 0 + +/-- The executable interpolation system attached to the multiplicity step. This +predicate intentionally talks about the concrete matrix built by +`interpolationMatrix`; connecting those rows to Hasse derivatives is the next +mathematical layer above this executable milestone. -/ +def SatisfiesInterpolationSystem + (domain received : Array F) (m d k : Nat) (q : Bivariate F) : Prop := + let monomials := monomialsBelowWeightedDegree d k + let matrix := interpolationMatrix domain received m d k + let solution := findKernelVector matrix monomials.size + q = Bivariate.fromMonomials monomials solution ∧ VectorInKernel matrix solution + +/-- The precise linear-algebra fact needed from `findKernelVector` for +interpolation correctness. -/ +def InterpolationKernelCorrect + (domain received : Array F) (m d k : Nat) : Prop := + let monomials := monomialsBelowWeightedDegree d k + let matrix := interpolationMatrix domain received m d k + VectorInKernel matrix (findKernelVector matrix monomials.size) + +/-- Soundness of the executable decoder's final filter. This theorem does not +claim that the heuristic root finder found every possible root; it says every +candidate that survives the decoder is genuinely accepted by the executable +agreement predicate. -/ +theorem gsListDecode_sound + (code : ReedSolomonCode F) (received : Array F) (f : UniPoly F) + (hsize : received.size = code.n) + (hmem : f ∈ (gsListDecode code received).candidates) : + IsAcceptedCandidate code received f := by + simp only [gsListDecode] at hmem + have hneq : (received.size != code.n) = false := by simp [hsize] + simp [hneq] at hmem + exact ⟨hmem.2.1, (Nat.sub_le_iff_le_add).2 hmem.2.2⟩ + +/-- Interpolation correctness reduced to the executable linear-algebra kernel +claim. This theorem has no proof placeholder and states exactly the remaining +obligation: the row-reduction routine must return a vector annihilating the +interpolation matrix. -/ +theorem interpolateWithMultiplicity_correct_of_kernel + (domain received : Array F) (m d k : Nat) + (hkernel : InterpolationKernelCorrect domain received m d k) : + SatisfiesInterpolationSystem domain received m d k + (interpolateWithMultiplicity domain received m d k) := by + simpa [SatisfiesInterpolationSystem, InterpolationKernelCorrect, + interpolateWithMultiplicity] using hkernel + +/-- Completeness assumption for the current heuristic root finder at one +particular interpolation polynomial. This is false in general, but useful as the +explicit hypothesis under which the rest of GS completeness can be stated. -/ +def RootFinderCompleteFor + (q : Bivariate F) (maxDegree : Nat) (hintValues domain : Array F) : Prop := + ∀ f : UniPoly F, + UniPoly.degree f < maxDegree → + UniPoly.isZero (Bivariate.evaluateYPolynomial q f) → + f ∈ findPolynomialRootsWithDomain q maxDegree hintValues domain + +/-- Soundness assumption for one concrete call to the heuristic root finder: +every returned candidate is really a `Y`-root of the bivariate polynomial and +has the requested degree bound. -/ +def RootFinderSoundFor + (q : Bivariate F) (maxDegree : Nat) (hintValues domain : Array F) : Prop := + ∀ f : UniPoly F, + f ∈ findPolynomialRootsWithDomain q maxDegree hintValues domain → + UniPoly.degree f < maxDegree ∧ + UniPoly.isZero (Bivariate.evaluateYPolynomial q f) + +/-- Exact per-case adequacy of the heuristic root finder. This is the assumption +that makes the executable decoder exact for the particular interpolation +polynomial used in this decoding run. -/ +def RootFinderAdequateFor + (q : Bivariate F) (maxDegree : Nat) (hintValues domain : Array F) : Prop := + RootFinderSoundFor q maxDegree hintValues domain ∧ + RootFinderCompleteFor q maxDegree hintValues domain + +/-- Conditional completeness of the executable decoder after the algebraic GS +step has established that the message polynomial is a root of the interpolated +`Q`, and assuming the heuristic root finder is complete for that `Q`. -/ +theorem gsListDecode_complete_of_rootFinderComplete + (code : ReedSolomonCode F) (received : Array F) (f : UniPoly F) + (hsize : received.size = code.n) + (m d : Nat) + (hparams : chooseParameters code.n code.k = (m, d)) + (hrootComplete : + let q := interpolateWithMultiplicity code.domain received m d code.k + RootFinderCompleteFor q code.k received code.domain) + (hroot : + let q := interpolateWithMultiplicity code.domain received m d code.k + UniPoly.isZero (Bivariate.evaluateYPolynomial q f)) + (hdeg : UniPoly.degree f < code.k) + (hagree : + agreement received code.domain f ≥ + code.n - gsDecodingRadius code.n code.k) : + f ∈ (gsListDecode code received).candidates := by + simp only [gsListDecode] + have hneq : (received.size != code.n) = false := by simp [hsize] + simp [hneq, hparams] + refine ⟨?_, hdeg, ?_⟩ + · exact hrootComplete f hdeg hroot + · exact (Nat.sub_le_iff_le_add).1 hagree + +/-- Under exact per-case root-finder adequacy, the executable decoder returns +exactly the low-degree roots of the interpolated polynomial that also pass the +agreement threshold. This is the strongest correctness statement available for +the current decoder without replacing the heuristic root finder or proving the +separate algebraic GS theorem that an actually-close message is a root of `Q`. -/ +theorem gsListDecode_exact_of_rootFinderAdequate + (code : ReedSolomonCode F) (received : Array F) (f : UniPoly F) + (hsize : received.size = code.n) + (m d : Nat) + (hparams : chooseParameters code.n code.k = (m, d)) + (hrootAdequate : + let q := interpolateWithMultiplicity code.domain received m d code.k + RootFinderAdequateFor q code.k received code.domain) : + f ∈ (gsListDecode code received).candidates ↔ + IsAcceptedCandidate code received f ∧ + (let q := interpolateWithMultiplicity code.domain received m d code.k + UniPoly.isZero (Bivariate.evaluateYPolynomial q f)) := by + simp only [gsListDecode] + have hneq : (received.size != code.n) = false := by simp [hsize] + simp [hneq, hparams, IsAcceptedCandidate] + constructor + · intro hmem + exact ⟨⟨hmem.2.1, hmem.2.2⟩, + (hrootAdequate.1 f hmem.1).2⟩ + · intro h + exact ⟨hrootAdequate.2 f h.1.1 h.2, h.1.1, + h.1.2⟩ + +end Generic +end GuruswamiSudan +end CodingTheory +end CompPoly diff --git a/CompPoly/CodingTheory/GuruswamiSudan/Counterexamples.lean b/CompPoly/CodingTheory/GuruswamiSudan/Counterexamples.lean new file mode 100644 index 00000000..9807c95f --- /dev/null +++ b/CompPoly/CodingTheory/GuruswamiSudan/Counterexamples.lean @@ -0,0 +1,55 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: OpenAI Codex +-/ + +import CompPoly.CodingTheory.GuruswamiSudan.Correctness +import CompPoly.Fields.KoalaBear + +/-! +# Counterexamples for the Heuristic Guruswami-Sudan Root Finder + +This file machine-checks a small obstruction to proving full completeness for +the current Lambdaworks-matching root finder: the root finder is heuristic, so it +does not enumerate every polynomial root of a bivariate polynomial. +-/ + +namespace CompPoly +namespace CodingTheory +namespace GuruswamiSudan +namespace Generic +namespace Counterexamples + +abbrev F := KoalaBear.Field + +def heuristicRootCounterexampleQ : Bivariate F := + Bivariate.fromMonomials + #[(0, 0), (0, 1), (0, 2)] + #[(5000 : F) * (5001 : F), -((5000 : F) + (5001 : F)), (1 : F)] + +def heuristicRootCounterexampleF : UniPoly F := #[(5000 : F)] + +theorem heuristicRootCounterexample_is_root : + UniPoly.isZero + (Bivariate.evaluateYPolynomial + heuristicRootCounterexampleQ heuristicRootCounterexampleF) := by + native_decide + +theorem heuristicRootCounterexample_not_returned : + heuristicRootCounterexampleF ∉ + findPolynomialRootsWithDomain heuristicRootCounterexampleQ 1 #[] #[] := by + native_decide + +theorem heuristicRootFinder_incomplete : + ¬ RootFinderCompleteFor heuristicRootCounterexampleQ 1 #[] #[] := by + intro hcomplete + exact heuristicRootCounterexample_not_returned + (hcomplete heuristicRootCounterexampleF (by native_decide) + heuristicRootCounterexample_is_root) + +end Counterexamples +end Generic +end GuruswamiSudan +end CodingTheory +end CompPoly diff --git a/CompPoly/CodingTheory/GuruswamiSudan/Generic.lean b/CompPoly/CodingTheory/GuruswamiSudan/Generic.lean new file mode 100644 index 00000000..d4949f32 --- /dev/null +++ b/CompPoly/CodingTheory/GuruswamiSudan/Generic.lean @@ -0,0 +1,823 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: OpenAI Codex +-/ + +import CompPoly.Fields.Basic + +/-! +# Generic Executable Guruswami-Sudan Decoder + +This module is a proof-light executable port of the Lambdaworks +`examples/reed-solomon-codes` Guruswami-Sudan decoder. It intentionally mirrors +the Lambdaworks educational implementation, including the same parameter search, +kernel-vector interpolation, Roth-Ruckenstein root-search heuristics, and small +brute-force fallbacks. It is parameterized over the executable field +representation so canonical and native-word KoalaBear backends share one +implementation. +-/ + +namespace CompPoly +namespace CodingTheory +namespace GuruswamiSudan +namespace Generic + +variable (F : Type) + +abbrev UniPoly := Array F + +structure Bivariate where + coeffs : Array (UniPoly F) +deriving BEq, Inhabited + +structure ReedSolomonCode where + n : Nat + k : Nat + domain : Array F +deriving BEq, Inhabited + +structure DecodeResult where + candidates : Array (UniPoly F) + multiplicity : Nat + degreeBound : Nat + errorBound : Nat +deriving BEq, Inhabited + +def maxCandidatesPerDepth : Nat := 15 +def maxTotalRoots : Nat := 10 + +namespace NatUtil + +def floorSqrt (n : Nat) : Nat := Id.run do + let mut r := 0 + while (r + 1) * (r + 1) <= n do + r := r + 1 + pure r + +def ceilSqrt (n : Nat) : Nat := + let r := floorSqrt n + if r * r = n then r else r + 1 + +def ceilSqrtDiv (num den : Nat) : Nat := Id.run do + if den = 0 then + pure 0 + else + let mut r := 0 + while r * r * den < num do + r := r + 1 + pure r + +def binomial (n k : Nat) : Nat := Id.run do + if k > n then + pure 0 + else if k = 0 || k = n then + pure 1 + else + let kk := min k (n - k) + let mut result := 1 + for i in [0:kk] do + result := result * (n - i) / (i + 1) + pure result + +end NatUtil + +variable {F : Type} [Field F] [BEq F] [Inhabited F] + +namespace UniPoly + +def zero : UniPoly F := #[0] + +def one : UniPoly F := #[1] + +def isZero (p : UniPoly F) : Bool := + p.all (fun c => c == 0) + +def trim (p : UniPoly F) : UniPoly F := Id.run do + let mut last := p.size + while 1 < last && p.getD (last - 1) 0 == 0 do + last := last - 1 + if last = 0 then + pure zero + else + pure (p.extract 0 last) + +def coeff (p : UniPoly F) (i : Nat) : F := + p.getD i 0 + +def degree (p : UniPoly F) : Nat := + let q := trim p + if q.isEmpty || isZero q then 0 else q.size - 1 + +def leadingCoeff (p : UniPoly F) : F := + let q := trim p + q.getD (q.size - 1) 0 + +def ofNatArray (xs : Array Nat) : UniPoly F := + let ys : Array F := xs.map (fun x : Nat => (x : F)) + trim ys + +def neg (p : UniPoly F) : UniPoly F := + trim (p.map fun c => -c) + +def add (p q : UniPoly F) : UniPoly F := Id.run do + let n := max p.size q.size + let mut out := Array.replicate n (0 : F) + for i in [0:n] do + out := out.set! i (p.getD i 0 + q.getD i 0) + pure (trim out) + +def sub (p q : UniPoly F) : UniPoly F := + add p (neg q) + +def scale (a : F) (p : UniPoly F) : UniPoly F := + trim (p.map fun c => a * c) + +def mul (p q : UniPoly F) : UniPoly F := Id.run do + if isZero p || isZero q then + pure zero + else + let n := p.size + q.size - 1 + let mut out := Array.replicate n (0 : F) + for i in [0:p.size] do + for j in [0:q.size] do + let value := out.getD (i + j) 0 + p.getD i 0 * q.getD j 0 + out := out.set! (i + j) value + pure (trim out) + +def mulXPow (p : UniPoly F) (k : Nat) : UniPoly F := + if isZero p then + zero + else + trim (Array.replicate k (0 : F) ++ p) + +def evaluate (p : UniPoly F) (x : F) : F := + p.foldr (fun c acc => acc * x + c) 0 + +def pow (p : UniPoly F) (e : Nat) : UniPoly F := Id.run do + let mut result := one + for _ in [0:e] do + result := mul result p + pure result + +def monomial (i : Nat) (c : F) : UniPoly F := + if c == 0 then + zero + else + (Array.replicate i (0 : F)).push c + +def longDivRem? (num den : UniPoly F) : Option (Prod (UniPoly F) (UniPoly F)) := Id.run do + let den := trim den + if isZero den then + pure none + else + let mut rem := trim num + let denDeg := degree den + let denLead := leadingCoeff den + let mut quot := Array.replicate + (if degree rem < denDeg then 1 else degree rem - denDeg + 1) (0 : F) + while !isZero rem && denDeg <= degree rem do + let shift := degree rem - denDeg + let coeff := leadingCoeff rem / denLead + if quot.size <= shift then + quot := quot ++ Array.replicate (shift + 1 - quot.size) (0 : F) + quot := quot.set! shift (quot.getD shift 0 + coeff) + rem := sub rem (mulXPow (scale coeff den) shift) + pure (some (trim quot, trim rem)) + +def coeffsEq (p q : UniPoly F) : Bool := + trim p == trim q + +end UniPoly + +namespace Bivariate + +def zero : Bivariate F := { coeffs := #[UniPoly.zero] } + +def trimCoeffs (coeffs : Array (UniPoly F)) : Array (UniPoly F) := Id.run do + let mut last := coeffs.size + while 1 < last && UniPoly.isZero (coeffs.getD (last - 1) UniPoly.zero) do + last := last - 1 + if last = 0 then + pure #[UniPoly.zero] + else + pure (coeffs.extract 0 last) + +def ofCoeffs (coeffs : Array (UniPoly F)) : Bivariate F := + { coeffs := trimCoeffs coeffs } + +def isZero (q : Bivariate F) : Bool := + q.coeffs.size = 1 && UniPoly.isZero (q.coeffs.getD 0 UniPoly.zero) + +def yDegree (q : Bivariate F) : Nat := + if isZero q then 0 else q.coeffs.size - 1 + +def maxXDegree (q : Bivariate F) : Nat := + q.coeffs.foldl (fun acc p => max acc (UniPoly.degree p)) 0 + +def coeff (q : Bivariate F) (i j : Nat) : F := + UniPoly.coeff (q.coeffs.getD j UniPoly.zero) i + +def evaluate (q : Bivariate F) (x y : F) : F := Id.run do + let mut result : F := 0 + let mut yPow := (1 : F) + for j in [0:q.coeffs.size] do + let xEval := UniPoly.evaluate (q.coeffs.getD j UniPoly.zero) x + result := result + xEval * yPow + yPow := yPow * y + pure result + +def weightedDegree (q : Bivariate F) (w : Nat) : Nat := Id.run do + let mut maxDeg := 0 + for j in [0:q.coeffs.size] do + let p := q.coeffs.getD j UniPoly.zero + for i in [0:p.size] do + if p.getD i 0 != 0 then + maxDeg := max maxDeg (i + w * j) + pure maxDeg + +def evaluateYPolynomial (q : Bivariate F) (f : UniPoly F) : UniPoly F := Id.run do + let mut result := UniPoly.zero + let mut fPow := UniPoly.one + for j in [0:q.coeffs.size] do + let term := UniPoly.mul (q.coeffs.getD j UniPoly.zero) fPow + result := UniPoly.add result term + fPow := UniPoly.mul fPow f + pure (UniPoly.trim result) + +def fromMonomials (monomials : Array (Prod Nat Nat)) (coefficients : Array F) : Bivariate F := Id.run do + let maxJ := monomials.foldl (fun acc pair => max acc pair.2) 0 + let mut rows : Array (UniPoly F) := Array.replicate (maxJ + 1) (#[] : UniPoly F) + for idx in [0:monomials.size] do + let (i, j) := monomials.getD idx (0, 0) + let c := coefficients.getD idx 0 + let mut row := rows.getD j #[] + while row.size <= i do + row := row.push 0 + row := row.set! i c + rows := rows.set! j row + let polys := rows.map fun row => + if row.isEmpty then UniPoly.zero else UniPoly.trim row + pure (ofCoeffs polys) + +end Bivariate + +namespace ReedSolomonCode + +def withDomain (domain : Array F) (k : Nat) : ReedSolomonCode F := + { n := domain.size, k := k, domain := domain } + +def consecutiveDomain (n : Nat) : Array F := Id.run do + let mut out := #[] + for i in [0:n] do + out := out.push (i : F) + pure out + +def withConsecutiveDomain (n k : Nat) : ReedSolomonCode F := + withDomain (consecutiveDomain n) k + +def pow2Domain (logN : Nat) (omega : F) : Array F := Id.run do + let n := 2 ^ logN + let mut out := #[] + let mut cur := (1 : F) + for _ in [0:n] do + out := out.push cur + cur := cur * omega + pure out + +def withRootsOfUnityDomain (logN k : Nat) (omega : F) : ReedSolomonCode F := + withDomain (pow2Domain logN omega) k + +def encodePolynomial (code : ReedSolomonCode F) (poly : UniPoly F) : Array F := + code.domain.map fun x => UniPoly.evaluate poly x + +def encode (code : ReedSolomonCode F) (message : Array F) : Array F := + encodePolynomial code (UniPoly.trim message) + +end ReedSolomonCode + +def gsDecodingRadius (n k : Nat) : Nat := + let s := NatUtil.floorSqrt (n * k) + if n <= s then 0 else n - s - 1 + +def johnsonListBound? (n k t : Nat) : Option Float := + let s := NatUtil.floorSqrt (n * k) + let denominatorInt := n - t + if denominatorInt <= s then + none + else + some ((Float.ofNat n) / (Float.ofNat (denominatorInt - s))) + +def countMonomials (d w : Nat) : Nat := Id.run do + if w = 0 then + pure d + else + let maxJ := d / w + 1 + let mut count := 0 + for j in [0:maxJ + 1] do + let maxI := d - w * j + count := count + maxI + pure count + +def chooseParameters (n k : Nat) : Prod Nat Nat := Id.run do + let targetRadius := gsDecodingRadius n k + for m in [1:21] do + let radiusWithM := n - NatUtil.ceilSqrtDiv (n * k * (m + 1)) m + if radiusWithM >= targetRadius then + let constraintsPerPoint := m * (m + 1) / 2 + let totalConstraints := n * constraintsPerPoint + let d := NatUtil.ceilSqrt (2 * totalConstraints * (k - 1)) + k + let numMonomials := countMonomials d (k - 1) + if numMonomials > totalConstraints then + return (m, d) + for m in [2:21] do + let constraintsPerPoint := m * (m + 1) / 2 + let totalConstraints := n * constraintsPerPoint + let d := NatUtil.ceilSqrt (2 * totalConstraints * (k - 1)) + k + let numMonomials := countMonomials d (k - 1) + if numMonomials > totalConstraints then + return (m, d) + pure (4, n + k) + +def monomialsBelowWeightedDegree (d k : Nat) : Array (Prod Nat Nat) := Id.run do + let w := k - 1 + let maxJ := if w = 0 then 0 else d / w + 1 + let mut monomials := #[] + for j in [0:maxJ + 1] do + let maxI := d - w * j + for i in [0:maxI] do + monomials := monomials.push (i, j) + pure monomials + +def findKernelVector (matrix : Array (Array F)) (numCols : Nat) : Array F := Id.run do + let m := matrix.size + let n := numCols + if m = 0 then + let mut result := Array.replicate n (0 : F) + if 0 < n then + result := result.set! 0 1 + pure result + else + let mut mat := matrix.map fun row => + if row.size < n then row ++ Array.replicate (n - row.size) (0 : F) else row + let mut pivotCols : Array Nat := #[] + let mut pivotRow := 0 + for col in [0:n] do + if pivotRow < m then + let mut found : Option Nat := none + for row in [pivotRow:m] do + if found.isNone && mat[row]!.getD col 0 != 0 then + found := some row + match found with + | none => pure () + | some row => + let pivotData := mat[pivotRow]! + let rowData := mat[row]! + mat := (mat.set! pivotRow rowData).set! row pivotData + pivotCols := pivotCols.push col + let pivot := mat[pivotRow]!.getD col 0 + let pivotInv := pivot⁻¹ + let mut pivotRowData := mat[pivotRow]! + for j in [col:n] do + pivotRowData := pivotRowData.set! j (pivotRowData.getD j 0 * pivotInv) + mat := mat.set! pivotRow pivotRowData + let pivotSlice := mat[pivotRow]!.extract col n + for row2 in [0:m] do + if row2 != pivotRow && mat[row2]!.getD col 0 != 0 then + let factor := mat[row2]!.getD col 0 + let mut rowData := mat[row2]! + for offset in [0:pivotSlice.size] do + let j := col + offset + let sub := factor * pivotSlice.getD offset 0 + rowData := rowData.set! j (rowData.getD j 0 - sub) + mat := mat.set! row2 rowData + pivotRow := pivotRow + 1 + let mut freeCol : Option Nat := none + for col in [0:n] do + if freeCol.isNone && !(pivotCols.contains col) then + freeCol := some col + let mut kernel := Array.replicate n (0 : F) + match freeCol with + | some fc => + kernel := kernel.set! fc 1 + for row in [0:pivotCols.size] do + let pc := pivotCols[row]! + if row < m then + kernel := kernel.set! pc (-(mat[row]!.getD fc 0)) + | none => + if 0 < n then + kernel := kernel.set! (n - 1) 1 + pure kernel + +def interpolationMatrix + (domain received : Array F) (m d k : Nat) : Array (Array F) := Id.run do + let monomials := monomialsBelowWeightedDegree d k + let numMonomials := monomials.size + let constraintsPerPoint := m * (m + 1) / 2 + let totalConstraints := domain.size * constraintsPerPoint + let mut matrix : Array (Array F) := Array.mkEmpty totalConstraints + for idx in [0:domain.size] do + let alpha := domain.getD idx 0 + let y := received.getD idx 0 + for totalOrder in [0:m] do + for b in [0:totalOrder + 1] do + let a := totalOrder - b + let mut row := Array.replicate numMonomials (0 : F) + for monIdx in [0:monomials.size] do + let (i, j) := monomials[monIdx]! + if !(i < a || j < b) then + let coeffScalar := NatUtil.binomial i a * NatUtil.binomial j b + let coeff : F := (coeffScalar : F) * (alpha ^ (i - a)) * (y ^ (j - b)) + row := row.set! monIdx coeff + matrix := matrix.push row + pure matrix + +def interpolateWithMultiplicity + (domain received : Array F) (m d k : Nat) : Bivariate F := + let monomials := monomialsBelowWeightedDegree d k + let numMonomials := monomials.size + let matrix := interpolationMatrix domain received m d k + let solution := findKernelVector matrix numMonomials + Bivariate.fromMonomials monomials solution + +def extractSmallValue? (fe : F) : Option Nat := Id.run do + for i in [0:101] do + if fe == (i : F) then + return some i + pure none + +def lagrangeInterpolateAtZeroWithPoints? (points : Array (Prod Nat F)) : Option F := Id.run do + if points.isEmpty then + pure none + else + let n := points.size + let mut result : F := 0 + for i in [0:n] do + let (xi, yi) := points[i]! + let xiFe : F := xi + let mut numerator := (1 : F) + let mut denominator := (1 : F) + for j in [0:n] do + if j != i then + let (xj, _) := points[j]! + let xjFe : F := xj + numerator := numerator * (-xjFe) + denominator := denominator * (xiFe - xjFe) + if denominator == 0 then + return none + result := result + yi * (numerator / denominator) + pure (some result) + +def lagrangeInterpolatePolynomial? (points : Array (Prod Nat F)) (maxDegree : Nat) : + Option (UniPoly F) := Id.run do + if points.isEmpty || points.size > maxDegree then + pure none + else + let n := points.size + let mut coeffs := Array.replicate n (0 : F) + for i in [0:n] do + let (xi, yi) := points[i]! + let xiFe : F := xi + let mut basis : UniPoly F := #[1] + let mut denominator := (1 : F) + for j in [0:n] do + if j != i then + let (xj, _) := points[j]! + let xjFe : F := xj + let mut next := Array.replicate (basis.size + 1) (0 : F) + for bIdx in [0:basis.size] do + let c := basis[bIdx]! + next := next.set! (bIdx + 1) (next.getD (bIdx + 1) 0 + c) + next := next.set! bIdx (next.getD bIdx 0 - c * xjFe) + basis := next + denominator := denominator * (xiFe - xjFe) + if denominator == 0 then + return none + let denomInv := 1 / denominator + for bIdx in [0:basis.size] do + if bIdx < coeffs.size then + coeffs := coeffs.set! bIdx (coeffs.getD bIdx 0 + yi * (basis[bIdx]! * denomInv)) + pure (some (UniPoly.trim coeffs)) + +def substituteAndDivide (q : Bivariate F) (c : F) : Bivariate F := Id.run do + let yDeg := Bivariate.yDegree q + let maxXDeg := Bivariate.maxXDegree q + yDeg + let maxYDeg := yDeg + let mut result : Array (UniPoly F) := + Array.replicate (maxYDeg + 1) (Array.replicate (maxXDeg + 2) (0 : F)) + for j in [0:q.coeffs.size] do + let qj := q.coeffs[j]! + for kk in [0:j + 1] do + let binom := NatUtil.binomial j kk + let scale : F := (binom : F) * (c ^ (j - kk)) + let mut row := result[kk]! + for i in [0:qj.size] do + let xPower := i + kk + if xPower <= maxXDeg + 1 && kk <= maxYDeg then + row := row.set! xPower (row.getD xPower 0 + qj.getD i 0 * scale) + result := result.set! kk row + let divided := result.map fun row => + if row.size <= 1 then UniPoly.zero else UniPoly.trim (row.extract 1 row.size) + pure (Bivariate.ofCoeffs divided) + +def findRootsLinearY (q : Bivariate F) (maxDegree : Nat) : Array (UniPoly F) := + if q.coeffs.size < 2 then + #[] + else + let a := q.coeffs.getD 0 UniPoly.zero + let b := q.coeffs.getD 1 UniPoly.zero + let negA := UniPoly.neg a + match UniPoly.longDivRem? negA b with + | none => #[] + | some (quot, rem) => + if !UniPoly.isZero rem then + #[] + else if UniPoly.degree quot >= maxDegree then + #[] + else + #[quot] + +def appendUniquePoly (roots : Array (UniPoly F)) (candidate : UniPoly F) : Array (UniPoly F) := + if roots.any (fun p => UniPoly.coeffsEq p candidate) then roots else roots.push candidate + +def appendUniqueField (roots : Array F) (candidate : F) : Array F := + if roots.contains candidate then roots else roots.push candidate + +partial def enumerateSmallPolys + (q : Bivariate F) (maxDegree adjustedMaxCoeff totalCandidates idx : Nat) + (roots : Array (UniPoly F)) : Array (UniPoly F) := + if idx >= totalCandidates then + roots + else + Id.run do + let mut coeffs : Array F := #[] + let mut value := idx + for _ in [0:maxDegree] do + coeffs := coeffs.push ((value % (adjustedMaxCoeff + 1) : Nat) : F) + value := value / (adjustedMaxCoeff + 1) + let candidate := UniPoly.trim coeffs + let roots := + if UniPoly.isZero candidate then + roots + else if UniPoly.isZero (Bivariate.evaluateYPolynomial q candidate) then + appendUniquePoly roots candidate + else + roots + pure (enumerateSmallPolys q maxDegree adjustedMaxCoeff totalCandidates (idx + 1) roots) + +def trySmallIntegerPolynomials + (q : Bivariate F) (maxDegree maxCoeff : Nat) (roots : Array (UniPoly F)) : Array (UniPoly F) := + let adjustedMaxCoeff := + if maxDegree <= 2 then maxCoeff + else if maxDegree <= 3 then min maxCoeff 30 + else if maxDegree <= 4 then min maxCoeff 15 + else min maxCoeff 8 + let totalCandidates := (adjustedMaxCoeff + 1) ^ maxDegree + if totalCandidates > 200000 then + roots + else + enumerateSmallPolys q maxDegree adjustedMaxCoeff totalCandidates 0 roots + +def tryDirectRoots + (q : Bivariate F) (maxDegree : Nat) (hintValues : Array F) + (roots : Array (UniPoly F)) : Array (UniPoly F) := Id.run do + let mut roots := roots + if maxDegree <= 4 then + roots := trySmallIntegerPolynomials q maxDegree 20 roots + for hint in hintValues do + let candidate := #[hint] + if UniPoly.isZero (Bivariate.evaluateYPolynomial q candidate) then + roots := appendUniquePoly roots candidate + pure roots + +def tryInterpolatedCandidates + (q : Bivariate F) (maxDegree : Nat) (hintValues domain : Array F) + (roots : Array (UniPoly F)) : Array (UniPoly F) := Id.run do + let mut roots := roots + if hintValues.size < maxDegree || domain.size < maxDegree then + pure roots + else + let n := min hintValues.size domain.size + let limit := min (n - maxDegree + 1) 20 + for start in [0:limit] do + if start + maxDegree <= n then + let mut points : Array (Prod Nat F) := #[] + for idx in [start:start + maxDegree] do + match extractSmallValue? (domain[idx]!) with + | some alpha => points := points.push (alpha, hintValues[idx]!) + | none => pure () + if points.size = maxDegree then + match lagrangeInterpolatePolynomial? points maxDegree with + | some candidate => + if UniPoly.isZero (Bivariate.evaluateYPolynomial q candidate) then + roots := appendUniquePoly roots candidate + if roots.size >= maxTotalRoots then + return roots + | none => pure () + pure roots + +partial def findUnivariateRootsWithHints (coeffs hintValues : Array F) : Array F := + if coeffs.isEmpty || coeffs.all (fun c => c == 0) then + Id.run do + let mut roots : Array F := #[] + for i in [0:20] do + roots := appendUniqueField roots (i : F) + for hint in hintValues do + roots := appendUniqueField roots hint + pure roots + else + let poly := UniPoly.trim coeffs + let deg := UniPoly.degree poly + if deg = 0 then + #[] + else if deg = 1 then + let a := poly.getD 1 0 + let b := poly.getD 0 0 + if a == 0 then #[] else #[(-b) / a] + else + Id.run do + let maxRoots := deg + let mut roots : Array F := #[] + for hint in hintValues do + if UniPoly.evaluate poly hint == 0 then + roots := appendUniqueField roots hint + if roots.size >= maxRoots then + return roots + for i in [0:2000] do + let elem : F := i + if UniPoly.evaluate poly elem == 0 then + roots := appendUniqueField roots elem + if roots.size >= maxRoots then + return roots + for i in [1:2000] do + let elem : F := -(i : F) + if UniPoly.evaluate poly elem == 0 then + roots := appendUniqueField roots elem + if roots.size >= maxRoots then + return roots + pure roots + +partial def findUnivariateRootsWithHintsAndDomain + (coeffs hintValues domain : Array F) : Array F := + if coeffs.isEmpty || coeffs.all (fun c => c == 0) then + Id.run do + let mut roots : Array F := #[] + if hintValues.size >= 3 && domain.size >= 3 then + let minLen := min hintValues.size domain.size + for start in [0:min minLen 10] do + for size in [3:min minLen 6 + 1] do + if start + size <= minLen then + let mut subset : Array (Prod Nat F) := #[] + for idx in [start:start + size] do + match extractSmallValue? (domain[idx]!) with + | some alpha => subset := subset.push (alpha, hintValues[idx]!) + | none => pure () + if subset.size = size then + match lagrangeInterpolateAtZeroWithPoints? subset with + | some root => roots := appendUniqueField roots root + | none => pure () + for i in [0:20] do + roots := appendUniqueField roots (i : F) + pure roots + else + findUnivariateRootsWithHints coeffs hintValues +partial def rrSearchWithDomain + (q : Bivariate F) (maxDegree : Nat) (currentCoeffs : Array F) + (roots : Array (UniPoly F)) (hintValues domain : Array F) (depth : Nat) : Array (UniPoly F) := + if roots.size >= maxTotalRoots then + roots + else + Id.run do + let qAtZero := q.coeffs.map fun p => UniPoly.evaluate p 0 + let yRootsAll := findUnivariateRootsWithHintsAndDomain qAtZero hintValues domain + let yRoots := yRootsAll.extract 0 (min yRootsAll.size maxCandidatesPerDepth) + let mut roots := roots + for yRoot in yRoots do + if roots.size >= maxTotalRoots then + return roots + let newCoeffs := currentCoeffs.push yRoot + if newCoeffs.size <= maxDegree then + let qTransformed := substituteAndDivide q yRoot + if Bivariate.isZero qTransformed then + roots := appendUniquePoly roots (UniPoly.trim newCoeffs) + else if newCoeffs.size < maxDegree then + let mut transformedHints : Array F := #[] + let mut filteredDomain : Array F := #[] + for idx in [0:min hintValues.size domain.size] do + let alpha := domain[idx]! + if alpha != 0 then + transformedHints := transformedHints.push ((hintValues[idx]! - yRoot) / alpha) + filteredDomain := filteredDomain.push alpha + roots := rrSearchWithDomain qTransformed maxDegree newCoeffs roots + transformedHints filteredDomain (depth + 1) + let candidate := UniPoly.trim newCoeffs + if !candidate.isEmpty && + UniPoly.isZero (Bivariate.evaluateYPolynomial q candidate) then + roots := appendUniquePoly roots candidate + pure roots + +def findPolynomialRootsWithDomain + (q : Bivariate F) (maxDegree : Nat) (hintValues domain : Array F) : Array (UniPoly F) := + /- + This root finder intentionally mirrors the Lambdaworks educational heuristic: + it tries a linear special case, interpolated candidates from hints, bounded + small-coefficient search, and a Roth-Ruckenstein-style search with capped + candidates. + + It is not complete as a polynomial root finder. For example over KoalaBear, + the bivariate polynomial `Q(Y) = (Y - 5000) * (Y - 5001)` has the constant + polynomial `f(X) = 5000` as a root, but with empty hints + + `findPolynomialRootsWithDomain Q 1 #[] #[] = #[]`. + + A full GS completeness theorem therefore has to assume completeness of this + root-finder call, or use a separate verified complete root finder. + -/ + if Bivariate.yDegree q <= 1 then + findRootsLinearY q maxDegree + else + let roots := tryInterpolatedCandidates q maxDegree hintValues domain #[] + let roots := + if maxDegree <= 10 && roots.size < maxTotalRoots then + tryDirectRoots q maxDegree hintValues roots + else + roots + if roots.size < maxTotalRoots then + rrSearchWithDomain q maxDegree #[] roots hintValues domain 0 + else + roots + +def agreement (received domain : Array F) (poly : UniPoly F) : Nat := Id.run do + let mut count := 0 + for i in [0:min received.size domain.size] do + if UniPoly.evaluate poly domain[i]! == received[i]! then + count := count + 1 + pure count + +def gsListDecodeWithMultiplicity + (code : ReedSolomonCode F) (received : Array F) (multiplicity : Nat) : DecodeResult F := + let n := code.n + let k := code.k + let m := multiplicity + let constraintsPerPoint := m * (m + 1) / 2 + let totalConstraints := n * constraintsPerPoint + let d := NatUtil.ceilSqrt (2 * totalConstraints * (k - 1)) + k + let q := interpolateWithMultiplicity code.domain received m d k + let allRoots := findPolynomialRootsWithDomain q k received code.domain + let errorBound := gsDecodingRadius n k + let agreementThreshold := n - errorBound + let candidates := + allRoots.filter fun f => + UniPoly.degree f < k && agreement received code.domain f >= agreementThreshold + { + candidates := candidates + multiplicity := m + degreeBound := d + errorBound := errorBound + } + +def gsListDecode (code : ReedSolomonCode F) (received : Array F) : DecodeResult F := + if received.size != code.n then + panic! "received word length must equal code length" + else + let params := chooseParameters code.n code.k + let m := params.1 + let d := params.2 + let q := interpolateWithMultiplicity code.domain received m d code.k + let allRoots := findPolynomialRootsWithDomain q code.k received code.domain + let errorBound := gsDecodingRadius code.n code.k + let agreementThreshold := code.n - errorBound + let candidates := + allRoots.filter fun f => + UniPoly.degree f < code.k && agreement received code.domain f >= agreementThreshold + { + candidates := candidates + multiplicity := m + degreeBound := d + errorBound := errorBound + } + +def introduceErrors (codeword : Array F) (positions values : Array Nat) : Array F := Id.run do + let mut out := codeword + for i in [0:min positions.size values.size] do + let pos := positions[i]! + if pos < out.size then + out := out.set! pos (out[pos]! + (values[i]! : F)) + pure out + +def introduceErrorsAtPositions (codeword : Array F) (positions : Array Nat) : Array F := + introduceErrors codeword positions (positions.mapIdx fun i _ => i + 1) + +def polyToNats (toNat : F → Nat) (p : UniPoly F) : Array Nat := + (UniPoly.trim p).map toNat + +def valuesToNats (toNat : F → Nat) (xs : Array F) : Array Nat := + xs.map toNat + +def candidatesToNats (toNat : F → Nat) (xs : Array (UniPoly F)) : Array (Array Nat) := + xs.map (polyToNats toNat) + +end Generic +end GuruswamiSudan +end CodingTheory +end CompPoly diff --git a/CompPoly/CodingTheory/GuruswamiSudanFast.lean b/CompPoly/CodingTheory/GuruswamiSudanFast.lean new file mode 100644 index 00000000..b978a125 --- /dev/null +++ b/CompPoly/CodingTheory/GuruswamiSudanFast.lean @@ -0,0 +1,106 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: OpenAI Codex +-/ + +import CompPoly.CodingTheory.GuruswamiSudan.Generic +import CompPoly.Fields.KoalaBear + +/-! +# Executable Guruswami-Sudan Decoder (Fast KoalaBear) + +Native-word `KoalaBear.Fast.Field` instantiation of the generic executable +Guruswami-Sudan decoder. +-/ + +namespace CompPoly +namespace CodingTheory +namespace GuruswamiSudanFast + +abbrev F := KoalaBear.Fast.Field +abbrev UniPoly := CompPoly.CodingTheory.GuruswamiSudan.Generic.UniPoly F +abbrev Bivariate := CompPoly.CodingTheory.GuruswamiSudan.Generic.Bivariate F +abbrev ReedSolomonCode := CompPoly.CodingTheory.GuruswamiSudan.Generic.ReedSolomonCode F +abbrev DecodeResult := CompPoly.CodingTheory.GuruswamiSudan.Generic.DecodeResult F + +instance : Inhabited F := ⟨0⟩ + +def maxCandidatesPerDepth : Nat := CompPoly.CodingTheory.GuruswamiSudan.Generic.maxCandidatesPerDepth +def maxTotalRoots : Nat := CompPoly.CodingTheory.GuruswamiSudan.Generic.maxTotalRoots + +namespace NatUtil + +export CompPoly.CodingTheory.GuruswamiSudan.Generic.NatUtil (floorSqrt ceilSqrt ceilSqrtDiv binomial) + +end NatUtil + +namespace UniPoly + +export CompPoly.CodingTheory.GuruswamiSudan.Generic.UniPoly + (zero one isZero trim coeff degree leadingCoeff ofNatArray neg add sub scale mul + mulXPow evaluate pow monomial longDivRem? coeffsEq) + +end UniPoly + +namespace Bivariate + +export CompPoly.CodingTheory.GuruswamiSudan.Generic.Bivariate + (zero trimCoeffs ofCoeffs isZero yDegree maxXDegree coeff evaluate weightedDegree + evaluateYPolynomial fromMonomials) + +end Bivariate + +namespace ReedSolomonCode + +def withDomain (domain : Array F) (k : Nat) : ReedSolomonCode := + CompPoly.CodingTheory.GuruswamiSudan.Generic.ReedSolomonCode.withDomain domain k + +def consecutiveDomain (n : Nat) : Array F := + CompPoly.CodingTheory.GuruswamiSudan.Generic.ReedSolomonCode.consecutiveDomain (F := F) n + +def withConsecutiveDomain (n k : Nat) : ReedSolomonCode := + CompPoly.CodingTheory.GuruswamiSudan.Generic.ReedSolomonCode.withConsecutiveDomain (F := F) n k + +def pow2Domain (logN : Nat) : Array F := + let omega := KoalaBear.Fast.ofField + (KoalaBear.twoAdicGenerators.toArray.getD logN (1 : KoalaBear.Field)) + CompPoly.CodingTheory.GuruswamiSudan.Generic.ReedSolomonCode.pow2Domain (F := F) logN omega + +def withRootsOfUnityDomain (logN k : Nat) : ReedSolomonCode := + let omega := KoalaBear.Fast.ofField + (KoalaBear.twoAdicGenerators.toArray.getD logN (1 : KoalaBear.Field)) + CompPoly.CodingTheory.GuruswamiSudan.Generic.ReedSolomonCode.withRootsOfUnityDomain (F := F) logN k omega + +def encodePolynomial (code : ReedSolomonCode) (poly : UniPoly) : Array F := + CompPoly.CodingTheory.GuruswamiSudan.Generic.ReedSolomonCode.encodePolynomial code poly + +def encode (code : ReedSolomonCode) (message : Array F) : Array F := + CompPoly.CodingTheory.GuruswamiSudan.Generic.ReedSolomonCode.encode code message + +end ReedSolomonCode + +export CompPoly.CodingTheory.GuruswamiSudan.Generic + (gsDecodingRadius johnsonListBound? countMonomials chooseParameters + monomialsBelowWeightedDegree findKernelVector interpolateWithMultiplicity + extractSmallValue? lagrangeInterpolateAtZeroWithPoints? + lagrangeInterpolatePolynomial? substituteAndDivide findRootsLinearY + appendUniquePoly appendUniqueField enumerateSmallPolys + trySmallIntegerPolynomials tryDirectRoots tryInterpolatedCandidates + findUnivariateRootsWithHints findUnivariateRootsWithHintsAndDomain + rrSearchWithDomain findPolynomialRootsWithDomain agreement + gsListDecodeWithMultiplicity gsListDecode introduceErrors + introduceErrorsAtPositions) + +def polyToNats (p : UniPoly) : Array Nat := + CompPoly.CodingTheory.GuruswamiSudan.Generic.polyToNats (F := F) KoalaBear.Fast.toNat p + +def valuesToNats (xs : Array F) : Array Nat := + CompPoly.CodingTheory.GuruswamiSudan.Generic.valuesToNats (F := F) KoalaBear.Fast.toNat xs + +def candidatesToNats (xs : Array UniPoly) : Array (Array Nat) := + CompPoly.CodingTheory.GuruswamiSudan.Generic.candidatesToNats (F := F) KoalaBear.Fast.toNat xs + +end GuruswamiSudanFast +end CodingTheory +end CompPoly diff --git a/bench/CompPolyBench/CodingTheory.lean b/bench/CompPolyBench/CodingTheory.lean new file mode 100644 index 00000000..34f9aca8 --- /dev/null +++ b/bench/CompPolyBench/CodingTheory.lean @@ -0,0 +1,19 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: OpenAI Codex +-/ + +import CompPolyBench.CodingTheory.GuruswamiSudan + +/-! +# Coding-Theory Benchmarks +-/ + +namespace CompPolyBench + +/-- Metadata for all coding-theory benchmark modules. -/ +def codingTheoryGroupInfos : List BenchGroupInfo := + gsGroupInfos + +end CompPolyBench diff --git a/bench/CompPolyBench/CodingTheory/GuruswamiSudan.lean b/bench/CompPolyBench/CodingTheory/GuruswamiSudan.lean new file mode 100644 index 00000000..28133b16 --- /dev/null +++ b/bench/CompPolyBench/CodingTheory/GuruswamiSudan.lean @@ -0,0 +1,218 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: OpenAI Codex +-/ + +import CompPolyBench.Common +import CompPoly.CodingTheory.GuruswamiSudan +import CompPoly.CodingTheory.GuruswamiSudanFast + +/-! +# Guruswami-Sudan Decoder Benchmarks + +Benchmarks for the executable KoalaBear Guruswami-Sudan decoder. +-/ + +namespace CompPolyBench + +open CompPoly.CodingTheory.GuruswamiSudan + +private def gsWarmupIterations (preset : BenchPreset) : Nat := + preset.selectNat 2 1 0 + +private def gsMeasuredIterations (preset : BenchPreset) : Nat := + preset.selectNat 20 5 1 + +private def checksumUniPoly (p : UniPoly) : Nat := + checksumArray checksumKoalaBear p + +private def checksumDecodeResult (result : DecodeResult) : Nat := + let candidateChecksum := checksumArray checksumUniPoly result.candidates + mixChecksum + (mixChecksum + (mixChecksum candidateChecksum result.multiplicity) + result.degreeBound) + result.errorBound + +private def checksumFastUniPoly + (p : CompPoly.CodingTheory.GuruswamiSudanFast.UniPoly) : Nat := + checksumArray KoalaBear.Fast.toNat p + +private def checksumFastDecodeResult + (result : CompPoly.CodingTheory.GuruswamiSudanFast.DecodeResult) : Nat := + let candidateChecksum := checksumArray checksumFastUniPoly result.candidates + mixChecksum + (mixChecksum + (mixChecksum candidateChecksum result.multiplicity) + result.degreeBound) + result.errorBound + +private structure GsBenchInput where + code : ReedSolomonCode + received : Array F + inputShape : String + +private structure FastGsBenchInput where + code : CompPoly.CodingTheory.GuruswamiSudanFast.ReedSolomonCode + received : Array CompPoly.CodingTheory.GuruswamiSudanFast.F + inputShape : String + +private def consecutiveInput (withError : Bool) : GsBenchInput := + let code := ReedSolomonCode.withConsecutiveDomain 8 2 + let message : Array F := #[1, 2] + let codeword := ReedSolomonCode.encode code message + let received := + if withError then + introduceErrorsAtPositions codeword #[1] + else + codeword + { + code := code + received := received + inputShape := + if withError then "RS(8,2), consecutive domain, 1 error" + else "RS(8,2), consecutive domain, no errors" + } + +private def fastConsecutiveInput (withError : Bool) : FastGsBenchInput := + let code := CompPoly.CodingTheory.GuruswamiSudanFast.ReedSolomonCode.withConsecutiveDomain 8 2 + let message : Array CompPoly.CodingTheory.GuruswamiSudanFast.F := #[1, 2] + let codeword := CompPoly.CodingTheory.GuruswamiSudanFast.ReedSolomonCode.encode code message + let received := + if withError then + CompPoly.CodingTheory.GuruswamiSudanFast.introduceErrorsAtPositions codeword #[1] + else + codeword + { + code := code + received := received + inputShape := + if withError then "RS(8,2), consecutive domain, 1 error" + else "RS(8,2), consecutive domain, no errors" + } + +private def consecutiveRs16Input : GsBenchInput := + let code := ReedSolomonCode.withConsecutiveDomain 16 4 + let message : Array F := #[1, 2, 3, 4] + let codeword := ReedSolomonCode.encode code message + { + code := code + received := codeword + inputShape := "RS(16,4), consecutive domain, no errors" + } + +private def fastConsecutiveRs16Input : FastGsBenchInput := + let code := CompPoly.CodingTheory.GuruswamiSudanFast.ReedSolomonCode.withConsecutiveDomain 16 4 + let message : Array CompPoly.CodingTheory.GuruswamiSudanFast.F := #[1, 2, 3, 4] + let codeword := CompPoly.CodingTheory.GuruswamiSudanFast.ReedSolomonCode.encode code message + { + code := code + received := codeword + inputShape := "RS(16,4), consecutive domain, no errors" + } + +private def rootsOfUnityInput (withError : Bool) : GsBenchInput := + let code := ReedSolomonCode.withRootsOfUnityDomain 3 2 + let message : Array F := #[1, 2] + let codeword := ReedSolomonCode.encode code message + let received := + if withError then + introduceErrorsAtPositions codeword #[2] + else + codeword + { + code := code + received := received + inputShape := + if withError then "RS(8,2), roots-of-unity domain, 1 error" + else "RS(8,2), roots-of-unity domain, no errors" + } + +private def fastRootsOfUnityInput (withError : Bool) : FastGsBenchInput := + let code := CompPoly.CodingTheory.GuruswamiSudanFast.ReedSolomonCode.withRootsOfUnityDomain 3 2 + let message : Array CompPoly.CodingTheory.GuruswamiSudanFast.F := #[1, 2] + let codeword := CompPoly.CodingTheory.GuruswamiSudanFast.ReedSolomonCode.encode code message + let received := + if withError then + CompPoly.CodingTheory.GuruswamiSudanFast.introduceErrorsAtPositions codeword #[2] + else + codeword + { + code := code + received := received + inputShape := + if withError then "RS(8,2), roots-of-unity domain, 1 error" + else "RS(8,2), roots-of-unity domain, no errors" + } + +private def runGsDecodeGroup (info : BenchGroupInfo) (input : GsBenchInput) + (fastInput : FastGsBenchInput) + (preset : BenchPreset) (gen : StdGen) : IO (BenchGroup × StdGen) := do + let warmup := gsWarmupIterations preset + let measured := gsMeasuredIterations preset + let checksumIterations := measured + let canonicalRecord ← runTimed + (info.groupKey ++ "-canonical") "GuruswamiSudan" "gsListDecode" "KoalaBear.Field" + input.inputShape preset warmup measured + (fun _ => gsListDecode input.code input.received) + checksumDecodeResult (checksumIterations := checksumIterations) + let fastRecord ← runTimed + (info.groupKey ++ "-fast") "GuruswamiSudanFast" "gsListDecode" "KoalaBear.Fast.Field" + fastInput.inputShape preset warmup measured + (fun _ => CompPoly.CodingTheory.GuruswamiSudanFast.gsListDecode + fastInput.code fastInput.received) + checksumFastDecodeResult (checksumIterations := checksumIterations) + pure ({ + groupKey := info.groupKey + title := info.title + records := #[canonicalRecord, fastRecord] + }, gen) + +def gsConsecutiveNoErrorInfo : BenchGroupInfo := + ⟨"gs-koalabear-consecutive-no-error", + "Guruswami-Sudan decode, consecutive domain, no errors (KoalaBear)"⟩ + +def gsConsecutiveOneErrorInfo : BenchGroupInfo := + ⟨"gs-koalabear-consecutive-one-error", + "Guruswami-Sudan decode, consecutive domain, one error (KoalaBear)"⟩ + +def gsConsecutiveRs16NoErrorInfo : BenchGroupInfo := + ⟨"gs-koalabear-consecutive-rs16-no-error", + "Guruswami-Sudan decode, consecutive domain, RS(16,4), no errors (KoalaBear)"⟩ + +def gsRootsNoErrorInfo : BenchGroupInfo := + ⟨"gs-koalabear-roots-no-error", + "Guruswami-Sudan decode, roots-of-unity domain, no errors (KoalaBear)"⟩ + +def gsRootsOneErrorInfo : BenchGroupInfo := + ⟨"gs-koalabear-roots-one-error", + "Guruswami-Sudan decode, roots-of-unity domain, one error (KoalaBear)"⟩ + +def gsGroupInfos : List BenchGroupInfo := [ + gsConsecutiveNoErrorInfo, + gsConsecutiveOneErrorInfo, + gsConsecutiveRs16NoErrorInfo, + gsRootsNoErrorInfo, + gsRootsOneErrorInfo +] + +def codingTheoryTasks : List BenchTask := [ + BenchTask.fromGroupRunner gsConsecutiveNoErrorInfo + (runGsDecodeGroup gsConsecutiveNoErrorInfo (consecutiveInput false) + (fastConsecutiveInput false)), + BenchTask.fromGroupRunner gsConsecutiveOneErrorInfo + (runGsDecodeGroup gsConsecutiveOneErrorInfo (consecutiveInput true) + (fastConsecutiveInput true)), + BenchTask.fromGroupRunner gsConsecutiveRs16NoErrorInfo + (runGsDecodeGroup gsConsecutiveRs16NoErrorInfo consecutiveRs16Input + fastConsecutiveRs16Input), + BenchTask.fromGroupRunner gsRootsNoErrorInfo + (runGsDecodeGroup gsRootsNoErrorInfo (rootsOfUnityInput false) + (fastRootsOfUnityInput false)), + BenchTask.fromGroupRunner gsRootsOneErrorInfo + (runGsDecodeGroup gsRootsOneErrorInfo (rootsOfUnityInput true) + (fastRootsOfUnityInput true)) +] + +end CompPolyBench diff --git a/bench/CompPolyBench/Setup.lean b/bench/CompPolyBench/Setup.lean index 180aeab1..353aa8c3 100644 --- a/bench/CompPolyBench/Setup.lean +++ b/bench/CompPolyBench/Setup.lean @@ -5,6 +5,7 @@ Authors: Valerii Huhnin -/ import CompPolyBench.Bivariate.Basic +import CompPolyBench.CodingTheory import CompPolyBench.Fields.Binary.AdditiveNTT.Impl import CompPolyBench.Multilinear.Basic import CompPolyBench.Multivariate.CMvPolynomial @@ -20,7 +21,8 @@ namespace CompPolyBench /-- Runnable benchmark registry. -/ def allTasks : List BenchTask := - univariateTasks ++ multivariateTasks ++ multilinearTasks ++ bivariateTasks ++ additiveNttTasks + univariateTasks ++ multivariateTasks ++ multilinearTasks ++ bivariateTasks ++ additiveNttTasks ++ + codingTheoryTasks /-- Metadata for every benchmark group accepted by the command-line selector. -/ def allGroupInfos : List BenchGroupInfo := diff --git a/tests/CompPolyTests.lean b/tests/CompPolyTests.lean index 94f2d786..8ed7bd8a 100644 --- a/tests/CompPolyTests.lean +++ b/tests/CompPolyTests.lean @@ -8,6 +8,8 @@ import CompPolyTests.Bivariate.Basic import CompPolyTests.Bivariate.Degree import CompPolyTests.Bivariate.Multiplicity import CompPolyTests.Bivariate.WeightedDegree +import CompPolyTests.CodingTheory.GuruswamiSudan +import CompPolyTests.CodingTheory.GuruswamiSudanFast import CompPolyTests.Data.MvPolynomial.Notation import CompPolyTests.Fields.Binary.AdditiveNTT.NovelPolynomialBasis import CompPolyTests.Fields.Binary.BF128Ghash.Prelude diff --git a/tests/CompPolyTests/CodingTheory/GuruswamiSudan.lean b/tests/CompPolyTests/CodingTheory/GuruswamiSudan.lean new file mode 100644 index 00000000..9b42c671 --- /dev/null +++ b/tests/CompPolyTests/CodingTheory/GuruswamiSudan.lean @@ -0,0 +1,109 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: OpenAI Codex +-/ + +import CompPoly.CodingTheory.GuruswamiSudan + +/-! +# Guruswami-Sudan Decoder Tests + +Executable regression checks for the KoalaBear Guruswami-Sudan port. +-/ + +namespace CompPoly +namespace CodingTheory +namespace GuruswamiSudan + +def hasCandidate (result : DecodeResult) (message : Array F) : Bool := + result.candidates.any fun candidate => UniPoly.coeffsEq candidate message + +#guard NatUtil.floorSqrt 960 = 30 +#guard NatUtil.ceilSqrt 960 = 31 +#guard NatUtil.ceilSqrtDiv (16 * 4 * 5) 4 = 9 +#guard NatUtil.binomial 8 3 = 56 + +#guard countMonomials 10 3 = 22 +#guard gsDecodingRadius 8 2 = 3 +#guard chooseParameters 8 2 = (2, 9) +#guard gsDecodingRadius 16 4 = 7 +#guard chooseParameters 16 4 = (4, 35) + +#guard + valuesToNats (ReedSolomonCode.consecutiveDomain 8) = + #[0, 1, 2, 3, 4, 5, 6, 7] + +#guard + valuesToNats (ReedSolomonCode.pow2Domain 3) = + #[1, 1748172362, 2113994754, 391001680, 2130706432, 382534071, + 16711679, 1739704753] + +#guard + let code := ReedSolomonCode.withConsecutiveDomain 8 2 + let message : Array F := #[3, 5] + valuesToNats (ReedSolomonCode.encode code message) = + #[3, 8, 13, 18, 23, 28, 33, 38] + +#guard + let code := ReedSolomonCode.withConsecutiveDomain 8 2 + let message : Array F := #[1, 2] + let received := ReedSolomonCode.encode code message + let (multiplicity, degreeBound) := chooseParameters code.n code.k + let interpolation := + interpolateWithMultiplicity code.domain received multiplicity degreeBound code.k + UniPoly.isZero (Bivariate.evaluateYPolynomial interpolation message) + +#guard + let code := ReedSolomonCode.withConsecutiveDomain 8 2 + let message : Array F := #[1, 2] + let received := ReedSolomonCode.encode code message + let result := gsListDecode code received + result.multiplicity = 2 && result.degreeBound = 9 && result.errorBound = 3 && + candidatesToNats result.candidates = #[#[1, 2]] && hasCandidate result message + +#guard + let code := ReedSolomonCode.withConsecutiveDomain 8 2 + let message : Array F := #[1, 2] + let codeword := ReedSolomonCode.encode code message + let received := introduceErrorsAtPositions codeword #[1] + let result := gsListDecode code received + valuesToNats received = #[1, 4, 5, 7, 9, 11, 13, 15] && + candidatesToNats result.candidates = #[#[1, 2]] && hasCandidate result message + +#guard + let code := ReedSolomonCode.withRootsOfUnityDomain 3 2 + let message : Array F := #[1, 2] + valuesToNats (ReedSolomonCode.encode code message) = + #[3, 1365638292, 2097283076, 782003361, 2130706432, 765068143, 33423359, + 1348703074] + +#guard + let code := ReedSolomonCode.withRootsOfUnityDomain 3 2 + let message : Array F := #[1, 2] + let received := ReedSolomonCode.encode code message + let result := gsListDecode code received + result.multiplicity = 2 && result.degreeBound = 9 && result.errorBound = 3 && + candidatesToNats result.candidates = #[#[1, 2]] && hasCandidate result message + +#guard + let code := ReedSolomonCode.withRootsOfUnityDomain 3 2 + let message : Array F := #[1, 2] + let codeword := ReedSolomonCode.encode code message + let received := introduceErrorsAtPositions codeword #[2] + let result := gsListDecode code received + valuesToNats received = + #[3, 1365638292, 2097283077, 782003361, 2130706432, 765068143, 33423359, + 1348703074] && + candidatesToNats result.candidates = #[#[1, 2]] && hasCandidate result message + +#guard + let code := ReedSolomonCode.withConsecutiveDomain 8 2 + let message : Array F := #[1, 2] + let received := introduceErrorsAtPositions (ReedSolomonCode.encode code message) #[1] + let result := gsListDecodeWithMultiplicity code received 2 + result.multiplicity = 2 && result.degreeBound = 9 && hasCandidate result message + +end GuruswamiSudan +end CodingTheory +end CompPoly diff --git a/tests/CompPolyTests/CodingTheory/GuruswamiSudanFast.lean b/tests/CompPolyTests/CodingTheory/GuruswamiSudanFast.lean new file mode 100644 index 00000000..cc24bb01 --- /dev/null +++ b/tests/CompPolyTests/CodingTheory/GuruswamiSudanFast.lean @@ -0,0 +1,78 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: OpenAI Codex +-/ + +import CompPoly.CodingTheory.GuruswamiSudanFast + +/-! +# Fast KoalaBear Guruswami-Sudan Decoder Tests + +Executable regression checks for the native-word KoalaBear Guruswami-Sudan port. +-/ + +namespace CompPoly +namespace CodingTheory +namespace GuruswamiSudanFast + +def hasCandidate (result : DecodeResult) (message : Array F) : Bool := + result.candidates.any fun candidate => UniPoly.coeffsEq candidate message + +#guard countMonomials 10 3 = 22 +#guard gsDecodingRadius 8 2 = 3 +#guard chooseParameters 8 2 = (2, 9) +#guard gsDecodingRadius 16 4 = 7 +#guard chooseParameters 16 4 = (4, 35) + +#guard + valuesToNats (ReedSolomonCode.consecutiveDomain 8) = + #[0, 1, 2, 3, 4, 5, 6, 7] + +#guard + valuesToNats (ReedSolomonCode.pow2Domain 3) = + #[1, 1748172362, 2113994754, 391001680, 2130706432, 382534071, + 16711679, 1739704753] + +#guard + let code := ReedSolomonCode.withConsecutiveDomain 8 2 + let message : Array F := #[1, 2] + let received := ReedSolomonCode.encode code message + let result := gsListDecode code received + result.multiplicity = 2 && result.degreeBound = 9 && result.errorBound = 3 && + valuesToNats received = #[1, 3, 5, 7, 9, 11, 13, 15] && + candidatesToNats result.candidates = #[#[1, 2]] && hasCandidate result message + +#guard + let code := ReedSolomonCode.withConsecutiveDomain 8 2 + let message : Array F := #[1, 2] + let codeword := ReedSolomonCode.encode code message + let received := introduceErrorsAtPositions codeword #[1] + let result := gsListDecode code received + valuesToNats received = #[1, 4, 5, 7, 9, 11, 13, 15] && + candidatesToNats result.candidates = #[#[1, 2]] && hasCandidate result message + +#guard + let code := ReedSolomonCode.withRootsOfUnityDomain 3 2 + let message : Array F := #[1, 2] + let received := ReedSolomonCode.encode code message + let result := gsListDecode code received + valuesToNats received = + #[3, 1365638292, 2097283076, 782003361, 2130706432, 765068143, 33423359, + 1348703074] && + candidatesToNats result.candidates = #[#[1, 2]] && hasCandidate result message + +#guard + let code := ReedSolomonCode.withRootsOfUnityDomain 3 2 + let message : Array F := #[1, 2] + let codeword := ReedSolomonCode.encode code message + let received := introduceErrorsAtPositions codeword #[2] + let result := gsListDecode code received + valuesToNats received = + #[3, 1365638292, 2097283077, 782003361, 2130706432, 765068143, 33423359, + 1348703074] && + candidatesToNats result.candidates = #[#[1, 2]] && hasCandidate result message + +end GuruswamiSudanFast +end CodingTheory +end CompPoly