diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b19f7ed..38d7502 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -41,3 +41,13 @@ jobs: python3 scripts/gen_jubjub.py > CompElliptic/Fields/Jubjub.lean - name: Fail if the regenerated files differ from the committed ones run: git diff --exit-code CompElliptic/Fields/Pasta.lean CompElliptic/Fields/Jubjub.lean + + native-lane: + name: FastFieldNative lane is core-only + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@34e114876b0b11c390a56381ad16ebd13914f8d5 # v4.3.1 + with: + persist-credentials: false + - name: Check the precompiled lane's import closure + run: scripts/check_native_lane.sh diff --git a/CompElliptic/Curves/Pasta/Fast/Msm.lean b/CompElliptic/Curves/Pasta/Fast/Msm.lean new file mode 100644 index 0000000..a33c3eb --- /dev/null +++ b/CompElliptic/Curves/Pasta/Fast/Msm.lean @@ -0,0 +1,533 @@ +/- +Copyright (c) 2026 CompElliptic Contributors. All rights reserved. +Released under the Apache License, Version 2.0, or the MIT license, at your option, +as described in the files LICENSE-APACHE and LICENSE-MIT. +Authors: Gregor Mitscha-Baude +-/ +import Mathlib.Data.List.GetD +import CompElliptic.Fields.Pasta + +/-! +# Windowed Pippenger multi-scalar multiplication, proven equal to the naive MSM + +The naive MSM `(terms.map fun (n, x) => n • x).sum` is the mathematical core of a Pedersen-style +vector commitment, and with 255-bit scalars it dominates a verifying-key derivation. This module +gives the windowed Pippenger accelerator `pippenger c terms` — per-window bucket sums, the +suffix-sum trick, a Horner pass across windows — and proves it equal to the naive MSM +(`pippenger_eq_msm`). Measured at `n = 2048`: `242 s` naive against `22.5 s` at the optimal +window `c = 8` (`defaultWindow`); `pippengerFast` replaces the proof-layer per-bucket `filterMap` +by a single-pass `bucketScatter`. + +Everything is generic over `[AddCommMonoid M]`: scalars are `ℕ`, so nothing below needs negation. + +The equality goes through the base-`2 ^ c` digit decomposition `n • x = ∑ i, digit i n • +(base ^ i • x)` (`smul_eq_sum_digits`), an exchange of the sum over terms with the sum over +windows (`list_sum_finset_sum`), the identification of each window's contribution with its +weighted bucket sum (`naiveWindow_eq_buckets`, `foldr_accStep`), and `hornerList_eq`. + +The `Fast` interfaces are provisional: they are not guaranteed to remain public, and +may be folded into the existing API or otherwise changed incompatibly. +-/ + +namespace CompElliptic.Curves.Pasta.Fast.Msm + +open scoped BigOperators + +variable {M : Type*} [AddCommMonoid M] + +/-! ## The naive MSM (the object we accelerate) -/ + +/-- The naive multi-scalar multiplication: `∑ (n, x), n • x`. This is the shape the VK-commitment +`commitLagrangeWith` reduces to. -/ +def naiveMsm (terms : List (ℕ × M)) : M := + (terms.map fun t => t.1 • t.2).sum + +/-! ## Base-`2 ^ c` digit decomposition of a scalar -/ + +/-- The `i`-th base-`base` digit of `n`: `⌊n / base ^ i⌋ mod base`. -/ +def digit (base i n : ℕ) : ℕ := n / base ^ i % base + +/-- Reconstruct `n` from its first `W` base-`base` digits, provided `n < base ^ W`. A bespoke +windowed decomposition (cleaner here than `Nat.digits`). -/ +theorem sum_digit (base : ℕ) (hb : 0 < base) : + ∀ W n, n < base ^ W → n = ∑ i ∈ Finset.range W, digit base i n * base ^ i := by + intro W + induction W with + | zero => intro n hn; simpa using Nat.lt_one_iff.mp hn + | succ W ih => + intro n hn + rw [Finset.sum_range_succ'] + have hdiv : n / base < base ^ W := by + rw [Nat.div_lt_iff_lt_mul hb, ← pow_succ]; exact hn + have hstep : ∀ i ∈ Finset.range W, + digit base (i + 1) n * base ^ (i + 1) = base * (digit base i (n / base) * base ^ i) := by + intro i _ + have hd : digit base (i + 1) n = digit base i (n / base) := by + simp only [digit, pow_succ, Nat.div_div_eq_div_mul, Nat.mul_comm] + rw [hd, pow_succ] + ring + rw [Finset.sum_congr rfl hstep, ← Finset.mul_sum, ← ih (n / base) hdiv] + simp only [digit, pow_zero, Nat.div_one, Nat.mul_one] + exact (Nat.div_add_mod n base).symm + +/-- `n • x = ∑ i, digit i n • (base ^ i • x)` — the per-scalar decomposition that drives the +window exchange. -/ +theorem smul_eq_sum_digits (base : ℕ) (hb : 0 < base) (W n : ℕ) (x : M) (h : n < base ^ W) : + n • x = ∑ i ∈ Finset.range W, digit base i n • (base ^ i • x) := by + conv_lhs => rw [sum_digit base hb W n h] + rw [Finset.sum_smul] + simp_rw [mul_smul] + +/-! ## Buckets and the weighted bucket sum -/ + +/-- The bucket-`b` sum of a digit-tagged point list: `∑ {x | (b, x) ∈ dp}`. A single `filterMap` +pass, so the group additions across all buckets total the list length. -/ +def bucketOf (dp : List (ℕ × M)) (b : ℕ) : M := + (dp.filterMap fun p => if p.1 = b then some p.2 else none).sum + +/-- The direct window value of a digit-tagged list: `∑ (d, x), d • x`. -/ +def naiveWindow (dp : List (ℕ × M)) : M := + (dp.map fun p => p.1 • p.2).sum + +@[simp] theorem bucketOf_nil (b : ℕ) : bucketOf ([] : List (ℕ × M)) b = 0 := rfl + +/-- Prepending `(d, x)` adds `x` to bucket `d` and leaves the others fixed. -/ +theorem bucketOf_cons (d : ℕ) (x : M) (dp : List (ℕ × M)) (b : ℕ) : + bucketOf ((d, x) :: dp) b = (if d = b then x else 0) + bucketOf dp b := by + unfold bucketOf + rw [List.filterMap_cons] + by_cases h : d = b <;> simp [h] + +/-- A single spike: the weighted-bucket sum of a one-hot indicator at digit `d` is `d • x`, provided +`d < base`. -/ +theorem sum_ite_digit (base d : ℕ) (x : M) (hd : d < base) : + ∑ k ∈ Finset.range (base - 1), (k + 1) • (if d = k + 1 then x else (0 : M)) = d • x := by + cases d with + | zero => simp + | succ d' => + have hcongr : ∀ k ∈ Finset.range (base - 1), + (k + 1) • (if d' + 1 = k + 1 then x else (0 : M)) + = if d' = k then (d' + 1) • x else 0 := by + intro k _ + by_cases h : d' = k + · subst h; simp + · simp [h] + rw [Finset.sum_congr rfl hcongr, Finset.sum_ite_eq (Finset.range (base - 1)) d' (fun _ => (d' + 1) • x)] + have : d' ∈ Finset.range (base - 1) := by + rw [Finset.mem_range]; omega + simp [this] + +/-- **The bucket identity.** The direct window value equals the weighted sum of buckets +`∑_{b=1}^{base-1} b • bucketSum b`, provided every digit is `< base`. -/ +theorem naiveWindow_eq_buckets (base : ℕ) (dp : List (ℕ × M)) + (hlt : ∀ p ∈ dp, p.1 < base) : + naiveWindow dp = ∑ k ∈ Finset.range (base - 1), (k + 1) • bucketOf dp (k + 1) := by + induction dp with + | nil => simp [naiveWindow] + | cons p dp ih => + obtain ⟨d, x⟩ := p + have hd : d < base := hlt (d, x) (by simp) + have ih' := ih (fun q hq => hlt q (by simp [hq])) + simp only [naiveWindow, List.map_cons, List.sum_cons] at ih' ⊢ + rw [ih'] + have hb : ∀ k ∈ Finset.range (base - 1), + (k + 1) • bucketOf ((d, x) :: dp) (k + 1) + = (k + 1) • (if d = k + 1 then x else 0) + (k + 1) • bucketOf dp (k + 1) := by + intro k _ + rw [bucketOf_cons, smul_add] + rw [Finset.sum_congr rfl hb, Finset.sum_add_distrib, sum_ite_digit base d x hd] + +/-! ## The suffix-sum accumulation -/ + +/-- One step of the suffix-sum accumulator: carry `(running, total)`, add the next bucket to the +running sum, then fold the running sum into the total. Folding this from the top bucket down +computes `∑ b, b • bucket_b` in `2 · #buckets` adds. -/ +def accStep (a : M) (p : M × M) : M × M := (p.1 + a, p.2 + (p.1 + a)) + +/-- The mathematical spec of the accumulator's total: the position-weighted sum `∑_k (k+1) • vals[k]` +(`vals[0]` weighted `1`, …). Recursive so the `foldr` invariant is a one-line induction — this is a +*spec*, never executed (it recomputes `vals.sum`); the executable path is `accStep`. -/ +def weightedSum : List M → M + | [] => 0 + | a :: vals => a + vals.sum + weightedSum vals + +/-- The `foldr accStep` invariant: running is the plain sum, total is the position-weighted sum +(shifted by the initial carry). -/ +theorem foldr_accStep (r0 t0 : M) (vals : List M) : + List.foldr accStep (r0, t0) vals + = (r0 + vals.sum, t0 + vals.length • r0 + weightedSum vals) := by + induction vals with + | nil => simp [weightedSum] + | cons a vals ih => + rw [List.foldr_cons, ih, accStep] + simp only [List.sum_cons, List.length_cons, weightedSum, succ_nsmul] + refine Prod.ext ?_ ?_ <;> simp only <;> abel + +/-- Weighted sum of an appended singleton: the new element sits at the end with weight `length + 1`. -/ +theorem weightedSum_append (xs : List M) (y : M) : + weightedSum (xs ++ [y]) = weightedSum xs + (xs.length + 1) • y := by + induction xs with + | nil => simp [weightedSum] + | cons a xs ih => + simp only [List.cons_append, weightedSum, List.sum_append, List.sum_cons, List.sum_nil, + add_zero, ih, List.length_cons, succ_nsmul] + abel + +/-- Weighted sum of a `range`-indexed list is the finite weighted sum. -/ +theorem weightedSum_range_map (K : ℕ) (g : ℕ → M) : + weightedSum ((List.range K).map g) = ∑ k ∈ Finset.range K, (k + 1) • g k := by + induction K with + | zero => simp [weightedSum] + | succ K ih => + rw [List.range_succ, List.map_append, List.map_cons, List.map_nil, weightedSum_append, + List.length_map, List.length_range, ih, Finset.sum_range_succ] + +/-! ## The Horner fold across windows -/ + +/-- Horner across windows, least-significant first: `∑ i, base ^ i • vals[i]` in `#windows` +doublings. `foldr` processes the most-significant window first (innermost), each step doubles +(`base •`) and adds the window value. -/ +def hornerList (base : ℕ) (vals : List M) : M := + List.foldr (fun v acc => base • acc + v) 0 vals + +theorem hornerList_eq (base : ℕ) (vals : List M) : + hornerList base vals = ∑ k ∈ Finset.range vals.length, base ^ k • vals.getD k 0 := by + induction vals with + | nil => simp [hornerList] + | cons v vals ih => + rw [hornerList, List.foldr_cons] + show base • hornerList base vals + v = _ + rw [ih, List.length_cons, Finset.sum_range_succ'] + simp only [List.getD_cons_zero, List.getD_cons_succ, pow_zero, one_smul, pow_succ'] + rw [Finset.smul_sum] + simp only [mul_smul] + +/-! ## The full Pippenger MSM -/ + +/-- The digit-tagged term list for window `i`: each scalar replaced by its `i`-th base-`base` digit. -/ +def dpOf (base i : ℕ) (terms : List (ℕ × M)) : List (ℕ × M) := + terms.map fun t => (digit base i t.1, t.2) + +/-- The window `i` value, computed the fast way: build the `base − 1` buckets, then run the +suffix-sum accumulation (`2 · (base − 1)` adds). -/ +def windowValue (base i : ℕ) (terms : List (ℕ × M)) : M := + (List.foldr accStep (0, 0) + ((List.range (base - 1)).map fun k => bucketOf (dpOf base i terms) (k + 1))).2 + +/-- The largest scalar appearing in `terms` (`0` if empty). -/ +def maxScalar (terms : List (ℕ × M)) : ℕ := (terms.map Prod.fst).foldr max 0 + +omit [AddCommMonoid M] in +theorem le_maxScalar (terms : List (ℕ × M)) (t : ℕ × M) (ht : t ∈ terms) : + t.1 ≤ maxScalar terms := by + unfold maxScalar + have : t.1 ∈ terms.map Prod.fst := List.mem_map_of_mem ht + generalize terms.map Prod.fst = l at this + clear ht + induction l with + | nil => simp at this + | cons a l ih => + rcases List.mem_cons.mp this with h | h + · subst h; exact le_max_left _ _ + · exact le_trans (ih h) (le_max_right _ _) + +/-- The number of windows: enough base-`2 ^ c` digits to cover the largest scalar. -/ +def numWindows (c : ℕ) (terms : List (ℕ × M)) : ℕ := + Nat.log (2 ^ c) (maxScalar terms) + 1 + +/-- **Windowed Pippenger MSM.** Provably equal to `naiveMsm` (`pippenger_eq_msm`). -/ +def pippenger (c : ℕ) (terms : List (ℕ × M)) : M := + hornerList (2 ^ c) + ((List.range (numWindows c terms)).map fun i => windowValue (2 ^ c) i terms) + +/-! ## Auxiliary sum-exchange lemmas -/ + +/-- Push a scalar through a mapped list sum. -/ +theorem smul_list_sum (c : ℕ) (l : List M) : + c • l.sum = (l.map fun a => c • a).sum := by + induction l with + | nil => simp + | cons a l ih => simp [smul_add, ih] + +/-- Fubini between a list sum and a finite sum. -/ +theorem list_sum_finset_sum {α ι : Type*} (l : List α) (s : Finset ι) (F : α → ι → M) : + (l.map fun a => ∑ i ∈ s, F a i).sum = ∑ i ∈ s, (l.map fun a => F a i).sum := by + induction l with + | nil => simp + | cons a l ih => simp only [List.map_cons, List.sum_cons, ih, ← Finset.sum_add_distrib] + +/-- `getD` of a `range`-map at an in-range index. -/ +theorem getD_map_range (f : ℕ → M) (W k : ℕ) (h : k < W) : + ((List.range W).map f).getD k 0 = f k := by + rw [List.getD_eq_getElem?_getD] + simp [h] + +/-! ## The main equality -/ + +/-- The fast window value equals the direct window contribution `∑ (n, x), digit(n) • x`. -/ +theorem windowValue_eq (base i : ℕ) (hbase : 0 < base) (terms : List (ℕ × M)) : + windowValue base i terms = (terms.map fun t => digit base i t.1 • t.2).sum := by + unfold windowValue + rw [foldr_accStep] + simp only [smul_zero, add_zero, zero_add] + rw [weightedSum_range_map] + have hdp : ∀ p ∈ dpOf base i terms, p.1 < base := by + intro p hp + unfold dpOf at hp + rw [List.mem_map] at hp + obtain ⟨t, _, rfl⟩ := hp + exact Nat.mod_lt _ hbase + rw [← naiveWindow_eq_buckets base (dpOf base i terms) hdp] + unfold naiveWindow dpOf + rw [List.map_map] + rfl + +/-- **Pippenger equals the naive MSM.** For any window size `c ≥ 1` and any term list. -/ +theorem pippenger_eq_msm (c : ℕ) (hc : 0 < c) (terms : List (ℕ × M)) : + pippenger c terms = (terms.map fun t => t.1 • t.2).sum := by + have hb1 : 1 < 2 ^ c := Nat.one_lt_two_pow_iff.mpr hc.ne' + have hb0 : 0 < 2 ^ c := by positivity + set base := 2 ^ c + set W := numWindows c terms with hW + -- every scalar fits in `W` base-`base` digits + have hbound : ∀ t ∈ terms, t.1 < base ^ W := by + intro t ht + calc t.1 ≤ maxScalar terms := le_maxScalar terms t ht + _ < base ^ W := Nat.lt_pow_succ_log_self hb1 _ + -- the Pippenger value is the Horner sum of the true window contributions + have hpip : pippenger c terms + = ∑ i ∈ Finset.range W, base ^ i • (terms.map fun t => digit base i t.1 • t.2).sum := by + unfold pippenger + rw [hornerList_eq, List.length_map, List.length_range] + apply Finset.sum_congr rfl + intro k hk + rw [Finset.mem_range] at hk + rw [getD_map_range _ _ _ hk, windowValue_eq base k hb0] + rw [hpip] + -- the naive MSM is the same Horner sum, via digit decomposition + sum exchange + have e1 : (terms.map fun t => t.1 • t.2) + = terms.map fun t => ∑ i ∈ Finset.range W, digit base i t.1 • (base ^ i • t.2) := by + apply List.map_congr_left + intro t ht + exact smul_eq_sum_digits base hb0 W t.1 t.2 (hbound t ht) + rw [e1, list_sum_finset_sum] + apply Finset.sum_congr rfl + intro i _ + rw [smul_list_sum, List.map_map] + apply congrArg List.sum + apply List.map_congr_left + intro t _ + simp only [Function.comp] + rw [smul_comm] + +/-! ## Single-pass Array bucketing (constant factors) + +`windowValue` computes each bucket by its own `filterMap` pass — `(base − 1) · n` list steps per +window, which measurably competes with the group ops at `base = 2 ^ 8` and dominates beyond. +`bucketScatter` builds all buckets in ONE pass (an `Array.modify` per term) and is proven to +produce exactly the `bucketOf` list (`bucketScatter_toList`), so +`pippengerFast = pippenger = naiveMsm`. -/ + +/-- One scatter step: drop digit-`0` terms, otherwise add the point into bucket slot `d − 1` +(slot `k` holds bucket `k + 1`). -/ +def scatterStep (a : Array M) (p : ℕ × M) : Array M := + if p.1 = 0 then a else a.modify (p.1 - 1) (· + p.2) + +/-- Scatter a digit-tagged point list into its `base − 1` buckets in one pass. -/ +def bucketScatter (base : ℕ) (dp : List (ℕ × M)) : Array M := + dp.foldl scatterStep (Array.replicate (base - 1) 0) + +@[simp] theorem size_scatterStep (a : Array M) (p : ℕ × M) : + (scatterStep a p).size = a.size := by + unfold scatterStep; split <;> simp + +@[simp] theorem size_foldl_scatterStep (dp : List (ℕ × M)) (a : Array M) : + (dp.foldl scatterStep a).size = a.size := by + induction dp generalizing a with + | nil => rfl + | cons p dp ih => rw [List.foldl_cons, ih, size_scatterStep] + +/-- A single scatter step read back through `getElem?` (proof-free indexing, so the fold invariant +rewrites cleanly): slot `k` gains `x` exactly when the digit is `k + 1`. -/ +theorem getElem?_scatterStep (a : Array M) (d : ℕ) (x : M) (k : ℕ) : + (scatterStep a (d, x))[k]? + = a[k]?.map fun v => v + (if d = k + 1 then x else 0) := by + unfold scatterStep + by_cases h0 : d = 0 + · simp [h0, Option.map_id'] + · simp only [if_neg h0] + rw [Array.getElem?_modify] + by_cases hdk : d = k + 1 + · rw [if_pos (by omega), if_pos hdk] + · rw [if_neg (by omega), if_neg hdk] + simp [Option.map_id'] + +/-- The fold invariant: slot `k` accumulates bucket `k + 1` on top of its initial value. -/ +theorem getElem?_foldl_scatterStep (dp : List (ℕ × M)) : + ∀ (a : Array M) (k : ℕ), + (dp.foldl scatterStep a)[k]? = a[k]?.map fun v => v + bucketOf dp (k + 1) := by + induction dp with + | nil => intro a k; simp [Option.map_id'] + | cons p dp ih => + intro a k + obtain ⟨d, x⟩ := p + rw [List.foldl_cons, ih (scatterStep a (d, x)) k, getElem?_scatterStep a d x k, + Option.map_map, bucketOf_cons] + apply congrFun + apply congrArg + funext v + simp only [Function.comp_apply] + rw [add_assoc] + +/-- `bucketScatter` produces exactly the per-bucket list `windowValue` folds over. -/ +theorem bucketScatter_toList (base : ℕ) (dp : List (ℕ × M)) : + (bucketScatter base dp).toList + = (List.range (base - 1)).map fun k => bucketOf dp (k + 1) := by + apply List.ext_getElem? + intro k + rw [Array.getElem?_toList, bucketScatter, getElem?_foldl_scatterStep dp _ k, + List.getElem?_map, Array.getElem?_replicate] + by_cases hk : k < base - 1 + · rw [if_pos hk, List.getElem?_range hk] + simp + · rw [if_neg hk, List.getElem?_eq_none (by simp; omega)] + simp + +/-- The window value via the single-pass scatter. -/ +def windowValueFast (base i : ℕ) (terms : List (ℕ × M)) : M := + (List.foldr accStep (0, 0) (bucketScatter base (dpOf base i terms)).toList).2 + +theorem windowValueFast_eq (base i : ℕ) (terms : List (ℕ × M)) : + windowValueFast base i terms = windowValue base i terms := by + unfold windowValueFast windowValue + rw [bucketScatter_toList] + +/-- **Windowed Pippenger MSM, Array-bucketed** — the executable accelerator. Equal to `pippenger` +(`pippengerFast_eq`) and hence to the naive MSM (`pippengerFast_eq_msm`). -/ +def pippengerFast (c : ℕ) (terms : List (ℕ × M)) : M := + hornerList (2 ^ c) + ((List.range (numWindows c terms)).map fun i => windowValueFast (2 ^ c) i terms) + +theorem pippengerFast_eq (c : ℕ) (terms : List (ℕ × M)) : + pippengerFast c terms = pippenger c terms := by + unfold pippengerFast pippenger + congr 1 + apply List.map_congr_left + intro i _ + exact windowValueFast_eq (2 ^ c) i terms + +/-- The Array-bucketed Pippenger equals the naive MSM — the accelerator's end-to-end +correctness. -/ +theorem pippengerFast_eq_msm (c : ℕ) (hc : 0 < c) (terms : List (ℕ × M)) : + pippengerFast c terms = (terms.map fun t => t.1 • t.2).sum := by + rw [pippengerFast_eq, pippenger_eq_msm c hc] + +/-! ## Windows-parallel Pippenger + +The `numWindows c terms ≈ ⌈255/c⌉` window values are mutually independent (each scans the term +list, buckets, and runs its own suffix accumulation); only the final Horner fold consumes them in +order. `pippengerFastPar` evaluates the windows through the proven `parMap` below, so the +equality to the sequential accelerator is purely `parMap_eq_map` — the evaluation strategy is +the only difference. Measured (interpreted, `n = 2048`, `c = 8`, 12 cores): +`11.6 s → 4.3 s` per MSM. -/ + +/-- Parallel `List.map`: spawn a task per element, then collect the results in order. + +This is a *pure* function — `Task` is logically the one-field structure `⟨get : α⟩`, +`Task.spawn fn = ⟨fn ()⟩` and `Task.get ⟨a⟩ = a` are both definitional — so it is provably equal +to `List.map` (`parMap_eq_map`). The only thing it changes is the *evaluation strategy*: under +`native_decide` (compiled to C, run on the real task scheduler) the elementwise work runs on the +thread pool. Because the equality to `List.map` is definitional, swapping `List.map` for `parMap` +in a `native_decide`-backed definition changes nothing about the statement's meaning or trust +base: the compiled `Task` primitives are already inside the `native_decide` trust boundary, and +the kernel's view is literally `List.map`. -/ +def parMap {α : Type u} {β : Type v} (f : α → β) (xs : List α) : List β := + (xs.map (fun x => Task.spawn (fun _ => f x))).map Task.get + +/-- `parMap` is `List.map`: the task round-trip `Task.get (Task.spawn (fun _ => f x))` is `f x` +definitionally, so nothing but the evaluation strategy differs. -/ +@[simp] theorem parMap_eq_map {α : Type u} {β : Type v} (f : α → β) (xs : List α) : + parMap f xs = xs.map f := by + rw [parMap, List.map_map]; rfl + +/-- **Windows-parallel Array-bucketed Pippenger MSM**: `pippengerFast` with the independent +window values evaluated as parallel tasks. Equal to `pippengerFast` (`pippengerFastPar_eq`) +and hence to the naive MSM (`pippengerFastPar_eq_msm`). -/ +def pippengerFastPar (c : ℕ) (terms : List (ℕ × M)) : M := + hornerList (2 ^ c) + (parMap (fun i => windowValueFast (2 ^ c) i terms) (List.range (numWindows c terms))) + +/-- The windows-parallel Pippenger is the sequential one: `parMap` is `map`. -/ +theorem pippengerFastPar_eq (c : ℕ) (terms : List (ℕ × M)) : + pippengerFastPar c terms = pippengerFast c terms := by + rw [pippengerFastPar, parMap_eq_map, pippengerFast] + +/-- The windows-parallel Pippenger equals the naive MSM. -/ +theorem pippengerFastPar_eq_msm (c : ℕ) (hc : 0 < c) (terms : List (ℕ × M)) : + pippengerFastPar c terms = (terms.map fun t => t.1 • t.2).sum := by + rw [pippengerFastPar_eq, pippengerFast_eq_msm c hc] + +/-! ## The `commit_lagrange` wrapper -/ + +/-- The default window size, picked by measurement (see the module docstring's ladder): `c = 8` +is the optimum for the `n = 2048`-term MSMs of an Orchard-sized verifying key — `10.7×` over the +naive MSM. -/ +def defaultWindow : ℕ := 8 + +section Wrapper + +variable {G : Type*} [AddCommGroup G] [Inhabited G] + +/-- The scalar field the commitment coefficients live in: Vesta's scalar field (`= pallas::Base`), +which is what a Vesta-based Pedersen commitment takes its coefficients from. -/ +abbrev Fp := CompElliptic.Fields.Pasta.VestaScalarField + +/-- The naive `commit_lagrange` spec — a Pedersen vector commitment against a fixed Lagrange +basis, plus a blinding term: +`(∑ᵢ coeffsᵢ.val • basisᵢ) + blind`. -/ +def commitLagrangeSpec (blind : G) (basis : List G) (coeffs : List Fp) : G := + ((List.range coeffs.length).map + (fun i => (coeffs.getD i 0).val • basis.getD i 0)).sum + blind + +/-- The fast `commit_lagrange`: Array-bucketed windowed Pippenger over the `(coeffᵢ.val, basisᵢ)` +terms, plus the blind. The terms are zipped directly (`List.zip`), not indexed through `getD`, to +avoid the quadratic list-indexing of the spec's `range`/`getD` spelling. -/ +def commitLagrangeFastWith (c : ℕ) (blind : G) (basis : List G) (coeffs : List Fp) : G := + pippengerFast c + ((coeffs.zip (basis ++ List.replicate (coeffs.length - basis.length) 0)).map + fun t => (t.1.val, t.2)) + blind + +omit [Inhabited G] in +/-- The zipped term list is the spec's `range`/`getD` term list: the padded basis makes the zip +total at length `coeffs.length`, and in-range `getD` is `getElem`. -/ +theorem zip_terms_eq (basis : List G) (coeffs : List Fp) : + ((coeffs.zip (basis ++ List.replicate (coeffs.length - basis.length) 0)).map + fun t => (t.1.val, t.2)) + = (List.range coeffs.length).map + fun i => ((coeffs.getD i 0).val, basis.getD i 0) := by + apply List.ext_getElem + · simp [List.length_zip]; omega + · intro i h1 h2 + have hi : i < coeffs.length := by + simp only [List.length_map, List.length_zip] at h1 + omega + simp only [List.getElem_map, List.getElem_zip, List.getElem_range] + have hc : coeffs.getD i 0 = coeffs[i] := List.getD_eq_getElem coeffs 0 hi + rw [hc] + rcases Nat.lt_or_ge i basis.length with hib | hib + · rw [List.getElem_append_left hib, List.getD_eq_getElem basis 0 hib] + · rw [List.getElem_append_right hib, List.getD_eq_default basis 0 (by omega)] + simp + +omit [Inhabited G] in +/-- **The fast commitment equals the naive spec.** -/ +theorem commitLagrangeFastWith_eq (c : ℕ) (hc : 0 < c) + (blind : G) (basis : List G) (coeffs : List Fp) : + commitLagrangeFastWith c blind basis coeffs = commitLagrangeSpec blind basis coeffs := by + unfold commitLagrangeFastWith commitLagrangeSpec + rw [zip_terms_eq, pippengerFast_eq, pippenger_eq_msm c hc, List.map_map] + rfl + +end Wrapper + +end CompElliptic.Curves.Pasta.Fast.Msm diff --git a/CompElliptic/Curves/Pasta/Fast/MsmProj.lean b/CompElliptic/Curves/Pasta/Fast/MsmProj.lean new file mode 100644 index 0000000..0fd0481 --- /dev/null +++ b/CompElliptic/Curves/Pasta/Fast/MsmProj.lean @@ -0,0 +1,481 @@ +/- +Copyright (c) 2026 CompElliptic Contributors. All rights reserved. +Released under the Apache License, Version 2.0, or the MIT license, at your option, +as described in the files LICENSE-APACHE and LICENSE-MIT. +Authors: Gregor Mitscha-Baude +-/ +import CompElliptic.Curves.Pasta.Fast.Msm +import CompElliptic.Curves.Pasta.Fast.Projective + +/-! +# The windowed Pippenger MSM run entirely in projective coordinates + +Instantiated at the Vesta affine group, every group addition of `Msm.lean` pays a field inversion. +This module runs the whole Pippenger interior over `Projective.PVes` with `padd`/`pnsmulFast` — +one inversion per MSM instead of one per add — and proves the result equal to the affine +accelerator (`pippengerProj_eq`, `pippengerProjScatter_eq`, and `..._eq_msm` to the naive MSM). + +`PVes` is not a lawful monoid on the nose (associativity only holds after `toAffine`), so the +Pippenger equality is not re-run here: each fold of `Msm.lean` is mirrored by a `padd`/`pid` +version and shown to commute with `toAffine`, carrying `Valid` through every intermediate +(`psum_spec`, `pbucketOf_spec`, `foldr_paccStep_spec`, `pwindowValue_spec`, `phornerList_spec`). + +The `Fast` interfaces are provisional: they are not guaranteed to remain public, and +may be folded into the existing API or otherwise changed incompatibly. +-/ + +open CompElliptic.Curves.Pasta.Fast +open CompElliptic.Curves.Pasta.Fast.Msm (Fp) +open CompElliptic.Curves.Pasta.Fast.Projective +open CompElliptic.Curves.Pasta.Fast.Projective.PVes +open CompElliptic +open CompElliptic.CurveForms.ShortWeierstrass +open CompElliptic.Curves.Pasta + +namespace CompElliptic.Curves.Pasta.Fast.MsmProj + +/-- The Vesta affine group carries a default point, needed to instantiate the `Fast.Msm` wrapper +(`commitLagrangeSpec`, `commitLagrangeFastWith_eq`) whose section fixes `[Inhabited G]`. -/ +local instance : Inhabited Projective.G := ⟨0⟩ + +/-! ## Projective bucket sum (`foldr padd pid`), mirroring `Fast.Msm.bucketOf` -/ + +/-- The `padd`-fold sum of a list of projective points, from the projective identity `pid`. -/ +def psum (L : List PVes) : PVes := L.foldr padd pid + +theorem psum_cons (a : PVes) (xs : List PVes) : psum (a :: xs) = padd a (psum xs) := rfl + +/-- A projective `padd`-sum of valid points is valid, and `toAffine` carries it to the affine sum. -/ +theorem psum_spec (L : List PVes) (h : ∀ P ∈ L, Valid P) : + Valid (psum L) ∧ toAffine (psum L) = (L.map toAffine).sum := by + induction L with + | nil => exact ⟨valid_pid, by simp [psum, toAffine_pid]⟩ + | cons a xs ih => + have ha : Valid a := h a (by simp) + have hxs : ∀ P ∈ xs, Valid P := fun P hP => h P (by simp [hP]) + obtain ⟨hvxs, hxseq⟩ := ih hxs + rw [psum_cons] + exact ⟨valid_padd ha hvxs, by rw [toAffine_padd ha hvxs, hxseq, List.map_cons, List.sum_cons]⟩ + +/-- The projective bucket-`b` sum: `padd`-fold the points whose digit tag is `b`. -/ +def pbucketOf (dp : List (ℕ × PVes)) (b : ℕ) : PVes := + psum (dp.filterMap fun p => if p.1 = b then some p.2 else none) + +/-- The projective bucket sum is valid and matches `Fast.Msm.bucketOf` of the `toAffine`-mapped +digit-tagged list. -/ +theorem pbucketOf_spec (dp : List (ℕ × PVes)) (b : ℕ) (h : ∀ p ∈ dp, Valid p.2) : + Valid (pbucketOf dp b) + ∧ toAffine (pbucketOf dp b) + = Msm.bucketOf (dp.map fun t => (t.1, toAffine t.2)) b := by + have hval : ∀ P ∈ dp.filterMap (fun p => if p.1 = b then some p.2 else none), Valid P := by + intro P hP + rw [List.mem_filterMap] at hP + obtain ⟨q, hq, hqeq⟩ := hP + by_cases hc : q.1 = b + · rw [if_pos hc, Option.some.injEq] at hqeq; rw [← hqeq]; exact h q hq + · rw [if_neg hc] at hqeq; exact absurd hqeq (by simp) + obtain ⟨hv, heq⟩ := psum_spec _ hval + refine ⟨hv, ?_⟩ + rw [pbucketOf, heq, Msm.bucketOf, List.filterMap_map, List.map_filterMap] + refine congrArg List.sum (congrArg (dp.filterMap ·) ?_) + funext p + by_cases hc : p.1 = b <;> simp [hc, Function.comp] + +/-! ## Projective suffix-sum accumulator, mirroring `Fast.Msm.accStep` -/ + +/-- One step of the projective suffix-sum accumulator: carry `(running, total)`, `padd` the next +bucket into `running`, then `padd` `running` into `total`. -/ +def paccStep (a : PVes) (p : PVes × PVes) : PVes × PVes := + (padd p.1 a, padd p.2 (padd p.1 a)) + +/-- Folding `paccStep` over valid points keeps both accumulator components valid, and `toAffine` +carries the whole fold (componentwise) to the affine `Fast.Msm.accStep` fold. -/ +theorem foldr_paccStep_spec (L : List PVes) (h : ∀ P ∈ L, Valid P) : + Valid (L.foldr paccStep (pid, pid)).1 ∧ Valid (L.foldr paccStep (pid, pid)).2 + ∧ (toAffine (L.foldr paccStep (pid, pid)).1, toAffine (L.foldr paccStep (pid, pid)).2) + = List.foldr Msm.accStep ((0 : Projective.G), (0 : Projective.G)) (L.map toAffine) := by + induction L with + | nil => exact ⟨valid_pid, valid_pid, by simp [toAffine_pid]⟩ + | cons a xs ih => + have ha : Valid a := h a (by simp) + have hxs : ∀ P ∈ xs, Valid P := fun P hP => h P (by simp [hP]) + obtain ⟨hv1, hv2, heq⟩ := ih hxs + have hstep : (a :: xs).foldr paccStep (pid, pid) = paccStep a (xs.foldr paccStep (pid, pid)) := + rfl + rw [hstep] + have hv1' : Valid (padd (xs.foldr paccStep (pid, pid)).1 a) := valid_padd hv1 ha + refine ⟨hv1', valid_padd hv2 hv1', ?_⟩ + rw [List.map_cons, List.foldr_cons, ← heq, Msm.accStep] + simp only [paccStep, toAffine_padd hv1 ha, toAffine_padd hv2 hv1'] + +/-! ## Projective window value, mirroring `Fast.Msm.windowValue` -/ + +/-- The digit-tagged projective term list for window `i`. -/ +def pdpOf (base i : ℕ) (pterms : List (ℕ × PVes)) : List (ℕ × PVes) := + pterms.map fun t => (Msm.digit base i t.1, t.2) + +/-- The projective window value: build the `base − 1` projective buckets, then run the projective +suffix-sum accumulation. -/ +def pwindowValue (base i : ℕ) (pterms : List (ℕ × PVes)) : PVes := + (List.foldr paccStep (pid, pid) + ((List.range (base - 1)).map fun k => pbucketOf (pdpOf base i pterms) (k + 1))).2 + +/-- The projective window value is valid and matches `Fast.Msm.windowValue` of the `toAffine`-mapped +terms. -/ +theorem pwindowValue_spec (base i : ℕ) (pterms : List (ℕ × PVes)) + (h : ∀ p ∈ pterms, Valid p.2) : + Valid (pwindowValue base i pterms) + ∧ toAffine (pwindowValue base i pterms) + = Msm.windowValue base i (pterms.map fun t => (t.1, toAffine t.2)) := by + have hdp : ∀ p ∈ pdpOf base i pterms, Valid p.2 := by + intro p hp + rw [pdpOf, List.mem_map] at hp + obtain ⟨t, ht, rfl⟩ := hp + exact h t ht + simp only [pwindowValue] + set buckets := (List.range (base - 1)).map fun k => pbucketOf (pdpOf base i pterms) (k + 1) + with hbuck + have hbval : ∀ P ∈ buckets, Valid P := by + intro P hP + rw [hbuck, List.mem_map] at hP + obtain ⟨k, _, rfl⟩ := hP + exact (pbucketOf_spec _ _ hdp).1 + obtain ⟨_, hv2, heq⟩ := foldr_paccStep_spec buckets hbval + refine ⟨hv2, ?_⟩ + have hmap : buckets.map toAffine + = (List.range (base - 1)).map fun k => + Msm.bucketOf (Msm.dpOf base i (pterms.map fun t => (t.1, toAffine t.2))) (k + 1) := by + rw [hbuck, List.map_map] + apply List.map_congr_left + intro k _ + rw [Function.comp_apply, (pbucketOf_spec _ _ hdp).2] + congr 1 + simp only [pdpOf, Msm.dpOf, List.map_map, Function.comp_def] + rw [Msm.windowValue, ← hmap, ← heq] + +/-! ## Projective Horner fold, mirroring `Fast.Msm.hornerList` -/ + +/-- Horner across windows in projective coordinates: each `base •` doubling is the fast binary +scalar multiplication `pnsmulFast base`, and the window value is folded in with `padd`. -/ +def phornerList (base : ℕ) (vals : List PVes) : PVes := + List.foldr (fun v acc => padd (pnsmulFast base acc) v) pid vals + +/-- The projective Horner fold is valid and matches `Fast.Msm.hornerList` of the `toAffine`-mapped +window values; each `2^c`-fold doubling transports via `pnsmulFast_spec`. -/ +theorem phornerList_spec (base : ℕ) (vals : List PVes) (h : ∀ P ∈ vals, Valid P) : + Valid (phornerList base vals) + ∧ toAffine (phornerList base vals) = Msm.hornerList base (vals.map toAffine) := by + induction vals with + | nil => exact ⟨valid_pid, by simp [phornerList, Msm.hornerList, toAffine_pid]⟩ + | cons v xs ih => + have hv : Valid v := h v (by simp) + have hxs : ∀ P ∈ xs, Valid P := fun P hP => h P (by simp [hP]) + obtain ⟨hvacc, heq⟩ := ih hxs + have hstep : phornerList base (v :: xs) + = padd (pnsmulFast base (phornerList base xs)) v := rfl + rw [hstep] + refine ⟨valid_padd (pnsmulFast_spec hvacc base).1 hv, ?_⟩ + rw [toAffine_padd (pnsmulFast_spec hvacc base).1 hv, (pnsmulFast_spec hvacc base).2, heq, + List.map_cons] + simp only [Msm.hornerList, List.foldr_cons] + +/-! ## The projective Pippenger MSM and its equality to the affine accelerator -/ + +/-- **Windowed Pippenger MSM run in projective coordinates.** Each affine point is lifted via +`ofAffine`; the same windowed algorithm (digits, buckets, suffix-sum, Horner doublings) runs over +`PVes` with `padd`/`pnsmulFast`; a single `toAffine` (one field inversion) is taken at the end. +Proven equal to the affine `Fast.Msm.pippenger` (`pippengerProj_eq`). -/ +def pippengerProj (c : ℕ) (terms : List (ℕ × Projective.G)) : Projective.G := + toAffine (phornerList (2 ^ c) + ((List.range (Msm.numWindows c terms)).map fun i => + pwindowValue (2 ^ c) i (terms.map fun t => (t.1, ofAffine t.2)))) + +/-- **The projective Pippenger equals the affine Pippenger.** Lifting by `ofAffine` and reading back +by `toAffine` is the identity on terms (`toAffine_ofAffine`); each window and the Horner fold commute +with `toAffine` (`pwindowValue_spec`, `phornerList_spec`), so the whole projective interior collapses +to the affine `Fast.Msm.pippenger`. -/ +theorem pippengerProj_eq (c : ℕ) (terms : List (ℕ × Projective.G)) : + pippengerProj c terms = Msm.pippenger c terms := by + rw [pippengerProj] + set pterms := terms.map (fun t => (t.1, ofAffine t.2)) with hpterms + set windows := (List.range (Msm.numWindows c terms)).map fun i => pwindowValue (2 ^ c) i pterms + with hwin + have hval : ∀ p ∈ pterms, Valid p.2 := by + intro p hp + rw [hpterms, List.mem_map] at hp + obtain ⟨t, _, rfl⟩ := hp + exact valid_ofAffine t.2 + have hf : pterms.map (fun t => (t.1, toAffine t.2)) = terms := by + rw [hpterms, List.map_map] + conv_rhs => rw [← List.map_id terms] + apply List.map_congr_left + intro t _ + simp only [Function.comp_apply, toAffine_ofAffine, id_eq] + have hwval : ∀ P ∈ windows, Valid P := by + intro P hP + rw [hwin, List.mem_map] at hP + obtain ⟨i, _, rfl⟩ := hP + exact (pwindowValue_spec (2 ^ c) i pterms hval).1 + rw [(phornerList_spec (2 ^ c) windows hwval).2, Msm.pippenger] + congr 1 + rw [hwin, List.map_map] + apply List.map_congr_left + intro i _ + rw [Function.comp_apply, (pwindowValue_spec (2 ^ c) i pterms hval).2, hf] + +/-- **The projective Pippenger equals the naive MSM.** End-to-end correctness for `c ≥ 1`. -/ +theorem pippengerProj_eq_msm (c : ℕ) (hc : 0 < c) (terms : List (ℕ × Projective.G)) : + pippengerProj c terms = (terms.map fun t => t.1 • t.2).sum := by + rw [pippengerProj_eq, Msm.pippenger_eq_msm c hc] + +/-! ## Single-pass Array bucketing for the projective interior + +`pwindowValue`'s per-bucket `filterMap` dominates the group work at `base = 2 ^ 8` (~5× slower +than the scatter below). `pbucketScatter` builds all buckets in one pass; its invariant cannot +reuse `Msm.bucketScatter_toList` (that one uses `add_assoc`, and `padd` is not associative on the +nose), so it too is transported through `toAffine` (`foldl_pscatterStep_spec`). -/ + +/-- One projective scatter step: drop digit-`0` terms, otherwise `padd` the point into bucket +slot `d − 1` (slot `k` holds bucket `k + 1`). -/ +def pscatterStep (a : Array PVes) (p : ℕ × PVes) : Array PVes := + if p.1 = 0 then a else a.modify (p.1 - 1) (fun v => padd v p.2) + +/-- Scatter a digit-tagged projective point list into its `base − 1` buckets in one pass. -/ +def pbucketScatter (base : ℕ) (dp : List (ℕ × PVes)) : Array PVes := + dp.foldl pscatterStep (Array.replicate (base - 1) pid) + +@[simp] theorem size_pscatterStep (a : Array PVes) (p : ℕ × PVes) : + (pscatterStep a p).size = a.size := by + unfold pscatterStep; split <;> simp + +@[simp] theorem size_foldl_pscatterStep (dp : List (ℕ × PVes)) (a : Array PVes) : + (dp.foldl pscatterStep a).size = a.size := by + induction dp generalizing a with + | nil => rfl + | cons p dp ih => rw [List.foldl_cons, ih, size_pscatterStep] + +/-- A single projective scatter step through `getElem?`: slot `k` gains a right-`padd` of `x` +exactly when the digit is `k + 1`. -/ +theorem getElem?_pscatterStep (a : Array PVes) (d : ℕ) (x : PVes) (k : ℕ) : + (pscatterStep a (d, x))[k]? + = if d = k + 1 then a[k]?.map (fun v => padd v x) else a[k]? := by + unfold pscatterStep + by_cases h0 : d = 0 + · subst h0 + rw [if_pos rfl, if_neg (by omega)] + · rw [if_neg h0, Array.getElem?_modify] + by_cases hdk : d = k + 1 + · rw [if_pos (by omega), if_pos hdk] + · rw [if_neg (by omega), if_neg hdk] + +/-- **The scatter fold invariant, transported through `toAffine`.** Starting from a valid slot +value `v`, slot `k` ends valid and reads — in `G` — as `toAffine v` plus the affine bucket sum +of the processed digit-tagged list. -/ +private theorem foldl_pscatterStep_spec (dp : List (ℕ × PVes)) (hdp : ∀ p ∈ dp, Valid p.2) : + ∀ (a : Array PVes) (k : ℕ) (v : PVes), a[k]? = some v → Valid v → + ∃ w, (dp.foldl pscatterStep a)[k]? = some w ∧ Valid w ∧ + toAffine w = toAffine v + + Msm.bucketOf (dp.map fun t => (t.1, toAffine t.2)) (k + 1) := by + induction dp with + | nil => + intro a k v hv hval + exact ⟨v, hv, hval, by simp⟩ + | cons p dp ih => + intro a k v hv hval + obtain ⟨d, x⟩ := p + have hx : Valid x := hdp (d, x) (by simp) + have hdp' : ∀ q ∈ dp, Valid q.2 := fun q hq => hdp q (by simp [hq]) + rw [List.foldl_cons] + by_cases hdk : d = k + 1 + · have hv' : (pscatterStep a (d, x))[k]? = some (padd v x) := by + rw [getElem?_pscatterStep, if_pos hdk, hv, Option.map_some] + obtain ⟨w, hw, hwval, hweq⟩ := + ih hdp' (pscatterStep a (d, x)) k (padd v x) hv' (valid_padd hval hx) + refine ⟨w, hw, hwval, ?_⟩ + rw [hweq, toAffine_padd hval hx, List.map_cons, Msm.bucketOf_cons, if_pos hdk, _root_.add_assoc] + · have hv' : (pscatterStep a (d, x))[k]? = some v := by + rw [getElem?_pscatterStep, if_neg hdk, hv] + obtain ⟨w, hw, hwval, hweq⟩ := ih hdp' (pscatterStep a (d, x)) k v hv' hval + refine ⟨w, hw, hwval, ?_⟩ + rw [hweq, List.map_cons, Msm.bucketOf_cons, if_neg hdk, _root_.zero_add] + +/-- **The scattered buckets are the affine bucket list, slotwise.** Every slot is valid, and +`toAffine` maps the slot list to exactly the bucket list `Msm.windowValue` folds over. -/ +theorem pbucketScatter_spec (base : ℕ) (dp : List (ℕ × PVes)) (hdp : ∀ p ∈ dp, Valid p.2) : + (∀ P ∈ (pbucketScatter base dp).toList, Valid P) + ∧ (pbucketScatter base dp).toList.map toAffine + = (List.range (base - 1)).map fun k => + Msm.bucketOf (dp.map fun t => (t.1, toAffine t.2)) (k + 1) := by + have hsize : (pbucketScatter base dp).size = base - 1 := by + rw [pbucketScatter, size_foldl_pscatterStep, Array.size_replicate] + have hslot : ∀ k, k < base - 1 → + ∃ w, (pbucketScatter base dp)[k]? = some w ∧ Valid w ∧ + toAffine w = Msm.bucketOf (dp.map fun t => (t.1, toAffine t.2)) (k + 1) := by + intro k hk + have hinit : (Array.replicate (base - 1) pid)[k]? = some pid := by + rw [Array.getElem?_replicate, if_pos hk] + obtain ⟨w, hw, hwval, hweq⟩ := + foldl_pscatterStep_spec dp hdp (Array.replicate (base - 1) pid) k pid hinit valid_pid + exact ⟨w, hw, hwval, by rw [hweq, toAffine_pid, _root_.zero_add]⟩ + constructor + · intro P hP + obtain ⟨k, hk, hPk⟩ := List.mem_iff_getElem.mp hP + have hk' : k < base - 1 := by + rw [Array.length_toList, hsize] at hk + exact hk + obtain ⟨w, hw, hwval, -⟩ := hslot k hk' + have hPk? : (pbucketScatter base dp)[k]? = some P := by + rw [← Array.getElem?_toList, List.getElem?_eq_getElem hk, hPk] + rw [hPk?] at hw + obtain rfl : P = w := Option.some.inj hw + exact hwval + · apply List.ext_getElem? + intro k + rw [List.getElem?_map, List.getElem?_map, Array.getElem?_toList] + by_cases hk : k < base - 1 + · obtain ⟨w, hw, -, hweq⟩ := hslot k hk + rw [hw, List.getElem?_range hk, Option.map_some, Option.map_some, hweq] + · rw [Array.getElem?_eq_none (by rw [hsize]; omega), + List.getElem?_eq_none (by rw [List.length_range]; omega)] + rfl + +/-- The projective window value via the single-pass scatter (mirror of +`Msm.windowValueFast`). -/ +def pwindowValueFast (base i : ℕ) (pterms : List (ℕ × PVes)) : PVes := + (List.foldr paccStep (pid, pid) (pbucketScatter base (pdpOf base i pterms)).toList).2 + +/-- The scatter-bucketed projective window value is valid and matches `Msm.windowValue` of the +`toAffine`-mapped terms — the scatter twin of `pwindowValue_spec`. -/ +theorem pwindowValueFast_spec (base i : ℕ) (pterms : List (ℕ × PVes)) + (h : ∀ p ∈ pterms, Valid p.2) : + Valid (pwindowValueFast base i pterms) + ∧ toAffine (pwindowValueFast base i pterms) + = Msm.windowValue base i (pterms.map fun t => (t.1, toAffine t.2)) := by + have hdp : ∀ p ∈ pdpOf base i pterms, Valid p.2 := by + intro p hp + rw [pdpOf, List.mem_map] at hp + obtain ⟨t, ht, rfl⟩ := hp + exact h t ht + obtain ⟨hbval, hbmap⟩ := pbucketScatter_spec base (pdpOf base i pterms) hdp + have hdpof : (pdpOf base i pterms).map (fun t => (t.1, toAffine t.2)) + = Msm.dpOf base i (pterms.map fun t => (t.1, toAffine t.2)) := by + simp only [pdpOf, Msm.dpOf, List.map_map, Function.comp_def] + rw [hdpof] at hbmap + obtain ⟨-, hv2, heq⟩ := foldr_paccStep_spec _ hbval + simp only [pwindowValueFast] + refine ⟨hv2, ?_⟩ + rw [Msm.windowValue, ← hbmap, ← heq] + +/-! ## The scatter-bucketed projective Pippenger MSM -/ + +/-- **Windowed Pippenger MSM in projective coordinates with single-pass Array bucketing** — +the fast serial form of `pippengerProj`. Proven equal to the affine `Msm.pippenger` +(`pippengerProjScatter_eq`). -/ +def pippengerProjScatter (c : ℕ) (terms : List (ℕ × Projective.G)) : Projective.G := + toAffine (phornerList (2 ^ c) + ((List.range (Msm.numWindows c terms)).map fun i => + pwindowValueFast (2 ^ c) i (terms.map fun t => (t.1, ofAffine t.2)))) + +/-- **The scatter-bucketed projective Pippenger equals the affine Pippenger** — same transport +as `pippengerProj_eq`, window values via `pwindowValueFast_spec`. -/ +theorem pippengerProjScatter_eq (c : ℕ) (terms : List (ℕ × Projective.G)) : + pippengerProjScatter c terms = Msm.pippenger c terms := by + rw [pippengerProjScatter] + set pterms := terms.map (fun t => (t.1, ofAffine t.2)) with hpterms + set windows := (List.range (Msm.numWindows c terms)).map fun i => + pwindowValueFast (2 ^ c) i pterms + with hwin + have hval : ∀ p ∈ pterms, Valid p.2 := by + intro p hp + rw [hpterms, List.mem_map] at hp + obtain ⟨t, _, rfl⟩ := hp + exact valid_ofAffine t.2 + have hf : pterms.map (fun t => (t.1, toAffine t.2)) = terms := by + rw [hpterms, List.map_map] + conv_rhs => rw [← List.map_id terms] + apply List.map_congr_left + intro t _ + simp only [Function.comp_apply, toAffine_ofAffine, id_eq] + have hwval : ∀ P ∈ windows, Valid P := by + intro P hP + rw [hwin, List.mem_map] at hP + obtain ⟨i, _, rfl⟩ := hP + exact (pwindowValueFast_spec (2 ^ c) i pterms hval).1 + rw [(phornerList_spec (2 ^ c) windows hwval).2, Msm.pippenger] + congr 1 + rw [hwin, List.map_map] + apply List.map_congr_left + intro i _ + rw [Function.comp_apply, (pwindowValueFast_spec (2 ^ c) i pterms hval).2, hf] + +/-- The scatter-bucketed projective Pippenger equals the naive MSM. -/ +theorem pippengerProjScatter_eq_msm (c : ℕ) (hc : 0 < c) (terms : List (ℕ × Projective.G)) : + pippengerProjScatter c terms = (terms.map fun t => t.1 • t.2).sum := by + rw [pippengerProjScatter_eq, Msm.pippenger_eq_msm c hc] + +/-! ## Windows-parallel projective Pippenger + +The window values are mutually independent; evaluating them through the proven `Msm.parMap` +changes only the evaluation strategy, so the equality is purely `Msm.parMap_eq_map`. Measured +(interpreted, `n = 2048`, `c = 8`, 12 cores): `4.3 s → 0.8 s` per MSM on top of the scatter +port. -/ + +/-- **Windows-parallel scatter-bucketed projective Pippenger**: `pippengerProjScatter` with the +window values evaluated as parallel tasks. -/ +def pippengerProjScatterPar (c : ℕ) (terms : List (ℕ × Projective.G)) : Projective.G := + toAffine (phornerList (2 ^ c) + (Msm.parMap (fun i => + pwindowValueFast (2 ^ c) i (terms.map fun t => (t.1, ofAffine t.2))) + (List.range (Msm.numWindows c terms)))) + +/-- The windows-parallel projective Pippenger is the sequential one: `parMap` is `map`. -/ +theorem pippengerProjScatterPar_eq (c : ℕ) (terms : List (ℕ × Projective.G)) : + pippengerProjScatterPar c terms = pippengerProjScatter c terms := by + rw [pippengerProjScatterPar, Msm.parMap_eq_map, pippengerProjScatter] + +/-- The windows-parallel projective Pippenger equals the naive MSM. -/ +theorem pippengerProjScatterPar_eq_msm (c : ℕ) (hc : 0 < c) + (terms : List (ℕ × Projective.G)) : + pippengerProjScatterPar c terms = (terms.map fun t => t.1 • t.2).sum := by + rw [pippengerProjScatterPar_eq, pippengerProjScatter_eq_msm c hc] + +/-! ## The `commit_lagrange` wrapper -/ + +/-- The fast `commit_lagrange` run in projective coordinates: projective windowed Pippenger over the +`(coeffᵢ.val, basisᵢ)` terms (single final inversion), plus the blind. The terms are zipped +identically to `Fast.Msm.commitLagrangeFastWith` so the wrapper equality reuses `zip_terms_eq`. -/ +def commitLagrangeProjWith (c : ℕ) (blind : Projective.G) (basis : List Projective.G) + (coeffs : List Fp) : Projective.G := + pippengerProj c + ((coeffs.zip (basis ++ List.replicate (coeffs.length - basis.length) 0)).map + fun t => (t.1.val, t.2)) + blind + +/-- **The projective commitment equals the naive `commit_lagrange` spec** +(`Fast.Msm.commitLagrangeSpec`, verbatim from `Keygen.commitLagrangeWith`). -/ +theorem commitLagrangeProjWith_eq (c : ℕ) (hc : 0 < c) + (blind : Projective.G) (basis : List Projective.G) (coeffs : List Fp) : + commitLagrangeProjWith c blind basis coeffs = Msm.commitLagrangeSpec blind basis coeffs := by + unfold commitLagrangeProjWith Msm.commitLagrangeSpec + rw [pippengerProj_eq, Msm.zip_terms_eq, Msm.pippenger_eq_msm c hc, List.map_map] + rfl + +/-- The fast `commit_lagrange` with scatter-bucketed projective Pippenger — the drop-in +committer for the certificate's per-column pass (serial per column; the 44-column outer +`parMap` already saturates the cores). -/ +def commitLagrangeProjScatterWith (c : ℕ) (blind : Projective.G) (basis : List Projective.G) + (coeffs : List Fp) : Projective.G := + pippengerProjScatter c + ((coeffs.zip (basis ++ List.replicate (coeffs.length - basis.length) 0)).map + fun t => (t.1.val, t.2)) + blind + +/-- **The scatter-bucketed projective commitment equals the naive `commit_lagrange` spec** +(`Fast.Msm.commitLagrangeSpec`). -/ +theorem commitLagrangeProjScatterWith_eq (c : ℕ) (hc : 0 < c) + (blind : Projective.G) (basis : List Projective.G) (coeffs : List Fp) : + commitLagrangeProjScatterWith c blind basis coeffs + = Msm.commitLagrangeSpec blind basis coeffs := by + unfold commitLagrangeProjScatterWith Msm.commitLagrangeSpec + rw [pippengerProjScatter_eq, Msm.zip_terms_eq, Msm.pippenger_eq_msm c hc, List.map_map] + rfl + +end CompElliptic.Curves.Pasta.Fast.MsmProj diff --git a/CompElliptic/Curves/Pasta/Fast/Projective.lean b/CompElliptic/Curves/Pasta/Fast/Projective.lean new file mode 100644 index 0000000..5d00624 --- /dev/null +++ b/CompElliptic/Curves/Pasta/Fast/Projective.lean @@ -0,0 +1,633 @@ +/- +Copyright (c) 2026 CompElliptic Contributors. All rights reserved. +Released under the Apache License, Version 2.0, or the MIT license, at your option, +as described in the files LICENSE-APACHE and LICENSE-MIT. +Authors: Gregor Mitscha-Baude +-/ +import CompElliptic.Curves.Pasta + +/-! +# Projective (Renes–Costello–Batina) point arithmetic for Vesta + +The affine group law pays a field inversion per addition. This module replaces it, for the hot +`n • point` loop, by the complete projective addition formulas of Renes, Costello and Batina +(EUROCRYPT 2016; the EFD `add-2015-rcb` sequence at `a = 0`): single branchless formulas in +`(X : Y : Z)` valid for *all* input pairs, with the inversion paid once in `toAffine`. Their +completeness needs the curve to have no 2-torsion, which Vesta's odd prime order gives +(`Vesta.no_onCurve_y_zero`); that is what discharges `Z₃ ≠ 0` outside the genuine identity cases. + +Specialized to `a = 0`, `b = 5`, `b3 = 15`, the closed forms are + +* `X₃ = X₁Y₁Y₂² − 15X₁Y₁Z₂² − 30X₁Z₁Y₂Z₂ + Y₁²X₂Y₂ − 15Z₁²X₂Y₂ − 30Y₁Z₁X₂Z₂` +* `Y₃ = Y₁²Y₂² + 45X₁²X₂Z₂ + 45X₁Z₁X₂² − 225Z₁²Z₂²` +* `Z₃ = Y₁²Y₂Z₂ + Y₁Z₁Y₂² + 3X₁²X₂Y₂ + 3X₁Y₁X₂² + 15Y₁Z₁Z₂² + 15Z₁²Y₂Z₂` + +The `Fast` interfaces are provisional: they are not guaranteed to remain public, and +may be folded into the existing API or otherwise changed incompatibly. +-/ + +open CompElliptic +open CompElliptic.CurveForms.ShortWeierstrass +open CompElliptic.Curves.Pasta + +namespace CompElliptic.Curves.Pasta.Fast.Projective + +/-- The Vesta base field `𝔽_q` (= `PallasScalarField`), over which the Vesta curve is defined. -/ +abbrev Fq := CompElliptic.Fields.Pasta.VestaBaseField + +/-- The downstream affine group: on-curve points of Vesta with the complete affine group law. -/ +abbrev G := SWPoint Vesta.curve + +/-- A projective point in `(X : Y : Z)` coordinates over `𝔽_q`. -/ +structure PVes where + X : Fq + Y : Fq + Z : Fq +deriving DecidableEq + +namespace PVes + +/-- Renes–Costello–Batina complete addition (`add-2015-rcb`, `a = 0`, `b3 = 15`). -/ +def padd (P Q : PVes) : PVes where + X := P.X*P.Y*Q.Y^2 - 15*P.X*P.Y*Q.Z^2 - 30*P.X*P.Z*Q.Y*Q.Z + + P.Y^2*Q.X*Q.Y - 15*P.Z^2*Q.X*Q.Y - 30*P.Y*P.Z*Q.X*Q.Z + Y := P.Y^2*Q.Y^2 + 45*P.X^2*Q.X*Q.Z + 45*P.X*P.Z*Q.X^2 - 225*P.Z^2*Q.Z^2 + Z := P.Y^2*Q.Y*Q.Z + P.Y*P.Z*Q.Y^2 + 3*P.X^2*Q.X*Q.Y + 3*P.X*P.Y*Q.X^2 + + 15*P.Y*P.Z*Q.Z^2 + 15*P.Z^2*Q.Y*Q.Z + +/-- Scalar (representative) rescaling `(X : Y : Z) ↦ (uX : uY : uZ)`. -/ +def smul (u : Fq) (P : PVes) : PVes := ⟨u*P.X, u*P.Y, u*P.Z⟩ + +/-- The projective identity `𝒪 = (0 : 1 : 0)`. -/ +def pid : PVes := ⟨0, 1, 0⟩ + +/-- The homogeneous projective curve equation `Y²Z = X³ + 5Z³` (Vesta, `a = 0`, `b = 5`). -/ +def OnCurveP (P : PVes) : Prop := P.Y^2*P.Z = P.X^3 + 5*P.Z^3 + +/-- `OnCurveP` unfolds definitionally to a field equation, hence is decidable. This is what keeps +`toAffine` — and everything downstream of it, `smulFast` in particular — computable. -/ +instance (P : PVes) : Decidable (OnCurveP P) := + inferInstanceAs (Decidable (P.Y^2*P.Z = P.X^3 + 5*P.Z^3)) + +/-- A representable projective point: on the projective curve and not the zero vector. -/ +def Valid (P : PVes) : Prop := OnCurveP P ∧ (P.X ≠ 0 ∨ P.Y ≠ 0 ∨ P.Z ≠ 0) + +/-- Affine interpretation as a raw coordinate pair: `Z ≠ 0 ↦ (X/Z, Y/Z)`, `Z = 0 ↦ 𝒪 = (0,0)`. -/ +def aff (P : PVes) : Fq × Fq := if P.Z = 0 then (0, 0) else (P.X / P.Z, P.Y / P.Z) + +/-! ## `ring`/`linear_combination` certificates + +Each identity below is a polynomial identity over `𝔽_q`, valid modulo the two curve equations +`e1 : y1² = x1³ + 5` and `e2 : y2² = x2³ + 5`. The cofactors were computed by an exact-rational +Gröbner/linear-algebra search over `ℚ[x1,y1,x2,y2]` (offline) and are checked here by `ring` +inside `linear_combination`. The `padd`-of-`Z=1`-points evaluation is spelled by `simp only [padd]` +which reduces the coordinate projections. -/ + +variable {x1 y1 x2 y2 : Fq} + +/-- Reading off the coordinates of `padd` on `Z = 1` representatives. -/ +@[simp] theorem padd_z1_X (x1 y1 x2 y2 : Fq) : + (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩).X + = y1^2*x2*y2 + x1*y1*y2^2 - 15*x2*y2 - 30*y1*x2 - 30*x1*y2 - 15*x1*y1 := by + simp only [padd]; ring + +@[simp] theorem padd_z1_Y (x1 y1 x2 y2 : Fq) : + (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩).Y + = y1^2*y2^2 + 45*x1*x2^2 + 45*x1^2*x2 - 225 := by + simp only [padd]; ring + +@[simp] theorem padd_z1_Z (x1 y1 x2 y2 : Fq) : + (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩).Z + = 3*x1*y1*x2^2 + 3*x1^2*x2*y2 + y1*y2^2 + y1^2*y2 + 15*y2 + 15*y1 := by + simp only [padd]; ring + +/-- `X₃·(x₂−x₁)² = x3num·Z₃`, the numerator match for the affine `x`-coordinate (distinct `x`). -/ +theorem keyx_dist (e1 : y1^2 = x1^3 + 5) (e2 : y2^2 = x2^3 + 5) : + (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩).X * (x2 - x1)^2 + = ((y2 - y1)^2 - (x1 + x2)*(x2 - x1)^2) * (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩).Z := by + simp only [padd_z1_X, padd_z1_Z] + linear_combination + ((2)*x2^3*y2 + (3)*x1*x2^2*y2 + (-3)*x1*y1*x2^2 + (-3)*x1^2*x2*y2 + y2^3 + y1*y2^2 - y1^2*y2 + (10)*y2 + (-15)*y1) * e1 + + ((-3)*x1*y1*x2^2 + (-3)*x1^2*x2*y2 + (3)*x1^2*y1*x2 + x1^3*y2 + (3)*x1^3*y1 - y1*y2^2 + (-10)*y2 + (15)*y1) * e2 + +/-- `Y₃·(x₂−x₁)³ = y3num·Z₃`, the numerator match for the affine `y`-coordinate (distinct `x`). -/ +theorem keyy_dist (e1 : y1^2 = x1^3 + 5) (e2 : y2^2 = x2^3 + 5) : + (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩).Y * (x2 - x1)^3 + = ((y2 - y1)*(x1*(x2 - x1)^2 - ((y2 - y1)^2 - (x1 + x2)*(x2 - x1)^2)) - y1*(x2 - x1)^3) + * (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩).Z := by + simp only [padd_z1_Y, padd_z1_Z] + linear_combination + ((6)*x1*x2^5 + (-9)*x1^2*x2^4 + (2)*x2^3*y2^2 + (2)*y1*x2^3*y2 + (-15)*x1*x2^2*y2^2 + (6)*x1*y1*x2^2*y2 + (-3)*x1*y1^2*x2^2 + (15)*x1^2*x2*y2^2 + (-3)*x1^2*y1*x2*y2 + (-2)*y2^4 + (2)*y1^2*y2^2 - y1^3*y2 + (30)*x2^3 + (-60)*x1*x2^2 + (10)*y2^2 + (25)*y1*y2 + (-15)*y1^2 + (-75)) * e1 + + ((-6)*x1^4*x2^2 + (9)*x1^5*x2 + (3)*x1*y1*x2^2*y2 + (3)*x1^2*x2*y2^2 + (-6)*x1^2*y1*x2*y2 + (-2)*x1^3*y2^2 + (-2)*x1^3*y1*y2 + y1*y2^3 + (-75)*x1*x2^2 + (135)*x1^2*x2 + (-30)*x1^3 + (5)*y2^2 + (-25)*y1*y2 + (75)) * e2 + +/-- Curve preservation on `Z = 1` representatives: `padd` lands on the projective curve. -/ +theorem oncurveP_padd_z1 (e1 : y1^2 = x1^3 + 5) (e2 : y2^2 = x2^3 + 5) : + OnCurveP (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩) := by + simp only [OnCurveP, padd_z1_X, padd_z1_Y, padd_z1_Z] + linear_combination + (-x1^6*x2^3*y2^3 - x1^3*y1^2*x2^3*y2^3 + x1^6*y2^5 - y1^4*x2^3*y2^3 + (-135)*x1^3*y1*x2^6 + x1^3*y1^2*y2^5 + (-405)*x1^4*x2^5*y2 + (-135)*x1^5*x2^4*y2 + y1^3*y2^6 + y1^4*y2^5 + (135)*x1^2*y1*x2^4*y2^2 + (-135)*x1^2*y1^2*x2^4*y2 + (35)*x1^3*x2^3*y2^3 + (90)*x1^3*y1*x2^3*y2^2 + (405)*x1^4*x2^2*y2^3 + (135)*x1^5*x2*y2^3 + (-5)*x1^6*y2^3 + (40)*y1^2*x2^3*y2^3 + (90)*y1^3*x2^3*y2^2 + (135)*x1*y1*x2^2*y2^4 + (270)*x1*y1^2*x2^2*y2^3 + (270)*x1^2*y1*x2*y2^4 + (135)*x1^2*y1^2*x2*y2^3 + (100)*x1^3*y2^5 + (45)*x1^3*y1*y2^4 + (-5)*x1^3*y1^2*y2^3 + (5)*y1^2*y2^5 + (-5)*y1^4*y2^3 + (-675)*x1^2*x2^4*y2 + (-2025)*x1^2*y1*x2^4 + (-2700)*x1^3*x2^3*y2 + (-2025)*x1^4*x2^2*y2 + (-675)*x1^5*x2*y2 + (-475)*x2^3*y2^3 + (-2250)*y1*x2^3*y2^2 + (-2700)*y1^2*x2^3*y2 + (-4050)*x1*x2^2*y2^3 + (-12150)*x1*y1*x2^2*y2^2 + (-4050)*x1*y1^2*x2^2*y2 + (-11475)*x1^2*x2*y2^3 + (-5400)*x1^2*y1*x2*y2^2 + (-675)*x1^2*y1^2*x2*y2 + (-3875)*x1^3*y2^3 + (-900)*x1^3*y1*y2^2 + (-200)*y2^5 + (-1125)*y1*y2^4 + (-1150)*y1^2*y2^3 + (-225)*y1^3*y2^2 + (27000)*x2^3*y2 + (27000)*y1*x2^3 + (60750)*x1*x2^2*y2 + (30375)*x1*y1*x2^2 + (57375)*x1^2*x2*y2 + (20250)*x1^2*y1*x2 + (16875)*x1^3*y2 + (3375)*x1^3*y1 + (-22625)*y2^3 + (-18000)*y1*y2^2 + (-3375)*y1^2*y2 + (-16875)*y2 + (-16875)*y1) * e1 + + (x1^9*y2^3 + (135)*x1^6*y1*x2^3 + (405)*x1^7*x2^2*y2 + (135)*x1^8*x2*y2 + (270)*x1^5*y1*x2*y2^2 + (105)*x1^6*y2^3 + (45)*x1^6*y1*y2^2 + (-5400)*x1^3*y1*x2^3 + (-4050)*x1^4*x2^2*y2 + (-12150)*x1^4*y1*x2^2 + (-10800)*x1^5*x2*y2 + (-4050)*x1^5*y1*x2 + (-3375)*x1^6*y2 + (-675)*x1^6*y1 + (-2700)*x1^2*y1*x2*y2^2 + (300)*x1^3*y2^3 + (-3600)*x1^3*y1*y2^2 + (-27000)*x1^2*x2*y2 + (40500)*x1^2*y1*x2 + (-13500)*x1^3*y2 + (-1000)*y2^3 + (-9000)*y1*y2^2 + (-135000)*y2 + (-135000)*y1) * e2 + +/-- **Distinct-`x` completeness (2-torsion-free).** If `Z₃ = 0` for two on-curve `Z = 1` points +with distinct `x`, then `w = wnum/wden` (the `x`-coordinate of `P − Q`) is a cube root of `−5`, +i.e. `(w, 0)` is a 2-torsion point — impossible on Vesta. Hence `Z₃ ≠ 0`. -/ +theorem Z_ne_zero_dist (e1 : y1^2 = x1^3 + 5) (e2 : y2^2 = x2^3 + 5) (hne : x1 ≠ x2) : + (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩).Z ≠ 0 := by + intro hZ + have hd : x2 - x1 ≠ 0 := sub_ne_zero.mpr (Ne.symm hne) + set wnum : Fq := (y1 + y2)^2 - (x1 + x2)*(x2 - x1)^2 with hwnumE + set wden : Fq := (x2 - x1)^2 with hwdenE + have hwdenne : wden ≠ 0 := by rw [hwdenE]; exact pow_ne_zero 2 hd + have hZpoly : (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩).Z + = 3*x1*y1*x2^2 + 3*x1^2*x2*y2 + y1*y2^2 + y1^2*y2 + 15*y2 + 15*y1 := padd_z1_Z x1 y1 x2 y2 + rw [hZpoly] at hZ + -- wnum^3 = -5 * wden^3 on the exceptional locus + have hcube : wnum^3 = -5 * wden^3 := by + rw [hwnumE, hwdenE] + linear_combination + (-x2^3*y2 + (-2)*y1*x2^3 + (3)*x1*y1*x2^2 + (3)*x1^2*x2*y2 + (-2)*x1^3*y2 + (-4)*x1^3*y1 + y2^3 + (3)*y1*y2^2 + (3)*y1^2*y2 + (4)*y1^3 + (-15)*y1) * hZ + + ((8)*x2^6 + (-6)*x1*x2^5 + (-6)*x1^2*x2^4 + (8)*x1^3*x2^3 + (-3)*x1^5*x2 + x1^6 + (-12)*x2^3*y2^2 + (-6)*y1*x2^3*y2 + (-4)*y1^2*x2^3 + (12)*x1*x2^2*y2^2 + (-9)*x1*y1^2*x2^2 + (3)*x1^2*y1^2*x2 + (-3)*x1^3*y2^2 + (-6)*x1^3*y1*y2 + (-2)*x1^3*y1^2 + (3)*y2^4 + (10)*y1*y2^3 + (9)*y1^2*y2^2 + (2)*y1^3*y2 + y1^4 + (60)*x2^3 + (-75)*x1*x2^2 + (45)*x1^2*x2 + (-10)*x1^3 + (-15)*y2^2 + (-60)*y1*y2 + (-60)*y1^2 + (50)) * e1 + + (x2^6 + (-3)*x1*x2^5 + (-2)*x2^3*y2^2 + (-6)*y1*x2^3*y2 + (5)*y1^2*x2^3 + (3)*x1*x2^2*y2^2 + (9)*x1*y1*x2^2*y2 + (-6)*x1*y1^2*x2^2 + (6)*x1^2*y1^2*x2 + y2^4 + (5)*y1*y2^3 + (8)*y1^2*y2^2 + (4)*y1^3*y2 - y1^4 + (-50)*x2^3 + (75)*x1*x2^2 + (-45)*x1^2*x2 + (5)*y2^2 + (15)*y1*y2 + (25)*y1^2 + (-50)) * e2 + -- so (wnum/wden) is a cube root of -5, contradicting no 2-torsion + have hw : (wnum / wden)^3 = -5 := by + rw [div_pow, hcube, mul_div_assoc, div_self (pow_ne_zero 3 hwdenne), mul_one] + exact Vesta.no_onCurve_y_zero (wnum / wden) (by + show (0 : Fq)^2 = (wnum / wden)^3 + Vesta.a * (wnum / wden) + Vesta.b + rw [hw]; simp [Vesta.a, Vesta.b]) + +/-- **Doubling `Z`-coordinate.** On `Z = 1` representatives with equal points, `Z₃ = 8y³`. -/ +theorem Z_doubling (e1 : y1^2 = x1^3 + 5) : + (padd ⟨x1, y1, 1⟩ ⟨x1, y1, 1⟩).Z = 8 * y1^3 := by + simp only [padd_z1_Z] + linear_combination ((-6)*y1) * e1 + +/-- `X₃·4y² = x3numd·Z₃`, doubling `x`-coordinate numerator match. -/ +theorem keyx_dbl (e1 : y1^2 = x1^3 + 5) : + (padd ⟨x1, y1, 1⟩ ⟨x1, y1, 1⟩).X * (4*y1^2) + = (9*x1^4 - 8*x1*y1^2) * (padd ⟨x1, y1, 1⟩ ⟨x1, y1, 1⟩).Z := by + simp only [padd_z1_X, padd_z1_Z] + linear_combination ((54)*x1^4*y1 + (24)*x1*y1^3) * e1 + +/-- `Y₃·8y³ = y3numd·Z₃`, doubling `y`-coordinate numerator match. -/ +theorem keyy_dbl (e1 : y1^2 = x1^3 + 5) : + (padd ⟨x1, y1, 1⟩ ⟨x1, y1, 1⟩).Y * (8*y1^3) + = (3*x1^2*(12*x1*y1^2 - 9*x1^4) - 8*y1^4) * (padd ⟨x1, y1, 1⟩ ⟨x1, y1, 1⟩).Z := by + simp only [padd_z1_Y, padd_z1_Z] + linear_combination ((-162)*x1^6*y1 + (24)*y1^5 + (360)*y1^3) * e1 + +/-- **Inverse-case identity representative is nonzero (2-torsion-free).** `padd P (−P) = (0 : Y₃ : 0)` +with `Y₃ ≠ 0`: if `Y₃ = 0` then `w = wnum2/wden2` (the `x`-coordinate of `2P`) is a cube root of +`−5`, impossible on Vesta. -/ +theorem Y_ne_zero_inv (e1 : y1^2 = x1^3 + 5) (hy : y1 ≠ 0) : + (padd ⟨x1, y1, 1⟩ ⟨x1, -y1, 1⟩).Y ≠ 0 := by + intro hY + have hYpoly : (padd ⟨x1, y1, 1⟩ ⟨x1, -y1, 1⟩).Y = y1^4 + 90*x1^3 - 225 := by + simp only [padd]; ring + rw [hYpoly] at hY + set wnum2 : Fq := 9*x1^4 - 8*x1*y1^2 with hwnum2 + set wden2 : Fq := 4*y1^2 with hwden2E + have hwden2ne : wden2 ≠ 0 := by rw [hwden2E]; exact mul_ne_zero (by decide) (pow_ne_zero 2 hy) + have hcube : wnum2^3 = -5 * wden2^3 := by + rw [hwnum2, hwden2E] + linear_combination + (x1^6 + (100)*x1^3 - 200) * hY + + ((-729)*x1^9 + (1215)*x1^6*y1^2 + (-512)*x1^3*y1^4 + (3735)*x1^6 + (-2340)*x1^3*y1^2 + (320)*y1^4 + (-9900)*x1^3 + (1800)*y1^2 + (9000)) * e1 + have hw : (wnum2 / wden2)^3 = -5 := by + rw [div_pow, hcube, mul_div_assoc, div_self (pow_ne_zero 3 hwden2ne), mul_one] + exact Vesta.no_onCurve_y_zero (wnum2 / wden2) (by + show (0 : Fq)^2 = (wnum2 / wden2)^3 + Vesta.a * (wnum2 / wden2) + Vesta.b + rw [hw]; simp [Vesta.a, Vesta.b]) + +/-- Inverse `X`-coordinate vanishes: `padd P (−P) = (0 : Y₃ : 0)` (pure identity). -/ +theorem X_zero_inv (x1 y1 : Fq) : (padd ⟨x1, y1, 1⟩ ⟨x1, -y1, 1⟩).X = 0 := by + simp only [padd]; ring + +/-- Inverse `Z`-coordinate vanishes (pure identity). -/ +theorem Z_zero_inv (x1 y1 : Fq) : (padd ⟨x1, y1, 1⟩ ⟨x1, -y1, 1⟩).Z = 0 := by + simp only [padd]; ring + +/-! ## Scaling (homogeneity) and affine interpretation -/ + +/-- `padd` is homogeneous of bidegree `(2,2)`: rescaling the inputs rescales the output. -/ +theorem smul_padd (u v : Fq) (P Q : PVes) : + padd (smul u P) (smul v Q) = smul (u^2*v^2) (padd P Q) := by + simp only [padd, smul, PVes.mk.injEq]; refine ⟨?_, ?_, ?_⟩ <;> ring + +/-- `aff` is invariant under nonzero rescaling. -/ +theorem aff_smul (u : Fq) (hu : u ≠ 0) (P : PVes) : aff (smul u P) = aff P := by + rcases eq_or_ne P.Z 0 with h | h + · simp [aff, smul, h] + · have huz : u * P.Z ≠ 0 := mul_ne_zero hu h + simp only [aff, smul, if_neg h, if_neg huz] + rw [Prod.mk.injEq] + exact ⟨mul_div_mul_left _ _ hu, mul_div_mul_left _ _ hu⟩ + +/-- `OnCurveP` is preserved by rescaling. -/ +theorem oncurveP_smul (u : Fq) {P : PVes} (h : OnCurveP P) : OnCurveP (smul u P) := by + simp only [OnCurveP, smul] at h ⊢; linear_combination (u^3) * h + +/-- The nonzero-vector condition is preserved by nonzero rescaling. -/ +theorem nezVec_smul (u : Fq) (hu : u ≠ 0) {P : PVes} + (h : P.X ≠ 0 ∨ P.Y ≠ 0 ∨ P.Z ≠ 0) : + (smul u P).X ≠ 0 ∨ (smul u P).Y ≠ 0 ∨ (smul u P).Z ≠ 0 := by + simp only [smul] + rcases h with h | h | h + · exact Or.inl (mul_ne_zero hu h) + · exact Or.inr (Or.inl (mul_ne_zero hu h)) + · exact Or.inr (Or.inr (mul_ne_zero hu h)) + +@[simp] theorem aff_z1 (x y : Fq) : aff ⟨x, y, 1⟩ = (x, y) := by + simp [aff] + +/-- `P = P.Z • (X/Z : Y/Z : 1)` when `Z ≠ 0`: every finite point is a rescaled `Z = 1` rep. -/ +theorem eq_smul_normalize {P : PVes} (h : P.Z ≠ 0) : + P = smul P.Z ⟨P.X / P.Z, P.Y / P.Z, 1⟩ := by + obtain ⟨X, Y, Z⟩ := P + simp only [smul, PVes.mk.injEq] + refine ⟨?_, ?_, ?_⟩ <;> field_simp + +/-! ## Unfolding the affine group law on `Z = 1` reps -/ + +/-- `(x₁,y₁)` and `(x₂,y₂)` on the curve are `≠ 𝒪 = (0,0)`. -/ +theorem onCurve_ne_zero {x y : Fq} (h : OnCurve 0 5 (x, y)) : (x, y) ≠ ((0 : Fq), (0 : Fq)) := by + intro he + rw [Prod.mk.injEq] at he + obtain ⟨hx, hy⟩ := he; subst hx; subst hy + exact not_onCurve_zero (show (5 : Fq) ≠ 0 by decide) h + +/-- The affine sum in the distinct-`x` branch, with cleared denominators. -/ +theorem add_dist {x1 y1 x2 y2 : Fq} (h1 : OnCurve 0 5 (x1, y1)) (h2 : OnCurve 0 5 (x2, y2)) + (hne : x1 ≠ x2) : + add 0 (x1, y1) (x2, y2) + = (((y2 - y1)^2 - (x1 + x2)*(x2 - x1)^2) / (x2 - x1)^2, + ((y2 - y1)*(x1*(x2 - x1)^2 - ((y2 - y1)^2 - (x1 + x2)*(x2 - x1)^2)) - y1*(x2 - x1)^3) + / (x2 - x1)^3) := by + have hp0 := onCurve_ne_zero h1 + have hq0 := onCurve_ne_zero h2 + have hd : x2 - x1 ≠ 0 := sub_ne_zero.mpr (Ne.symm hne) + unfold add + dsimp only + rw [if_neg hp0, if_neg hq0, if_neg hne] + rw [Prod.mk.injEq] + constructor <;> field_simp <;> ring + +/-- The affine sum in the doubling branch (`a = 0`), with cleared denominators. -/ +theorem add_dbl {x1 y1 : Fq} (h1 : OnCurve 0 5 (x1, y1)) (hy : y1 ≠ 0) : + add 0 (x1, y1) (x1, y1) + = ((9*x1^4 - 8*x1*y1^2) / (4*y1^2), + (3*x1^2*(12*x1*y1^2 - 9*x1^4) - 8*y1^4) / (8*y1^3)) := by + have hp0 := onCurve_ne_zero h1 + have h2y : y1 + y1 ≠ 0 := by rw [← two_mul]; exact mul_ne_zero (by decide) hy + have h2 : (2 : Fq) ≠ 0 := by decide + have h4 : (4 : Fq) ≠ 0 := by decide + have h8 : (8 : Fq) ≠ 0 := by decide + unfold add + dsimp only + rw [if_neg hp0, if_neg hp0, if_pos rfl, if_neg h2y] + rw [Prod.mk.injEq] + constructor <;> field_simp <;> ring + +/-! ## The core equivalence -/ + +/-- **Core equivalence on `Z = 1` representatives.** For on-curve `(x₁,y₁),(x₂,y₂)`, the RCB +projective sum interprets affinely as the complete affine sum, and stays a valid projective point. -/ +theorem padd_spec_z1 {x1 y1 x2 y2 : Fq} (h1 : OnCurve 0 5 (x1, y1)) (h2 : OnCurve 0 5 (x2, y2)) : + aff (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩) = add 0 (x1, y1) (x2, y2) + ∧ Valid (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩) := by + have e1 : y1^2 = x1^3 + 5 := by have h := h1; simp only [OnCurve] at h; linear_combination h + have e2 : y2^2 = x2^3 + 5 := by have h := h2; simp only [OnCurve] at h; linear_combination h + have hoc : OnCurveP (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩) := oncurveP_padd_z1 e1 e2 + have hy1 : y1 ≠ 0 := by + intro h; exact Vesta.no_onCurve_y_zero x1 (by rw [← h]; exact h1) + by_cases hx : x1 = x2 + · subst hx + by_cases hy : y1 + y2 = 0 + · -- inverse: y2 = -y1 + have hy2 : y2 = -y1 := by linear_combination hy + subst hy2 + refine ⟨?_, hoc, ?_⟩ + · have hz0 : (padd ⟨x1, y1, 1⟩ ⟨x1, -y1, 1⟩).Z = 0 := Z_zero_inv x1 y1 + rw [aff, if_pos hz0] + have hnp : (x1, -y1) = CompElliptic.CurveForms.ShortWeierstrass.neg (x1, y1) := by + simp [CompElliptic.CurveForms.ShortWeierstrass.neg] + rw [hnp, CompElliptic.CurveForms.ShortWeierstrass.add_neg] + · exact Or.inr (Or.inl (Y_ne_zero_inv e1 hy1)) + · -- doubling: y1 = y2 + have hyeq : y1 = y2 := by + have hz : (y1 - y2) * (y1 + y2) = 0 := by linear_combination e1 - e2 + rcases mul_eq_zero.mp hz with h | h + · exact sub_eq_zero.mp h + · exact absurd h hy + subst hyeq + have hZv : (padd ⟨x1, y1, 1⟩ ⟨x1, y1, 1⟩).Z = 8 * y1^3 := Z_doubling e1 + have hZne : (padd ⟨x1, y1, 1⟩ ⟨x1, y1, 1⟩).Z ≠ 0 := by + rw [hZv]; exact mul_ne_zero (by decide) (pow_ne_zero 3 hy1) + have hy2 : (4 : Fq) * y1^2 ≠ 0 := mul_ne_zero (by decide) (pow_ne_zero 2 hy1) + have hy3 : (8 : Fq) * y1^3 ≠ 0 := mul_ne_zero (by decide) (pow_ne_zero 3 hy1) + refine ⟨?_, hoc, ?_⟩ + · rw [aff, if_neg hZne, add_dbl h1 hy1, Prod.mk.injEq] + refine ⟨?_, ?_⟩ + · rw [div_eq_div_iff hZne hy2]; linear_combination keyx_dbl e1 + · rw [div_eq_div_iff hZne hy3]; linear_combination keyy_dbl e1 + · exact Or.inr (Or.inr hZne) + · -- distinct x + have hZne : (padd ⟨x1, y1, 1⟩ ⟨x2, y2, 1⟩).Z ≠ 0 := Z_ne_zero_dist e1 e2 hx + have hd : x2 - x1 ≠ 0 := sub_ne_zero.mpr (Ne.symm hx) + have hd2 : (x2 - x1)^2 ≠ 0 := pow_ne_zero 2 hd + have hd3 : (x2 - x1)^3 ≠ 0 := pow_ne_zero 3 hd + refine ⟨?_, hoc, ?_⟩ + · rw [aff, if_neg hZne, add_dist h1 h2 hx, Prod.mk.injEq] + refine ⟨?_, ?_⟩ + · rw [div_eq_div_iff hZne hd2]; linear_combination keyx_dist e1 e2 + · rw [div_eq_div_iff hZne hd3]; linear_combination keyy_dist e1 e2 + · exact Or.inr (Or.inr hZne) + +/-! ## Structural facts about `Valid` points and the identity-input reductions -/ + +/-- A valid point with `Z = 0` has `X = 0` (only the identity is at infinity). -/ +theorem X_zero_of_Z_zero {P : PVes} (h : OnCurveP P) (hz : P.Z = 0) : P.X = 0 := by + have hX3 : P.X ^ 3 = 0 := by rw [OnCurveP, hz] at h; simpa using h.symm + exact pow_eq_zero_iff (by norm_num : (3 : ℕ) ≠ 0) |>.mp hX3 + +/-- A valid point always has nonzero `Y` (no 2-torsion; the identity is `(0 : 1 : 0)`). -/ +theorem Valid.Y_ne_zero {P : PVes} (hP : Valid P) : P.Y ≠ 0 := by + rcases eq_or_ne P.Z 0 with hz | hz + · have hX := X_zero_of_Z_zero hP.1 hz + rcases hP.2 with h | h | h + · exact absurd hX h + · exact h + · exact absurd hz h + · intro hY + refine Vesta.no_onCurve_y_zero (P.X / P.Z) ?_ + have key : P.X^3 + 5 * P.Z^3 = 0 := by + have h := hP.1; rw [OnCurveP, hY] at h; linear_combination -h + show (0 : Fq)^2 = (P.X / P.Z)^3 + Vesta.a * (P.X / P.Z) + Vesta.b + rw [Vesta.a, Vesta.b]; field_simp; linear_combination -key + +/-- `padd` with the first argument an identity-type point `(0 : Y₁ : 0)` rescales the second. -/ +theorem padd_idL {P Q : PVes} (hX : P.X = 0) (hZ : P.Z = 0) : + padd P Q = smul (P.Y^2 * Q.Y) Q := by + simp only [padd, smul, hX, hZ, PVes.mk.injEq]; refine ⟨?_, ?_, ?_⟩ <;> ring + +/-- `padd` with the second argument an identity-type point `(0 : Y₂ : 0)` rescales the first. -/ +theorem padd_idR {P Q : PVes} (hX : Q.X = 0) (hZ : Q.Z = 0) : + padd P Q = smul (P.Y * Q.Y^2) P := by + simp only [padd, smul, hX, hZ, PVes.mk.injEq]; refine ⟨?_, ?_, ?_⟩ <;> ring + +/-- On-curve reading of the normalized affine coordinates of a finite valid point. -/ +theorem onCurve_norm {P : PVes} (h : OnCurveP P) (hz : P.Z ≠ 0) : + OnCurve 0 5 (P.X / P.Z, P.Y / P.Z) := by + have key : P.Y^2 * P.Z = P.X^3 + 5 * P.Z^3 := h + show (P.Y / P.Z)^2 = (P.X / P.Z)^3 + 0 * (P.X / P.Z) + 5 + field_simp; linear_combination key + +/-- **General equivalence.** For any two valid projective points, the RCB sum interprets +affinely as the complete affine sum, and stays valid. -/ +theorem padd_spec {P Q : PVes} (hP : Valid P) (hQ : Valid Q) : + aff (padd P Q) = add 0 (aff P) (aff Q) ∧ Valid (padd P Q) := by + have hPY := hP.Y_ne_zero + have hQY := hQ.Y_ne_zero + by_cases hPz : P.Z = 0 + · -- P is an identity-type point + have hPX := X_zero_of_Z_zero hP.1 hPz + have hlam : P.Y^2 * Q.Y ≠ 0 := mul_ne_zero (pow_ne_zero 2 hPY) hQY + have haffP : aff P = (0, 0) := by rw [aff, if_pos hPz] + rw [padd_idL hPX hPz, haffP, CompElliptic.CurveForms.ShortWeierstrass.zero_add] + refine ⟨aff_smul _ hlam Q, oncurveP_smul _ hQ.1, nezVec_smul _ hlam hQ.2⟩ + · by_cases hQz : Q.Z = 0 + · -- Q is an identity-type point + have hQX := X_zero_of_Z_zero hQ.1 hQz + have hlam : P.Y * Q.Y^2 ≠ 0 := mul_ne_zero hPY (pow_ne_zero 2 hQY) + have haffQ : aff Q = (0, 0) := by rw [aff, if_pos hQz] + rw [padd_idR hQX hQz, haffQ, CompElliptic.CurveForms.ShortWeierstrass.add_zero] + refine ⟨aff_smul _ hlam P, oncurveP_smul _ hP.1, nezVec_smul _ hlam hP.2⟩ + · -- both finite: normalize to `Z = 1` and use the core lemma + have hP1 := onCurve_norm hP.1 hPz + have hQ1 := onCurve_norm hQ.1 hQz + obtain ⟨haff, hvalid⟩ := padd_spec_z1 hP1 hQ1 + have hu : P.Z^2 * Q.Z^2 ≠ 0 := mul_ne_zero (pow_ne_zero 2 hPz) (pow_ne_zero 2 hQz) + have hpadd : padd P Q = smul (P.Z^2 * Q.Z^2) (padd ⟨P.X/P.Z, P.Y/P.Z, 1⟩ ⟨Q.X/Q.Z, Q.Y/Q.Z, 1⟩) := by + conv_lhs => rw [eq_smul_normalize hPz, eq_smul_normalize hQz] + rw [smul_padd] + have haffP : aff P = (P.X / P.Z, P.Y / P.Z) := by rw [aff, if_neg hPz] + have haffQ : aff Q = (Q.X / Q.Z, Q.Y / Q.Z) := by rw [aff, if_neg hQz] + refine ⟨?_, ?_⟩ + · rw [hpadd, aff_smul _ hu, haff, haffP, haffQ] + · rw [hpadd] + exact ⟨oncurveP_smul _ hvalid.1, nezVec_smul _ hu hvalid.2⟩ + +/-- Validity is preserved by `padd` (closure). -/ +theorem valid_padd {P Q : PVes} (hP : Valid P) (hQ : Valid Q) : Valid (padd P Q) := + (padd_spec hP hQ).2 + +/-- `aff` is a homomorphism from `padd` to the affine group law on valid inputs. -/ +theorem aff_padd {P Q : PVes} (hP : Valid P) (hQ : Valid Q) : + aff (padd P Q) = add 0 (aff P) (aff Q) := + (padd_spec hP hQ).1 + +/-- The identity `(0 : 1 : 0)` is valid. -/ +theorem valid_pid : Valid pid := by + refine ⟨?_, Or.inr (Or.inl one_ne_zero)⟩ + simp [OnCurveP, pid] + +/-- `aff pid = 𝒪`. -/ +@[simp] theorem aff_pid : aff pid = ((0 : Fq), (0 : Fq)) := by simp [aff, pid] + +/-! ## Bridge to the affine group `G = SWPoint Vesta.curve` -/ + +/-- Interpret a projective point as an element of the affine group: a finite on-curve point maps to +`(X/Z, Y/Z)`, everything else (points at infinity, off-curve junk) to the identity `𝒪`. -/ +def toAffine (P : PVes) : G := + if h : OnCurveP P ∧ P.Z ≠ 0 then + ⟨P.X / P.Z, P.Y / P.Z, Or.inl (onCurve_norm h.1 h.2)⟩ + else 0 + +/-- On valid points, `toAffine` agrees coordinatewise with `aff`. -/ +theorem toAffine_coords {P : PVes} (hP : Valid P) : + ((toAffine P).x, (toAffine P).y) = aff P := by + rcases eq_or_ne P.Z 0 with hz | hz + · have h0 : toAffine P = 0 := by rw [toAffine, dif_neg]; rintro ⟨_, h⟩; exact h hz + rw [h0, aff, if_pos hz]; rfl + · rw [toAffine, dif_pos ⟨hP.1, hz⟩, aff, if_neg hz] + +@[simp] theorem toAffine_pid : toAffine pid = 0 := by + rw [toAffine, dif_neg]; rintro ⟨_, h⟩; exact h rfl + +/-- **The homomorphism.** `toAffine` carries the projective RCB addition to the affine group law. -/ +theorem toAffine_padd {P Q : PVes} (hP : Valid P) (hQ : Valid Q) : + toAffine (padd P Q) = toAffine P + toAffine Q := by + apply SWPoint.ext_pair + have hL : ((toAffine (padd P Q)).x, (toAffine (padd P Q)).y) = aff (padd P Q) := + toAffine_coords (valid_padd hP hQ) + have hRc : ((toAffine P + toAffine Q).x, (toAffine P + toAffine Q).y) + = add Vesta.curve.A ((toAffine P).x, (toAffine P).y) ((toAffine Q).x, (toAffine Q).y) := rfl + rw [hL, hRc, toAffine_coords hP, toAffine_coords hQ, aff_padd hP hQ] + rfl + +/-- Materialize an affine group element as a projective representative: the identity as `(0 : 1 : 0)`, +a finite point as `(x : y : 1)`. -/ +def ofAffine (p : G) : PVes := + if p.x = 0 ∧ p.y = 0 then pid else ⟨p.x, p.y, 1⟩ + +theorem valid_ofAffine (p : G) : Valid (ofAffine p) := by + rw [ofAffine] + by_cases h : p.x = 0 ∧ p.y = 0 + · rw [if_pos h]; exact valid_pid + · rw [if_neg h] + have hoc : OnCurve 0 5 (p.x, p.y) := by + have hv := p.onCurve + exact hv.resolve_right (by rw [Prod.mk.injEq]; exact h) + have e : p.y^2 = p.x^3 + 5 := by have := hoc; simp only [OnCurve] at this; linear_combination this + refine ⟨?_, Or.inr (Or.inr one_ne_zero)⟩ + simp only [OnCurveP]; linear_combination e + +theorem toAffine_ofAffine (p : G) : toAffine (ofAffine p) = p := by + rw [ofAffine] + by_cases h : p.x = 0 ∧ p.y = 0 + · rw [if_pos h, toAffine_pid] + apply SWPoint.ext_pair + obtain ⟨hx, hy⟩ := h + show ((0 : G).x, (0 : G).y) = (p.x, p.y) + rw [hx, hy]; rfl + · rw [if_neg h] + have hoc : OnCurve 0 5 (p.x, p.y) := p.onCurve.resolve_right (by rw [Prod.mk.injEq]; exact h) + have hZ : ((⟨p.x, p.y, 1⟩ : PVes)).Z ≠ 0 := one_ne_zero + have hocp : OnCurveP (⟨p.x, p.y, 1⟩ : PVes) := by + have e : p.y^2 = p.x^3 + 5 := by have := hoc; simp only [OnCurve] at this; linear_combination this + simp only [OnCurveP]; linear_combination e + apply SWPoint.ext_pair + rw [toAffine, dif_pos ⟨hocp, hZ⟩] + show (p.x / 1, p.y / 1) = (p.x, p.y) + rw [div_one, div_one] + +/-! ## Fast (binary) scalar multiplication -/ + +/-- Binary double-and-add scalar multiplication over `PVes` using the RCB `padd`, matching the +`CompElliptic.binNsmul` recurrence used by the affine group's `nsmul`. -/ +def pnsmulFast (n : ℕ) (P : PVes) : PVes := binNsmul padd pid n P + +/-- **Fast scalar multiplication is the genuine group action.** By strong induction along the +`binNsmul` double-and-add recurrence: each intermediate stays valid (`valid_padd`) and `toAffine` +is a homomorphism (`toAffine_padd`), so the projective ladder computes `n • (toAffine P)`. -/ +theorem pnsmulFast_spec {P : PVes} (hP : Valid P) : ∀ n : ℕ, + Valid (pnsmulFast n P) ∧ toAffine (pnsmulFast n P) = n • toAffine P := by + intro n + induction n using Nat.strong_induction_on with + | _ n ih => + rw [pnsmulFast, binNsmul] + split + · rename_i h; subst h + exact ⟨valid_pid, by rw [toAffine_pid, zero_nsmul]⟩ + · rename_i hn + have hlt : n / 2 < n := Nat.div_lt_self (Nat.pos_of_ne_zero hn) (by decide) + obtain ⟨hvq, hq⟩ := ih (n / 2) hlt + rw [pnsmulFast] at hvq hq + set q := binNsmul padd pid (n / 2) P with hqdef + have hvd : Valid (padd q q) := valid_padd hvq hvq + have hdaff : toAffine (padd q q) = (n / 2) • toAffine P + (n / 2) • toAffine P := by + rw [toAffine_padd hvq hvq, hq] + have key : n • toAffine P + = (n / 2) • toAffine P + (n / 2) • toAffine P + (n % 2) • toAffine P := by + rw [← add_nsmul, ← add_nsmul]; congr 1; omega + split + · rename_i hodd + refine ⟨valid_padd hvd hP, ?_⟩ + rw [toAffine_padd hvd hP, hdaff, key, hodd, one_nsmul] + · rename_i heven + refine ⟨hvd, ?_⟩ + rw [hdaff, key, show n % 2 = 0 from by omega, zero_nsmul, _root_.add_zero] + +/-- The packaged fast scalar multiplication on `G`, going through projective coordinates. -/ +def smulFast (n : ℕ) (p : G) : G := toAffine (pnsmulFast n (ofAffine p)) + +/-- **Correctness of `smulFast`:** it equals the affine scalar multiplication `n • p`, with the +field inversion paid once (in the final `toAffine`) instead of once per addition. -/ +theorem smulFast_eq (n : ℕ) (p : G) : smulFast n p = n • p := by + rw [smulFast, (pnsmulFast_spec (valid_ofAffine p) n).2, toAffine_ofAffine] + +/-! ## Fast compiled spelling of `padd` (compiled-tier only) + +Compiled generically, every field operation in `padd` goes through boxed `CommRing (ZMod q)` +dictionary projections. `paddFast` is the same closed forms over raw `ℕ` representatives with +fused mul-mod, and `padd_eq_paddFast` is the proven equality — the proven counterpart of +`implemented_by` — so `padd` stays the statement surface while compiled call sites run fast. -/ + +/-- The Vesta base-field order, the modulus of the raw-`ℕ` fast path (reducible so `ZMod qv` +unifies with `Fq`). -/ +private abbrev qv : ℕ := CompElliptic.Fields.Pasta.PALLAS_SCALAR_CARD + +private theorem qv_pos : 0 < qv := by decide + +private instance : NeZero qv := ⟨qv_pos.ne'⟩ + +/-- Fused modular addition on canonical representatives. -/ +@[inline] private def fadd (a b : ℕ) : ℕ := (a + b) % qv + +/-- Fused modular multiplication: one GMP multiply, one GMP mod. -/ +@[inline] private def fmul (a b : ℕ) : ℕ := (a * b) % qv + +/-- Modular subtraction without hypotheses: `a + (q − b mod q) ≡ a − b`. -/ +@[inline] private def fsub (a b : ℕ) : ℕ := (a + (qv - b % qv)) % qv + +private theorem cast_fadd (a b : ℕ) : ((fadd a b : ℕ) : Fq) = (a : Fq) + (b : Fq) := by + rw [fadd, ZMod.natCast_mod, Nat.cast_add] + +private theorem cast_fmul (a b : ℕ) : ((fmul a b : ℕ) : Fq) = (a : Fq) * (b : Fq) := by + rw [fmul, ZMod.natCast_mod, Nat.cast_mul] + +private theorem cast_fsub (a b : ℕ) : ((fsub a b : ℕ) : Fq) = (a : Fq) - (b : Fq) := by + rw [fsub, ZMod.natCast_mod, Nat.cast_add, Nat.cast_sub (Nat.mod_lt _ qv_pos).le, + ZMod.natCast_mod, ZMod.natCast_self] + ring + +private theorem cast_val (a : Fq) : ((ZMod.val a : ℕ) : Fq) = a := + ZMod.natCast_rightInverse a + +/-- `padd` over raw `ℕ` representatives: same RCB closed forms, fused mul-mod, shared +subproducts. Statement-surface code should keep calling `padd`; the compiler substitutes this +body via `padd_eq_paddFast`. -/ +def paddFast (P Q : PVes) : PVes := + let x1 := P.X.val; let y1 := P.Y.val; let z1 := P.Z.val + let x2 := Q.X.val; let y2 := Q.Y.val; let z2 := Q.Z.val + let y2sq := fmul y2 y2 + let z2sq := fmul z2 z2 + let x1y1 := fmul x1 y1 + let y1sq := fmul y1 y1 + let z1sq := fmul z1 z1 + let x2y2 := fmul x2 y2 + let y1z1 := fmul y1 z1 + let x1z1 := fmul x1 z1 + let y2z2 := fmul y2 z2 + let x2z2 := fmul x2 z2 + let x1sq := fmul x1 x1 + let x2sq := fmul x2 x2 + let X3 := fsub (fadd (fsub (fsub (fmul x1y1 y2sq) (fmul 15 (fmul x1y1 z2sq))) + (fmul 30 (fmul x1z1 y2z2))) + (fsub (fmul y1sq x2y2) (fmul 15 (fmul z1sq x2y2)))) + (fmul 30 (fmul y1z1 x2z2)) + let Y3 := fsub (fadd (fadd (fmul y1sq y2sq) (fmul 45 (fmul x1sq x2z2))) + (fmul 45 (fmul x1z1 x2sq))) + (fmul 225 (fmul z1sq z2sq)) + let Z3 := fadd (fadd (fadd (fadd (fadd (fmul y1sq y2z2) (fmul y1z1 y2sq)) + (fmul 3 (fmul x1sq x2y2))) + (fmul 3 (fmul x1y1 x2sq))) + (fmul 15 (fmul y1z1 z2sq))) + (fmul 15 (fmul z1sq y2z2)) + ⟨(X3 : Fq), (Y3 : Fq), (Z3 : Fq)⟩ + +/-- **The fast spelling is `padd`.** NOT registered `@[csimp]`: in the current +`native_decide` tier, library calls execute through interpreted IR, where this raw-`ℕ` +spelling measured ~17× SLOWER than the typeclass path it replaces (1.5 ms vs 85 µs per +add — many small interpreted steps lose to few boxed GMP calls; the same effect that +makes interpreted Montgomery a 40× regression). The equality is kept for the compiled +tier: once the fast-field lib ships precompiled (`precompileModules` dylib), registering +this `@[csimp]` (or the Montgomery successor) makes every compiled call site (including +`native_decide` auxiliaries) run the fused raw-`ℕ` form. -/ +theorem padd_eq_paddFast : @padd = @paddFast := by + funext P Q + simp only [padd, paddFast, PVes.mk.injEq] + refine ⟨?_, ?_, ?_⟩ <;> + simp only [cast_fadd, cast_fmul, cast_fsub, cast_val, Nat.cast_ofNat] <;> ring + +end PVes +end CompElliptic.Curves.Pasta.Fast.Projective diff --git a/CompElliptic/Curves/Pasta/Fast/ProjectiveMontDefs.lean b/CompElliptic/Curves/Pasta/Fast/ProjectiveMontDefs.lean new file mode 100644 index 0000000..b14d89a --- /dev/null +++ b/CompElliptic/Curves/Pasta/Fast/ProjectiveMontDefs.lean @@ -0,0 +1,121 @@ +/- +Copyright (c) 2026 CompElliptic Contributors. All rights reserved. +Released under the Apache License, Version 2.0, or the MIT license, at your option, +as described in the files LICENSE-APACHE and LICENSE-MIT. +Authors: Gregor Mitscha-Baude +-/ +import CompElliptic.Vendor.CompPoly.Montgomery.Native64x8Defs + +/-! +# Vesta projective point arithmetic over Montgomery limbs: runtime definitions (core-only) + +The projective Vesta arithmetic of `CompElliptic.Curves.Pasta.Fast.Projective` and the Pippenger +schedules of `CompElliptic.Curves.Pasta.Fast.MsmProj`, transcribed operation for operation onto +eight-limb Montgomery residues so that they can be native-compiled (the `FastFieldNative` lane in +the lakefile — hence no imports beyond Lean core, and no proofs here). Every operation is proven +to compute its `𝔽_q`-valued counterpart in `CompElliptic.Curves.Pasta.Fast.ProjectiveMontEquiv`. + +The `Fast` interfaces are provisional: they are not guaranteed to remain public, and +may be folded into the existing API or otherwise changed incompatibly. +-/ + +namespace CompElliptic.Curves.Pasta.Fast.ProjectiveMont + +open Montgomery.Native64x8 (Limbs8) +open Montgomery.Native64x8.VestaFq + +/-- A projective point in `(X : Y : Z)` coordinates, each coordinate an eight-limb Montgomery +residue of the Vesta base field. -/ +structure PM where + /-- The `X` coordinate. -/ + X : Limbs8 + /-- The `Y` coordinate. -/ + Y : Limbs8 + /-- The `Z` coordinate. -/ + Z : Limbs8 +deriving DecidableEq + +namespace PM + +/-- `Array.get!`/`set!` need a default; the projective identity is the right one. -/ +instance : Inhabited PM := ⟨⟨zero, one, zero⟩⟩ + +/-- The Montgomery residue of `3`. -/ +def c3 : Limbs8 := ofNat 3 +/-- The Montgomery residue of `15` (`b3 = 3b` for Vesta). -/ +def c15 : Limbs8 := ofNat 15 +/-- The Montgomery residue of `30`. -/ +def c30 : Limbs8 := ofNat 30 +/-- The Montgomery residue of `45`. -/ +def c45 : Limbs8 := ofNat 45 +/-- The Montgomery residue of `225`. -/ +def c225 : Limbs8 := ofNat 225 + +/-- Renes–Costello–Batina complete addition (`add-2015-rcb`, `a = 0`, `b3 = 15`), on +Montgomery limbs. -/ +@[inline] def padd (P Q : PM) : PM where + X := + sub (sub (add (sub (sub (mul (mul (P.X) (P.Y)) (square (Q.Y))) (mul (mul (mul (c15) (P.X)) + (P.Y)) (square (Q.Z)))) (mul (mul (mul (mul (c30) (P.X)) (P.Z)) (Q.Y)) (Q.Z))) (mul (mul + (square (P.Y)) (Q.X)) (Q.Y))) (mul (mul (mul (c15) (square (P.Z))) (Q.X)) (Q.Y))) (mul + (mul (mul (mul (c30) (P.Y)) (P.Z)) (Q.X)) (Q.Z)) + Y := + sub (add (add (mul (square (P.Y)) (square (Q.Y))) (mul (mul (mul (c45) (square (P.X))) + (Q.X)) (Q.Z))) (mul (mul (mul (c45) (P.X)) (P.Z)) (square (Q.X)))) (mul (mul (c225) + (square (P.Z))) (square (Q.Z))) + Z := + add (add (add (add (add (mul (mul (square (P.Y)) (Q.Y)) (Q.Z)) (mul (mul (P.Y) (P.Z)) + (square (Q.Y)))) (mul (mul (mul (c3) (square (P.X))) (Q.X)) (Q.Y))) (mul (mul (mul (c3) + (P.X)) (P.Y)) (square (Q.X)))) (mul (mul (mul (c15) (P.Y)) (P.Z)) (square (Q.Z)))) (mul + (mul (mul (c15) (square (P.Z))) (Q.Y)) (Q.Z)) + +/-- The projective identity `𝒪 = (0 : 1 : 0)`. -/ +def pid : PM := ⟨zero, one, zero⟩ + +/-! ## Group kernels + +Mirroring `MsmProj`'s schedules step for step: only the interpretation of a coordinate changes, +never the schedule, so the correctness proofs stay structural inductions. -/ + +/-- Fixed 256-step LSB-first double-and-add ladder. -/ +def pnsmul (n : Nat) (p : PM) : PM := + (List.range 256).foldl + (fun (st : PM × PM) i => + let acc := if (n >>> i) &&& 1 = 1 then padd st.1 st.2 else st.1 + (acc, padd st.2 st.2)) + (pid, p) |>.1 + +/-- One Array-scatter step (mirrors `MsmProj.pscatterStep`). -/ +def scatterStep (a : Array PM) (p : Nat × PM) : Array PM := + if p.1 = 0 then a else a.modify (p.1 - 1) (fun v => padd v p.2) + +/-- Scatter a digit-tagged point list into its `base − 1` buckets in ONE pass. -/ +def bucketScatter (base : Nat) (dp : List (Nat × PM)) : Array PM := + dp.foldl scatterStep (Array.replicate (base - 1) pid) + +/-- One step of the bucket downsweep (mirrors `MsmProj.paccStep`). -/ +def accStep (a : PM) (p : PM × PM) : PM × PM := (padd p.1 a, padd p.2 (padd p.1 a)) + +/-- The window-`i` value in base `base` (mirrors `MsmProj.pwindowValueFast`). -/ +def windowValue (base i : Nat) (terms : List (Nat × PM)) : PM := + let scale := base ^ i + (List.foldr accStep (pid, pid) + (bucketScatter base (terms.map fun t => (t.1 / scale % base, t.2))).toList).2 + +/-- `c`-fold doubling — the `base •` Horner step between adjacent windows. -/ +def pdoublings (c : Nat) (p : PM) : PM := (List.range c).foldl (fun a _ => padd a a) p + +/-- Windowed Pippenger MSM, window `c` (mirrors `MsmProj.pippengerProjScatter`, with the window +count fixed at `⌈256 / c⌉` and the `base •` step spelled as `c` doublings). -/ +def msm (c : Nat) (terms : List (Nat × PM)) : PM := + let numWindows := (256 + c - 1) / c + let base := 2 ^ c + ((List.range numWindows).map fun i => windowValue base i terms).foldr + (fun v acc => padd (pdoublings c acc) v) pid + +/-- Projective negation: `-(X : Y : Z) = (X : −Y : Z)`. -/ +@[inline] def pneg (p : PM) : PM := ⟨p.X, sub zero p.Y, p.Z⟩ + +end PM + +end CompElliptic.Curves.Pasta.Fast.ProjectiveMont diff --git a/CompElliptic/Curves/Pasta/Fast/ProjectiveMontEquiv.lean b/CompElliptic/Curves/Pasta/Fast/ProjectiveMontEquiv.lean new file mode 100644 index 0000000..6586478 --- /dev/null +++ b/CompElliptic/Curves/Pasta/Fast/ProjectiveMontEquiv.lean @@ -0,0 +1,636 @@ +/- +Copyright (c) 2026 CompElliptic Contributors. All rights reserved. +Released under the Apache License, Version 2.0, or the MIT license, at your option, +as described in the files LICENSE-APACHE and LICENSE-MIT. +Authors: Gregor Mitscha-Baude +-/ +import CompElliptic.Curves.Pasta.Fast.ProjectiveMontDefs +import CompElliptic.Vendor.CompPoly.Montgomery.Pasta +import CompElliptic.Curves.Pasta.Fast.MsmProj + +/-! +# The Montgomery kernel is the proven projective arithmetic + +`CompElliptic.Curves.Pasta.Fast.ProjectiveMontDefs` transcribes the projective Vesta arithmetic +onto eight-limb Montgomery residues; this module proves that the transcription computes the +statement-surface functions of `Projective.lean` / `MsmProj.lean`. There are two tiers, not +three: the bridge is the coordinatewise `montVal`, which is the vendored field's ring isomorphism +`FastField.toField` read off raw limbs, correct on well-formed residues (`WF`). + +`toPVesM_padd` is the one commuting square of substance — the ~40 field operations of the +Renes–Costello–Batina formulas pushed through the ring homomorphism. Above it everything is +structural: the schedules mirror `MsmProj`'s, so they transport along +`RM p P := WFP p ∧ toPVesM p = P` (and, where the group law is needed, along `toGM`). + +## Main results + +* `toPVesM_padd`, `toPVesM_pid`, `toPVesM_pneg` — addition, identity and negation +* `pnsmulM_spec` — the 256-step ladder is `n • ·` in the affine group, for `n < 2 ^ 256` +* `msmM_spec` — the windowed Pippenger MSM is `Fast.Msm.pippenger` + +`RM`, `RM2`, `RA` and the fold combinators are public rather than private: lifting a *further* +shared schedule to the Montgomery tier needs the same inductions. + +The `Fast` interfaces are provisional: they are not guaranteed to remain public, and +may be folded into the existing API or otherwise changed incompatibly. +-/ + +namespace CompElliptic.Curves.Pasta.Fast.ProjectiveMont + +open Montgomery.Native64x8 +open CompElliptic.Curves.Pasta.Fast +open CompElliptic.Curves.Pasta.Fast.Projective +open CompElliptic.Curves.Pasta.Fast.Projective.PVes +open CompElliptic.Fields.Pasta (PALLAS_SCALAR_CARD) + +/-- The Vesta base field, the ambient field of the projective statement surface. -/ +local notation "Fq" => CompElliptic.Curves.Pasta.Fast.Projective.Fq + +/-! ## The field level + +Every `VestaFq` entry point is the corresponding `FastField` operation on the nose, so the whole +field tier is `rfl`-transport into the vendored proofs. -/ + +/-- A well-formed Montgomery residue: bounded limbs holding a value below `q`. This is exactly +the carrier property of `Montgomery.Native64x8.FastField PALLAS_SCALAR_CARD`. -/ +def WF (x : Limbs8) : Prop := x.Bounded ∧ x.toNat < PALLAS_SCALAR_CARD + +/-- The field value of a Montgomery residue: the residue divided by the radix `R = 2 ^ 256`. -/ +def montVal (x : Limbs8) : Fq := (x.toNat : Fq) * ((2 ^ 256 : ℕ) : Fq)⁻¹ + +/-- A well-formed residue, packaged as an element of the proven fast field. -/ +private def toFF (x : Limbs8) (h : WF x) : FastField PALLAS_SCALAR_CARD := ⟨x, h⟩ + +theorem montVal_eq_toField {x : Limbs8} (h : WF x) : + montVal x = FastField.toField (toFF x h) := by + rw [montVal] + exact (FastField.toField_eq (toFF x h)).symm + +/-! ### Per-operation cast lemmas -/ + +theorem wf_add {a b : Limbs8} (ha : WF a) (hb : WF b) : WF (VestaFq.add a b) := + (FastField.add (toFF a ha) (toFF b hb)).property + +theorem montVal_add {a b : Limbs8} (ha : WF a) (hb : WF b) : + montVal (VestaFq.add a b) = montVal a + montVal b := by + rw [montVal_eq_toField (wf_add ha hb), montVal_eq_toField ha, montVal_eq_toField hb] + exact FastField.toField_add (toFF a ha) (toFF b hb) + +theorem wf_sub {a b : Limbs8} (ha : WF a) (hb : WF b) : WF (VestaFq.sub a b) := + (FastField.sub (toFF a ha) (toFF b hb)).property + +theorem montVal_sub {a b : Limbs8} (ha : WF a) (hb : WF b) : + montVal (VestaFq.sub a b) = montVal a - montVal b := by + rw [montVal_eq_toField (wf_sub ha hb), montVal_eq_toField ha, montVal_eq_toField hb] + exact FastField.toField_sub (toFF a ha) (toFF b hb) + +theorem wf_neg {a : Limbs8} (ha : WF a) : WF (VestaFq.neg a) := + (FastField.neg (toFF a ha)).property + +theorem montVal_neg {a : Limbs8} (ha : WF a) : montVal (VestaFq.neg a) = -montVal a := by + rw [montVal_eq_toField (wf_neg ha), montVal_eq_toField ha] + exact FastField.toField_neg (toFF a ha) + +theorem wf_mul {a b : Limbs8} (ha : WF a) (hb : WF b) : WF (VestaFq.mul a b) := + (FastField.mul (toFF a ha) (toFF b hb)).property + +theorem montVal_mul {a b : Limbs8} (ha : WF a) (hb : WF b) : + montVal (VestaFq.mul a b) = montVal a * montVal b := by + rw [montVal_eq_toField (wf_mul ha hb), montVal_eq_toField ha, montVal_eq_toField hb] + exact FastField.toField_mul (toFF a ha) (toFF b hb) + +theorem wf_square {a : Limbs8} (ha : WF a) : WF (VestaFq.square a) := wf_mul ha ha + +theorem montVal_square {a : Limbs8} (ha : WF a) : + montVal (VestaFq.square a) = montVal a * montVal a := montVal_mul ha ha + +theorem wf_zero : WF VestaFq.zero := (FastField.zero PALLAS_SCALAR_CARD).property + +theorem montVal_zero : montVal VestaFq.zero = 0 := by + rw [montVal_eq_toField wf_zero] + exact FastField.toField_zero + +theorem wf_one : WF VestaFq.one := (FastField.one PALLAS_SCALAR_CARD).property + +theorem montVal_one : montVal VestaFq.one = 1 := by + rw [montVal_eq_toField wf_one] + exact FastField.toField_one + +/-- **Entering Montgomery form is the canonical cast**, for canonical naturals. This is +`FastField.ofCanonicalNat`'s correctness read off raw limbs. -/ +theorem wf_ofNat {k : ℕ} (h : k < PALLAS_SCALAR_CARD) : WF (VestaFq.ofNat k) := + (FastField.ofCanonicalNat k h).property + +theorem montVal_ofNat {k : ℕ} (h : k < PALLAS_SCALAR_CARD) : + montVal (VestaFq.ofNat k) = (k : Fq) := by + rw [montVal_eq_toField (wf_ofNat h)] + exact FastField.toField_ofCanonicalNat h + +/-! ### The RCB coefficients -/ + +theorem wf_c3 : WF PM.c3 := wf_ofNat (by decide) +theorem wf_c15 : WF PM.c15 := wf_ofNat (by decide) +theorem wf_c30 : WF PM.c30 := wf_ofNat (by decide) +theorem wf_c45 : WF PM.c45 := wf_ofNat (by decide) +theorem wf_c225 : WF PM.c225 := wf_ofNat (by decide) + +theorem montVal_c3 : montVal PM.c3 = 3 := by + rw [show PM.c3 = VestaFq.ofNat 3 from rfl, montVal_ofNat (by decide)]; norm_num + +theorem montVal_c15 : montVal PM.c15 = 15 := by + rw [show PM.c15 = VestaFq.ofNat 15 from rfl, montVal_ofNat (by decide)]; norm_num + +theorem montVal_c30 : montVal PM.c30 = 30 := by + rw [show PM.c30 = VestaFq.ofNat 30 from rfl, montVal_ofNat (by decide)]; norm_num + +theorem montVal_c45 : montVal PM.c45 = 45 := by + rw [show PM.c45 = VestaFq.ofNat 45 from rfl, montVal_ofNat (by decide)]; norm_num + +theorem montVal_c225 : montVal PM.c225 = 225 := by + rw [show PM.c225 = VestaFq.ofNat 225 from rfl, montVal_ofNat (by decide)]; norm_num + +/-! ## The point level -/ + +/-- Coordinatewise interpretation of a Montgomery triple as a projective point over `𝔽_q`. -/ +def toPVesM (p : PM) : PVes := ⟨montVal p.X, montVal p.Y, montVal p.Z⟩ + +/-- All three coordinates are well-formed Montgomery residues. -/ +def WFP (p : PM) : Prop := WF p.X ∧ WF p.Y ∧ WF p.Z + +/-- The affine reading of a Montgomery triple. -/ +def toGM (p : PM) : G := toAffine (toPVesM p) + +@[simp] theorem toPVesM_X (p : PM) : (toPVesM p).X = montVal p.X := rfl +@[simp] theorem toPVesM_Y (p : PM) : (toPVesM p).Y = montVal p.Y := rfl +@[simp] theorem toPVesM_Z (p : PM) : (toPVesM p).Z = montVal p.Z := rfl + +theorem wfp_pid : WFP PM.pid := ⟨wf_zero, wf_one, wf_zero⟩ + +theorem wfp_pneg {p : PM} (h : WFP p) : WFP (PM.pneg p) := + ⟨h.1, wf_sub wf_zero h.2.1, h.2.2⟩ + +set_option maxRecDepth 8000 in +/-- **Well-formedness is closed under the kernel's RCB addition.** Every node of the +coordinate expressions is one of the four field entry points, each of which ships its own +carrier proof. The closure search runs `with_reducible`, so a candidate rule is rejected on +the head symbol instead of unfolding the CIOS rounds underneath it. -/ +theorem wfp_padd {p r : PM} (hp : WFP p) (hr : WFP r) : WFP (PM.padd p r) := by + obtain ⟨hpx, hpy, hpz⟩ := hp + obtain ⟨hrx, hry, hrz⟩ := hr + refine ⟨?_, ?_, ?_⟩ <;> + · simp only [PM.padd] + repeat' first + | (with_reducible assumption) + | (with_reducible exact wf_c3) + | (with_reducible exact wf_c15) + | (with_reducible exact wf_c30) + | (with_reducible exact wf_c45) + | (with_reducible exact wf_c225) + | (with_reducible apply wf_square) + | (with_reducible apply wf_mul) + | (with_reducible apply wf_add) + | (with_reducible apply wf_sub) + +set_option maxRecDepth 8000 in +/-- **The kernel's addition is the projective addition.** Same Renes–Costello–Batina closed +forms, evaluated on Montgomery residues with the vendored field's operations. -/ +theorem toPVesM_padd {p r : PM} (hp : WFP p) (hr : WFP r) : + toPVesM (PM.padd p r) = PVes.padd (toPVesM p) (toPVesM r) := by + obtain ⟨hpx, hpy, hpz⟩ := hp + obtain ⟨hrx, hry, hrz⟩ := hr + simp only [toPVesM, PM.padd, PVes.padd, PVes.mk.injEq] + refine ⟨?_, ?_, ?_⟩ <;> + simp (maxDischargeDepth := 16) only [montVal_add, montVal_sub, montVal_mul, montVal_square, + montVal_c3, montVal_c15, montVal_c30, montVal_c45, montVal_c225, + wf_add, wf_sub, wf_mul, wf_square, wf_c3, wf_c15, wf_c30, wf_c45, wf_c225, + hpx, hpy, hpz, hrx, hry, hrz] <;> + ring + +@[simp] theorem toPVesM_pid : toPVesM PM.pid = PVes.pid := by + simp only [toPVesM, PM.pid, PVes.pid, PVes.mk.injEq] + exact ⟨montVal_zero, montVal_one, montVal_zero⟩ + +/-- **The kernel's negation is coordinate negation of `Y`**, the short-Weierstrass inverse. -/ +theorem toPVesM_pneg {p : PM} (h : WFP p) : + toPVesM (PM.pneg p) = ⟨(toPVesM p).X, -(toPVesM p).Y, (toPVesM p).Z⟩ := by + have hy : montVal (VestaFq.sub VestaFq.zero p.Y) = -montVal p.Y := by + rw [montVal_sub wf_zero h.2.1, montVal_zero, zero_sub] + simp only [toPVesM, PM.pneg, hy] + +/-! ### The affine reading + +`WV` bundles the two invariants every group schedule carries; under it the kernel's `padd` is the +affine group addition, which is all the ladder and the Horner recombination need. -/ + +/-- A well-formed Montgomery triple denoting a representable projective point. -/ +def WV (p : PM) : Prop := WFP p ∧ Valid (toPVesM p) + +theorem WV_pid : WV PM.pid := ⟨wfp_pid, by rw [toPVesM_pid]; exact valid_pid⟩ + +theorem WV_padd {p r : PM} (hp : WV p) (hr : WV r) : WV (PM.padd p r) := + ⟨wfp_padd hp.1 hr.1, by rw [toPVesM_padd hp.1 hr.1]; exact valid_padd hp.2 hr.2⟩ + +@[simp] theorem toGM_pid : toGM PM.pid = 0 := by rw [toGM, toPVesM_pid, toAffine_pid] + +/-- **The kernel's addition is the affine group addition** on representable points. -/ +theorem toGM_padd {p r : PM} (hp : WV p) (hr : WV r) : + toGM (PM.padd p r) = toGM p + toGM r := by + simp only [toGM] + rw [toPVesM_padd hp.1 hr.1, toAffine_padd hp.2 hr.2] + +/-! ## The correspondence with the projective statement surface + +`RM p P` says a Montgomery triple is well formed and denotes the projective point `P`; every +shared schedule preserves it. -/ + +/-- The kernel correspondence: a well-formed Montgomery triple denoting a given projective +point. -/ +def RM (p : PM) (P : PVes) : Prop := WFP p ∧ toPVesM p = P + +/-- The correspondence on ladder/downsweep state pairs. -/ +def RM2 (s : PM × PM) (s' : PVes × PVes) : Prop := RM s.1 s'.1 ∧ RM s.2 s'.2 + +theorem RM_self {p : PM} (h : WFP p) : RM p (toPVesM p) := ⟨h, rfl⟩ + +theorem RM_toG {p : PM} {P : PVes} (h : RM p P) : toGM p = toAffine P := by rw [toGM, h.2] + +/-- A corresponding pair inherits representability from the projective side. -/ +theorem RM_valid {p : PM} {P : PVes} (h : RM p P) (hv : Valid P) : WV p := + ⟨h.1, by rw [h.2]; exact hv⟩ + +theorem RM_padd {p r : PM} {P Q : PVes} (hp : RM p P) (hr : RM r Q) : + RM (PM.padd p r) (PVes.padd P Q) := by + refine ⟨wfp_padd hp.1 hr.1, ?_⟩ + rw [toPVesM_padd hp.1 hr.1, hp.2, hr.2] + +theorem RM_pid : RM PM.pid PVes.pid := ⟨wfp_pid, toPVesM_pid⟩ + +theorem RM_pneg {p : PM} {P : PVes} (h : RM p P) : + RM (PM.pneg p) ⟨P.X, -P.Y, P.Z⟩ := by + refine ⟨wfp_pneg h.1, ?_⟩ + rw [toPVesM_pneg h.1, h.2] + +/-! ### Relation-preserving folds -/ + +/-- A relation-preserving `foldl` whose step may use a property of the list elements. -/ +theorem foldl_rel₂ {σ σ' β : Type} (R : σ → σ' → Prop) (P : β → Prop) : + ∀ (l : List β) (f : σ → β → σ) (g : σ' → β → σ'), + (∀ s s' x, P x → R s s' → R (f s x) (g s' x)) → (∀ x ∈ l, P x) → + ∀ s s', R s s' → R (l.foldl f s) (l.foldl g s') := by + intro l + induction l with + | nil => intro _ _ _ _ s s' h; exact h + | cons x l ih => + intro f g hstep hP s s' h + exact ih f g hstep (fun y hy => hP y (by simp [hy])) _ _ + (hstep s s' x (hP x (by simp)) h) + +/-- A relation-preserving `foldr` whose step may use a property of the list elements. -/ +theorem foldr_rel₂ {σ σ' β : Type} (R : σ → σ' → Prop) (P : β → Prop) : + ∀ (l : List β) (f : β → σ → σ) (g : β → σ' → σ'), + (∀ x s s', P x → R s s' → R (f x s) (g x s')) → (∀ x ∈ l, P x) → + ∀ s s', R s s' → R (l.foldr f s) (l.foldr g s') := by + intro l + induction l with + | nil => intro _ _ _ _ s s' h; exact h + | cons x l ih => + intro f g hstep hP s s' h + exact hstep x _ _ (hP x (by simp)) + (ih f g hstep (fun y hy => hP y (by simp [hy])) s s' h) + +/-! ## The scalar ladder + +`PM.pnsmul` is LSB-first, unlike the statement-surface `pnsmulFast`, so it is proved where the +group law lives: through `toGM`, carrying the standard accumulator invariant. -/ + +/-- The ladder state after `i` steps: the accumulator holds `(n mod 2 ^ i) • A` and the base holds +`2 ^ i • A`, both well formed and representable. -/ +private def LadderInvM (A : G) (n i : ℕ) (st : PM × PM) : Prop := + WV st.1 ∧ WV st.2 ∧ toGM st.1 = (n % 2 ^ i) • A ∧ toGM st.2 = (2 ^ i : ℕ) • A + +private theorem ladderM_step {A : G} {n i : ℕ} {st : PM × PM} (h : LadderInvM A n i st) : + LadderInvM A n (i + 1) + ((if (n >>> i) &&& 1 = 1 then PM.padd st.1 st.2 else st.1), PM.padd st.2 st.2) := by + obtain ⟨hv1, hv2, ha1, ha2⟩ := h + have hbit : (n >>> i) &&& 1 = n / 2 ^ i % 2 := by + rw [Nat.shiftRight_eq_div_pow, Nat.and_one_is_mod] + have hmod : n % 2 ^ (i + 1) = n % 2 ^ i + 2 ^ i * (n / 2 ^ i % 2) := by + rw [pow_succ, Nat.mod_mul] + have hdouble : toGM (PM.padd st.2 st.2) = (2 ^ (i + 1) : ℕ) • A := by + rw [toGM_padd hv2 hv2, ha2, ← two_nsmul, smul_smul, pow_succ] + ring_nf + refine ⟨?_, WV_padd hv2 hv2, ?_, hdouble⟩ + · split + · exact WV_padd hv1 hv2 + · exact hv1 + · split + · rename_i hodd + rw [hbit] at hodd + rw [toGM_padd hv1 hv2, ha1, ha2, ← add_nsmul, hmod, hodd, Nat.mul_one] + · rename_i heven + rw [hbit] at heven + have h0 : n / 2 ^ i % 2 = 0 := by omega + rw [ha1, hmod, h0, Nat.mul_zero, Nat.add_zero] + +/-- **The kernel ladder computes `n • ·`** in the affine group, for scalars below `2 ^ 256` +(all Pasta scalars). Statement shape mirrors `PVes.pnsmulFast_spec`. -/ +theorem pnsmulM_spec {p : PM} (hwf : WFP p) (hp : Valid (toPVesM p)) (n : ℕ) (hn : n < 2 ^ 256) : + WFP (PM.pnsmul n p) ∧ Valid (toPVesM (PM.pnsmul n p)) ∧ + toAffine (toPVesM (PM.pnsmul n p)) = n • toAffine (toPVesM p) := by + have base : ∀ m : ℕ, m ≤ 256 → + LadderInvM (toGM p) n m ((List.range m).foldl + (fun (st : PM × PM) i => + (if (n >>> i) &&& 1 = 1 then PM.padd st.1 st.2 else st.1, PM.padd st.2 st.2)) + (PM.pid, p)) := by + intro m + induction m with + | zero => + intro _ + simp only [List.range_zero, List.foldl_nil] + refine ⟨WV_pid, ⟨hwf, hp⟩, ?_, ?_⟩ + · rw [toGM_pid, pow_zero, Nat.mod_one, zero_nsmul] + · rw [pow_zero, one_nsmul] + | succ k ih => + intro hk + rw [List.range_succ, List.foldl_append] + exact ladderM_step (ih (by omega)) + obtain ⟨⟨hw, hv⟩, -, hval, -⟩ := base 256 le_rfl + rw [Nat.mod_eq_of_lt hn] at hval + exact ⟨hw, hv, hval⟩ + +/-! ## Arrays of kernel points -/ + +private theorem forall₂_getElem {α β : Type} {R : α → β → Prop} : + ∀ {l₁ : List α} {l₂ : List β}, List.Forall₂ R l₁ l₂ → + ∀ (i : ℕ) (h₁ : i < l₁.length) (h₂ : i < l₂.length), R l₁[i] l₂[i] := by + intro l₁ l₂ h + induction h with + | nil => intro i h₁ _; simp at h₁ + | cons hab _ ih => + intro i h₁ h₂ + cases i with + | zero => exact hab + | succ k => exact ih k (by simpa using h₁) (by simpa using h₂) + +private theorem forall₂_set {α β : Type} {R : α → β → Prop} {x : α} {y : β} (hxy : R x y) : + ∀ {l₁ : List α} {l₂ : List β}, List.Forall₂ R l₁ l₂ → + ∀ i, List.Forall₂ R (l₁.set i x) (l₂.set i y) := by + intro l₁ l₂ h + induction h with + | nil => intro i; cases i <;> exact List.Forall₂.nil + | cons hab hl ih => + intro i + cases i with + | zero => exact List.Forall₂.cons hxy hl + | succ k => exact List.Forall₂.cons hab (ih k) + +private theorem forall₂_modify {α β : Type} {R : α → β → Prop} {f : α → α} {g : β → β} + (hfg : ∀ a b, R a b → R (f a) (g b)) : + ∀ {l₁ : List α} {l₂ : List β}, List.Forall₂ R l₁ l₂ → + ∀ i, List.Forall₂ R (l₁.modify i f) (l₂.modify i g) := by + intro l₁ l₂ h + induction h with + | nil => intro i; cases i <;> exact List.Forall₂.nil + | cons hab hl ih => + intro i + cases i with + | zero => exact List.Forall₂.cons (hfg _ _ hab) hl + | succ k => exact List.Forall₂.cons hab (ih k) + +private theorem forall₂_replicate {α β : Type} {R : α → β → Prop} {x : α} {y : β} (h : R x y) : + ∀ n, List.Forall₂ R (List.replicate n x) (List.replicate n y) + | 0 => List.Forall₂.nil + | n + 1 => List.Forall₂.cons h (forall₂_replicate h n) + +private theorem forall₂_map_eq {α β γ : Type} {R : α → β → Prop} {f : α → γ} {g : β → γ} + (hfg : ∀ a b, R a b → f a = g b) : + ∀ {l₁ : List α} {l₂ : List β}, List.Forall₂ R l₁ l₂ → l₁.map f = l₂.map g := by + intro l₁ l₂ h + induction h with + | nil => rfl + | cons hab _ ih => rw [List.map_cons, List.map_cons, hfg _ _ hab, ih] + +private theorem forall₂_map_self {α β : Type} {R : α → β → Prop} {P : α → Prop} {f : α → β} + (hf : ∀ a, P a → R a (f a)) : + ∀ (l : List α), (∀ a ∈ l, P a) → List.Forall₂ R l (l.map f) + | [], _ => List.Forall₂.nil + | a :: l, h => + List.Forall₂.cons (hf a (h a (by simp))) + (forall₂_map_self hf l fun b hb => h b (by simp [hb])) + +/-- `Array.get!`/`set!` need a default on the `PVes` side as well as on `PM`'s; the projective +identity is the right one on both, so an out-of-range read stays a corresponding pair. -/ +instance : Inhabited PVes := ⟨PVes.pid⟩ + +/-- The array-level correspondence: cellwise `RM`. -/ +def RA (a : Array PM) (b : Array PVes) : Prop := List.Forall₂ RM a.toList b.toList + +theorem RA.size {a : Array PM} {b : Array PVes} (h : RA a b) : a.size = b.size := by + have := List.Forall₂.length_eq h + rwa [Array.length_toList, Array.length_toList] at this + +theorem RA.get {a : Array PM} {b : Array PVes} (h : RA a b) (k : ℕ) : RM a[k]! b[k]! := by + by_cases hk : k < a.size + · have hk' : k < b.size := h.size ▸ hk + rw [getElem!_pos a k hk, getElem!_pos b k hk'] + have := forall₂_getElem h k (by rwa [Array.length_toList]) (by rwa [Array.length_toList]) + rwa [Array.getElem_toList, Array.getElem_toList] at this + · have hk' : ¬ k < b.size := by rw [← h.size]; exact hk + rw [getElem!_neg a k hk, getElem!_neg b k hk'] + exact RM_pid + +theorem RA.set {a : Array PM} {b : Array PVes} (h : RA a b) (k : ℕ) {p : PM} {P : PVes} + (hp : RM p P) : RA (a.set! k p) (b.set! k P) := by + simp only [RA, Array.set!_eq_setIfInBounds, Array.toList_setIfInBounds] + exact forall₂_set hp h k + +theorem RA.modify {a : Array PM} {b : Array PVes} (h : RA a b) (k : ℕ) + {f : PM → PM} {g : PVes → PVes} (hfg : ∀ p P, RM p P → RM (f p) (g P)) : + RA (a.modify k f) (b.modify k g) := by + simp only [RA, Array.toList_modify] + exact forall₂_modify hfg h k + +theorem RA.map_toG {a : Array PM} {b : Array PVes} (h : RA a b) : + a.map toGM = b.map toAffine := by + have := forall₂_map_eq (f := toGM) (g := toAffine) (fun _ _ hr => RM_toG hr) h + rw [← Array.toList_map, ← Array.toList_map] at this + exact Array.toList_inj.mp this + +theorem RA_self {a : Array PM} (h : ∀ p ∈ a, WFP p) : RA a (a.map toPVesM) := by + simp only [RA, Array.toList_map] + exact forall₂_map_self (P := WFP) (fun p hp => RM_self hp) a.toList + (fun p hp => h p (by simpa using hp)) + +/-! ## The windowed Pippenger MSM + +Window values transport by `RM` onto `MsmProj.pwindowValueFast` and get their affine meaning from +its spec. The Horner recombination is the one part that does not mirror `MsmProj` (`c` doublings +instead of `pnsmulFast (2 ^ c)`), so it goes through `toGM` directly, and +`hornerList_windows_eq_msm` reconciles the kernel's fixed `⌈256 / c⌉` windows with +`Msm.numWindows`. -/ + +private theorem RA_scatterStep {a : Array PM} {b : Array PVes} (h : RA a b) {t : ℕ × PM} + (ht : WFP t.2) : + RA (PM.scatterStep a t) (MsmProj.pscatterStep b (t.1, toPVesM t.2)) := by + simp only [PM.scatterStep, MsmProj.pscatterStep] + by_cases h0 : t.1 = 0 + · rw [if_pos h0, if_pos h0] + exact h + · rw [if_neg h0, if_neg h0] + exact h.modify _ (fun _ _ hpn => RM_padd hpn (RM_self ht)) + +private theorem RA_bucketScatter (base : ℕ) (dp : List (ℕ × PM)) (h : ∀ t ∈ dp, WFP t.2) : + RA (PM.bucketScatter base dp) + (MsmProj.pbucketScatter base (dp.map fun t => (t.1, toPVesM t.2))) := by + simp only [PM.bucketScatter, MsmProj.pbucketScatter, List.foldl_map] + refine foldl_rel₂ RA (fun t : ℕ × PM => WFP t.2) dp _ _ + (fun s s' x hx hs => RA_scatterStep (t := x) hs hx) h _ _ ?_ + simp only [RA, Array.toList_replicate] + exact forall₂_replicate RM_pid _ + +private theorem RM2_foldr_accStep : + ∀ {L : List PM} {L' : List PVes}, List.Forall₂ RM L L' → + RM2 (L.foldr PM.accStep (PM.pid, PM.pid)) + (L'.foldr MsmProj.paccStep (PVes.pid, PVes.pid)) := by + intro L L' h + induction h with + | nil => exact ⟨RM_pid, RM_pid⟩ + | cons hab _ ih => + simp only [List.foldr_cons, PM.accStep, MsmProj.paccStep] + exact ⟨RM_padd ih.1 hab, RM_padd ih.2 (RM_padd ih.1 hab)⟩ + +private theorem RM_windowValue (base i : ℕ) (terms : List (ℕ × PM)) + (h : ∀ t ∈ terms, WFP t.2) : + RM (PM.windowValue base i terms) + (MsmProj.pwindowValueFast base i (terms.map fun t => (t.1, toPVesM t.2))) := by + have hdp : ∀ t ∈ terms.map (fun t : ℕ × PM => (t.1 / base ^ i % base, t.2)), WFP t.2 := by + intro t ht + rw [List.mem_map] at ht + obtain ⟨s, hs, rfl⟩ := ht + exact h s hs + have hb := RA_bucketScatter base (terms.map fun t : ℕ × PM => (t.1 / base ^ i % base, t.2)) hdp + have halign : (terms.map fun t : ℕ × PM => (t.1 / base ^ i % base, t.2)).map + (fun t : ℕ × PM => (t.1, toPVesM t.2)) + = MsmProj.pdpOf base i (terms.map fun t : ℕ × PM => (t.1, toPVesM t.2)) := by + simp only [MsmProj.pdpOf, Msm.digit, List.map_map, Function.comp_def] + rw [halign] at hb + exact (RM2_foldr_accStep hb).2 + +/-- The `c`-fold doubling is `2 ^ c • ·` in the affine group. -/ +private theorem pdoublingsM_spec {p : PM} (hp : WV p) (c : ℕ) : + WV (PM.pdoublings c p) ∧ toGM (PM.pdoublings c p) = (2 ^ c : ℕ) • toGM p := by + induction c with + | zero => + rw [show PM.pdoublings 0 p = p from rfl] + exact ⟨hp, by rw [pow_zero, one_nsmul]⟩ + | succ k ih => + obtain ⟨hv, he⟩ := ih + have hstep : PM.pdoublings (k + 1) p = PM.padd (PM.pdoublings k p) (PM.pdoublings k p) := by + rw [PM.pdoublings, List.range_succ, List.foldl_append]; rfl + rw [hstep] + refine ⟨WV_padd hv hv, ?_⟩ + rw [toGM_padd hv hv, he, ← two_nsmul, smul_smul, pow_succ] + ring_nf + +/-- The kernel's Horner recombination across windows is `Msm.hornerList` after `toGM`. -/ +private theorem hfoldM_spec (c : ℕ) (vals : List PM) (h : ∀ v ∈ vals, WV v) : + WV (vals.foldr (fun v acc => PM.padd (PM.pdoublings c acc) v) PM.pid) ∧ + toGM (vals.foldr (fun v acc => PM.padd (PM.pdoublings c acc) v) PM.pid) + = Msm.hornerList (2 ^ c) (vals.map toGM) := by + induction vals with + | nil => + refine ⟨by rw [List.foldr_nil]; exact WV_pid, ?_⟩ + rw [List.foldr_nil, toGM_pid, List.map_nil, Msm.hornerList, List.foldr_nil] + | cons v xs ih => + have hv : WV v := h v (by simp) + obtain ⟨hacc, heq⟩ := ih fun w hw => h w (by simp [hw]) + obtain ⟨hdv, hde⟩ := pdoublingsM_spec hacc c + rw [List.foldr_cons] + refine ⟨WV_padd hdv hv, ?_⟩ + rw [toGM_padd hdv hv, hde, heq, List.map_cons] + simp only [Msm.hornerList, List.foldr_cons] + +/-- **Horner recombination of `W` windows is the naive MSM**, whenever `W` base-`2 ^ c` digits +cover every scalar. This is `Msm.pippenger_eq_msm` with the window count freed from +`Msm.numWindows` (the kernel fixes it at `⌈256 / c⌉`). -/ +private theorem hornerList_windows_eq_msm (c W : ℕ) (terms : List (ℕ × G)) + (hW : ∀ t ∈ terms, t.1 < (2 ^ c) ^ W) : + Msm.hornerList (2 ^ c) ((List.range W).map fun i => Msm.windowValue (2 ^ c) i terms) + = (terms.map fun t => t.1 • t.2).sum := by + have hb0 : 0 < 2 ^ c := by positivity + have hpip : Msm.hornerList (2 ^ c) + ((List.range W).map fun i => Msm.windowValue (2 ^ c) i terms) + = ∑ i ∈ Finset.range W, + (2 ^ c) ^ i • (terms.map fun t => Msm.digit (2 ^ c) i t.1 • t.2).sum := by + rw [Msm.hornerList_eq, List.length_map, List.length_range] + refine Finset.sum_congr rfl fun k hk => ?_ + rw [Finset.mem_range] at hk + rw [Msm.getD_map_range _ _ _ hk, Msm.windowValue_eq (2 ^ c) k hb0] + rw [hpip] + have e1 : (terms.map fun t => t.1 • t.2) + = terms.map fun t => + ∑ i ∈ Finset.range W, Msm.digit (2 ^ c) i t.1 • ((2 ^ c) ^ i • t.2) := + List.map_congr_left fun t ht => Msm.smul_eq_sum_digits (2 ^ c) hb0 W t.1 t.2 (hW t ht) + rw [e1, Msm.list_sum_finset_sum] + refine Finset.sum_congr rfl fun i _ => ?_ + rw [Msm.smul_list_sum, List.map_map] + refine congrArg List.sum (List.map_congr_left fun t _ => ?_) + simp only [Function.comp_apply] + rw [smul_comm] + +/-- **The kernel's windowed Pippenger MSM is the proven affine Pippenger.** Statement shape +mirrors `MsmProj.pippengerProjScatter_eq`. -/ +theorem msmM_spec (c : ℕ) (hc : 0 < c) (terms : List (ℕ × PM)) + (hwf : ∀ t ∈ terms, WFP t.2) (hv : ∀ t ∈ terms, Valid (toPVesM t.2)) + (hn : ∀ t ∈ terms, t.1 < 2 ^ 256) : + toAffine (toPVesM (PM.msm c terms)) + = Msm.pippenger c (terms.map fun t => (t.1, toAffine (toPVesM t.2))) := by + set W := (256 + c - 1) / c with hWdef + set pterms := terms.map fun t => (t.1, toPVesM t.2) with hpterms + set aterms := terms.map fun t => (t.1, toAffine (toPVesM t.2)) with haterms + have hptv : ∀ p ∈ pterms, Valid p.2 := by + intro p hp + rw [hpterms, List.mem_map] at hp + obtain ⟨t, ht, rfl⟩ := hp + exact hv t ht + have hmapaff : pterms.map (fun t => (t.1, toAffine t.2)) = aterms := by + rw [hpterms, haterms, List.map_map] + rfl + -- each window value: well formed, representable, and affinely the `Msm` window value + have hwin : ∀ i, WV (PM.windowValue (2 ^ c) i terms) ∧ + toGM (PM.windowValue (2 ^ c) i terms) = Msm.windowValue (2 ^ c) i aterms := by + intro i + have hR := RM_windowValue (2 ^ c) i terms hwf + rw [← hpterms] at hR + obtain ⟨hvw, hew⟩ := MsmProj.pwindowValueFast_spec (2 ^ c) i pterms hptv + exact ⟨RM_valid hR hvw, by rw [RM_toG hR, hew, hmapaff]⟩ + have hvals : ∀ v ∈ (List.range W).map (fun i => PM.windowValue (2 ^ c) i terms), WV v := by + intro v hvm + rw [List.mem_map] at hvm + obtain ⟨i, -, rfl⟩ := hvm + exact (hwin i).1 + have hmsm : PM.msm c terms + = ((List.range W).map fun i => PM.windowValue (2 ^ c) i terms).foldr + (fun v acc => PM.padd (PM.pdoublings c acc) v) PM.pid := rfl + have hmapwin : ((List.range W).map fun i => PM.windowValue (2 ^ c) i terms).map toGM + = (List.range W).map fun i => Msm.windowValue (2 ^ c) i aterms := by + simp only [List.map_map, Function.comp_def] + exact List.map_congr_left fun i _ => (hwin i).2 + show toGM (PM.msm c terms) = _ + rw [hmsm, (hfoldM_spec c _ hvals).2, hmapwin] + have hcW : 256 ≤ c * W := by + have h1 : c * W + (256 + c - 1) % c = 256 + c - 1 := by + rw [hWdef]; exact Nat.div_add_mod (256 + c - 1) c + have h2 : (256 + c - 1) % c < c := Nat.mod_lt _ hc + obtain ⟨x, hx⟩ : ∃ x, c * W = x := ⟨_, rfl⟩ + rw [hx] at h1 ⊢ + omega + have hbound : ∀ t ∈ aterms, t.1 < (2 ^ c) ^ W := by + intro t ht + rw [haterms, List.mem_map] at ht + obtain ⟨s, hs, rfl⟩ := ht + calc s.1 < 2 ^ 256 := hn s hs + _ ≤ (2 ^ c) ^ W := by rw [← pow_mul]; exact Nat.pow_le_pow_right (by norm_num) hcW + rw [hornerList_windows_eq_msm c W aterms hbound, ← Msm.pippenger_eq_msm c hc] + +end CompElliptic.Curves.Pasta.Fast.ProjectiveMont diff --git a/CompElliptic/TrustBoundary.lean b/CompElliptic/TrustBoundary.lean index 18a1439..cc1c508 100644 --- a/CompElliptic/TrustBoundary.lean +++ b/CompElliptic/TrustBoundary.lean @@ -5,6 +5,9 @@ as described in the files LICENSE-APACHE and LICENSE-MIT. Authors: Daira-Emma Hopwood -/ import CompElliptic.Curves.PastaOrder +import CompElliptic.Curves.Pasta.Fast.Projective +import CompElliptic.Curves.Pasta.Fast.Msm +import CompElliptic.Curves.Pasta.Fast.ProjectiveMontEquiv import CompElliptic.Fields.Sqrt import CompElliptic.Meta.AxiomCheck @@ -68,3 +71,19 @@ assert_axioms CompElliptic.Curves.Pasta.Pallas.q_nsmul_Gpt +native assert_axioms CompElliptic.Curves.Pasta.Vesta.p_nsmul_Gpt +native assert_axioms CompElliptic.Curves.Pasta.Pallas.card_eq +native assert_axioms CompElliptic.Curves.Pasta.Vesta.card_eq +native + +/-! ## Fast Vesta arithmetic — proven against the affine group law, standard axioms only + +The fast tier's headline correctness theorems: the RCB projective addition is the affine +group law (`toAffine_padd`), the projective scalar ladder and the windowed Pippenger MSM +compute the operations they replace, the raw-`ℕ` spelling is `padd`, and the Montgomery +kernel computes the same schedules. Completeness routes through `no_onCurve_y_zero` +(a kernel `decide`), so nothing here reaches `native_decide`. -/ + +assert_axioms CompElliptic.Curves.Pasta.Fast.Projective.PVes.toAffine_padd +assert_axioms CompElliptic.Curves.Pasta.Fast.Projective.PVes.smulFast_eq +assert_axioms CompElliptic.Curves.Pasta.Fast.Projective.PVes.padd_eq_paddFast +assert_axioms CompElliptic.Curves.Pasta.Fast.Msm.pippengerFastPar_eq_msm +assert_axioms CompElliptic.Curves.Pasta.Fast.Msm.commitLagrangeFastWith_eq +assert_axioms CompElliptic.Curves.Pasta.Fast.ProjectiveMont.pnsmulM_spec +assert_axioms CompElliptic.Curves.Pasta.Fast.ProjectiveMont.msmM_spec diff --git a/CompElliptic/Vendor/CompPoly/Montgomery/Basic.lean b/CompElliptic/Vendor/CompPoly/Montgomery/Basic.lean new file mode 100644 index 0000000..47336c8 --- /dev/null +++ b/CompElliptic/Vendor/CompPoly/Montgomery/Basic.lean @@ -0,0 +1,114 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Valerii Huhnin, Georgios Raikos +-/ + +import Mathlib.Data.Nat.ModEq +import Mathlib.Data.ZMod.Basic +import Mathlib.Algebra.Field.ZMod +import Mathlib.Tactic.Ring + +/-! +# Montgomery Reduction + +Radix-generic specification and correctness lemmas for single-word Montgomery reduction. +Word-specific implementations refine these results in sibling modules. +-/ + +namespace Montgomery + +/-- Natural-number Montgomery reduction used to specify the native-word reducer. -/ +def reduceNat (R p negInv x : ℕ) : ℕ := + let m := (x % R * negInv) % R + let u := (x + m * p) / R + if u < p then u else u - p + +/-- The quotient before the final conditional subtraction in Montgomery reduction. -/ +def reduceNatQuotient (R p negInv x : ℕ) : ℕ := + let m := (x % R * negInv) % R + (x + m * p) / R + +/-- The pre-subtraction quotient is below twice the modulus. -/ +theorem reduceNatQuotient_lt_two_mul (R p negInv x : ℕ) + (hR : 0 < R) (hp : 0 < p) (hx : x < p * R) : + reduceNatQuotient R p negInv x < 2 * p := by + let m := x % R * negInv % R + have hm_lt : m < R := Nat.mod_lt _ hR + change (x + m * p) / R < 2 * p + rw [Nat.div_lt_iff_lt_mul] + · have hprod_lt : m * p < R * p := Nat.mul_lt_mul_of_pos_right hm_lt hp + have hprod_lt' : m * p < p * R := by + simpa only [Nat.mul_comm] using hprod_lt + calc + x + m * p < p * R + p * R := Nat.add_lt_add hx hprod_lt' + _ = 2 * p * R := by ring + · exact hR + +/-- Montgomery reduction returns a canonical representative. -/ +theorem reduceNat_lt (R p negInv x : ℕ) + (hR : 0 < R) (hp : 0 < p) (hx : x < p * R) : + reduceNat R p negInv x < p := by + change (if reduceNatQuotient R p negInv x < p then + reduceNatQuotient R p negInv x else reduceNatQuotient R p negInv x - p) < p + have hu := reduceNatQuotient_lt_two_mul R p negInv x hR hp hx + by_cases h : reduceNatQuotient R p negInv x < p + · rw [if_pos h] + exact h + · rw [if_neg h] + omega + +/-- The Montgomery divisibility identity: if `(negInv * p) % R = R - 1` (i.e. +`negInv = -p⁻¹ mod R`), then `R ∣ x + ((x mod R)·negInv mod R)·p` for every `x`. -/ +theorem dvd_add (R p negInv : ℕ) (hR : 0 < R) + (hnegInv : negInv * p % R = R - 1) (x : ℕ) : + R ∣ x + ((x % R * negInv) % R) * p := by + rw [Nat.dvd_iff_mod_eq_zero] + rw [Nat.add_mod, Nat.mod_mul_mod (x % R * negInv) p R, Nat.mul_assoc] + rw [← Nat.mul_mod_mod, hnegInv, Nat.add_mod_mod, add_comm, ← Nat.mul_add_one] + rw [show R - 1 + 1 = R by omega, Nat.mul_mod_left] + +/-- The pre-subtraction quotient represents multiplication by `R⁻¹` in `ZMod p`. -/ +theorem reduceNatQuotient_cast (R p negInv : ℕ) [Fact (Nat.Prime p)] (hR : 0 < R) + (hnegInv : negInv * p % R = R - 1) (hRne : (R : ZMod p) ≠ 0) (x : ℕ) : + (reduceNatQuotient R p negInv x : ZMod p) = (x : ZMod p) * (R : ZMod p)⁻¹ := by + let m := x % R * negInv % R + let u := (x + m * p) / R + change (u : ZMod p) = (x : ZMod p) * (R : ZMod p)⁻¹ + have hdiv : R ∣ x + m * p := by + simpa [m] using dvd_add R p negInv hR hnegInv x + have hu_mul : u * R = x + m * p := Nat.div_mul_cancel hdiv + have hcast_mul : (u : ZMod p) * (R : ZMod p) = (x : ZMod p) := by + rw [← Nat.cast_mul, hu_mul, Nat.cast_add, Nat.cast_mul] + simp + rw [← hcast_mul] + exact (mul_inv_cancel_right₀ hRne (u : ZMod p)).symm + +/-- Montgomery reduction represents multiplication by `R⁻¹` in `ZMod p`. -/ +theorem reduceNat_cast (R p negInv : ℕ) [Fact (Nat.Prime p)] (hR : 0 < R) + (hnegInv : negInv * p % R = R - 1) (hRne : (R : ZMod p) ≠ 0) (x : ℕ) : + (reduceNat R p negInv x : ZMod p) = (x : ZMod p) * (R : ZMod p)⁻¹ := by + let m := x % R * negInv % R + let u := (x + m * p) / R + have hu_cast : (u : ZMod p) = (x : ZMod p) * (R : ZMod p)⁻¹ := by + simpa only [reduceNatQuotient, m, u] using + reduceNatQuotient_cast R p negInv hR hnegInv hRne x + change ((if u < p then u else u - p : ℕ) : ZMod p) = + (x : ZMod p) * (R : ZMod p)⁻¹ + by_cases hu : u < p + · rw [if_pos hu] + exact hu_cast + · have hfield : ((u - p : ℕ) : ZMod p) = (u : ZMod p) := by + rw [Nat.cast_sub (Nat.le_of_not_gt hu)] + simp + rw [if_neg hu] + rw [hfield] + exact hu_cast + +/-- Two naturals below `p` are equal once their `ZMod p` casts agree. -/ +theorem natCast_inj_of_lt {p a b : ℕ} (h : (a : ZMod p) = (b : ZMod p)) + (ha : a < p) (hb : b < p) : a = b := by + rw [ZMod.natCast_eq_natCast_iff] at h + exact Nat.ModEq.eq_of_lt_of_lt h ha hb + +end Montgomery diff --git a/CompElliptic/Vendor/CompPoly/Montgomery/Native64x8.lean b/CompElliptic/Vendor/CompPoly/Montgomery/Native64x8.lean new file mode 100644 index 0000000..be2a389 --- /dev/null +++ b/CompElliptic/Vendor/CompPoly/Montgomery/Native64x8.lean @@ -0,0 +1,429 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Gregor Mitscha-Baude +-/ +import CompElliptic.Vendor.CompPoly.Montgomery.Basic +import CompElliptic.Vendor.CompPoly.Montgomery.Native64x8Defs + +/-! +# Eight-limb Montgomery arithmetic: correctness of the raw operations + +Vendored from CompPoly branch `fast_multilimb_fields`; see `README.md` in this directory. +The runtime definitions live in `CompElliptic.Vendor.CompPoly.Montgomery.Native64x8Defs`. +-/ + +namespace Montgomery + +namespace Native64x8 + +/-! ## Word helpers -/ + +/-! ### Word-level specifications -/ + +private theorem land_mask32 (x : ℕ) : x &&& 4294967295 = x % 4294967296 := by + have h : (4294967295 : ℕ) = 2 ^ 32 - 1 := by norm_num + have h2 : (4294967296 : ℕ) = 2 ^ 32 := by norm_num + rw [h, h2, Nat.and_two_pow_sub_one_eq_mod] + +private theorem toNat_zero : (0 : UInt64).toNat = 0 := rfl + +/-- Add-with-carry, existential form: the low limb and the carry-out are named as fresh +naturals together with their defining equations, so that a chain of them is a linear system +over plain variables. -/ +theorem adc_spec (x y c : UInt64) (hx : x.toNat < 2 ^ 32) (hy : y.toNat < 2 ^ 32) + (hc : c.toNat ≤ 1) : + ∃ s co : ℕ, (adcLo x y c).toNat = s ∧ (adcCo x y c).toNat = co ∧ + s + 2 ^ 32 * co = x.toNat + y.toNat + c.toNat ∧ co ≤ 1 ∧ s < 2 ^ 32 := by + refine ⟨_, _, rfl, rfl, ?_⟩ + simp only [adcLo, adcCo, mask, UInt64.toNat_and, UInt64.toNat_add, + UInt64.toNat_shiftRight, UInt64.toNat_ofNat, Nat.reduceMod, Nat.reducePow, land_mask32] + omega + +/-- Add-with-carry when the first summand is one bit wider than a limb. -/ +theorem adc_spec_wide (x y c : UInt64) (hx : x.toNat < 2 ^ 33) (hy : y.toNat < 2 ^ 32) + (hc : c.toNat ≤ 1) : + ∃ s co : ℕ, (adcLo x y c).toNat = s ∧ (adcCo x y c).toNat = co ∧ + s + 2 ^ 32 * co = x.toNat + y.toNat + c.toNat ∧ co ≤ 2 ∧ s < 2 ^ 32 := by + refine ⟨_, _, rfl, rfl, ?_⟩ + simp only [adcLo, adcCo, mask, UInt64.toNat_and, UInt64.toNat_add, + UInt64.toNat_shiftRight, UInt64.toNat_ofNat, Nat.reduceMod, Nat.reducePow, land_mask32] + omega + +/-- Subtract-with-borrow, existential form. -/ +theorem sbb_spec (x y b : UInt64) (hx : x.toNat < 2 ^ 32) (hy : y.toNat < 2 ^ 32) + (hb : b.toNat ≤ 1) : + ∃ m bo : ℕ, (sbbLo x y b).toNat = m ∧ (sbbBo x y b).toNat = bo ∧ + m + y.toNat + b.toNat = x.toNat + 2 ^ 32 * bo ∧ bo ≤ 1 ∧ m < 2 ^ 32 := by + refine ⟨_, _, rfl, rfl, ?_⟩ + simp only [sbbLo, sbbBo, mask, UInt64.toNat_and, UInt64.toNat_sub, + UInt64.toNat_shiftRight, UInt64.toNat_ofNat, Nat.reduceMod, Nat.reducePow, land_mask32] + omega + +/-- Multiply-accumulate, existential form. The accumulator, both factors and the carry-in +all fit in 32 bits, so `t + x * y + c` cannot overflow a 64-bit word. -/ +theorem mac_spec (t x y c : UInt64) (ht : t.toNat < 2 ^ 32) (hx : x.toNat < 2 ^ 32) + (hy : y.toNat < 2 ^ 32) (hc : c.toNat < 2 ^ 32) : + ∃ lo hi : ℕ, (macLo t x y c).toNat = lo ∧ (macHi t x y c).toNat = hi ∧ + lo + 2 ^ 32 * hi = t.toNat + x.toNat * y.toNat + c.toNat ∧ hi < 2 ^ 32 ∧ + lo < 2 ^ 32 := by + have hp : x.toNat * y.toNat ≤ (2 ^ 32 - 1) * (2 ^ 32 - 1) := + Nat.mul_le_mul (by omega) (by omega) + refine ⟨_, _, rfl, rfl, ?_⟩ + simp only [macLo, macHi, mask, UInt64.toNat_and, UInt64.toNat_add, UInt64.toNat_mul, + UInt64.toNat_shiftRight, UInt64.toNat_ofNat, Nat.reduceMod, Nat.reducePow, land_mask32] + generalize x.toNat * y.toNat = p at hp ⊢ + omega + +/-- The Montgomery multiplier agrees with its natural-number specification. -/ +theorem montM_toNat (s negInv : UInt64) (hs : s.toNat < 2 ^ 32) (hn : negInv.toNat < 2 ^ 32) : + (montM s negInv).toNat = s.toNat * negInv.toNat % 2 ^ 32 := by + have hp : s.toNat * negInv.toNat ≤ (2 ^ 32 - 1) * (2 ^ 32 - 1) := + Nat.mul_le_mul (by omega) (by omega) + simp only [montM, mask, UInt64.toNat_and, UInt64.toNat_mul, UInt64.toNat_ofNat, + Nat.reducePow, Nat.reduceMod, land_mask32, + Nat.mod_eq_of_lt (show s.toNat < 4294967296 by omega)] + generalize s.toNat * negInv.toNat = p at hp ⊢ + omega + +theorem adcLo_lt (x y c : UInt64) : (adcLo x y c).toNat < 2 ^ 32 := by + simp only [adcLo, mask, UInt64.toNat_and, UInt64.toNat_ofNat, Nat.reducePow, + Nat.reduceMod, land_mask32] + omega + +theorem sbbLo_lt (x y b : UInt64) : (sbbLo x y b).toNat < 2 ^ 32 := by + simp only [sbbLo, mask, UInt64.toNat_and, UInt64.toNat_ofNat, Nat.reducePow, + Nat.reduceMod, land_mask32] + omega + +theorem macLo_lt (t x y c : UInt64) : (macLo t x y c).toNat < 2 ^ 32 := by + simp only [macLo, mask, UInt64.toNat_and, UInt64.toNat_ofNat, Nat.reducePow, + Nat.reduceMod, land_mask32] + omega + +theorem montM_lt (s negInv : UInt64) : (montM s negInv).toNat < 2 ^ 32 := by + simp only [montM, mask, UInt64.toNat_and, UInt64.toNat_ofNat, Nat.reducePow, + Nat.reduceMod, land_mask32] + omega + +/-! ## Eight-limb values -/ + +namespace Limbs8 + +theorem zero_bounded : zero.Bounded := by + simp only [Bounded, zero, toNat_zero] + norm_num + +theorem zero_toNat : zero.toNat = 0 := by + simp only [toNat, zero, toNat_zero] + +theorem one_bounded : one.Bounded := by decide + +theorem one_toNat : one.toNat = 1 := by decide + +theorem ofNat_bounded (n : ℕ) : (ofNat n).Bounded := by + simp only [Bounded, ofNat, UInt64.toNat_ofNat', land_mask32] + omega + +theorem ofNat_toNat (n : ℕ) : (ofNat n).toNat = n % 2 ^ 256 := by + simp only [toNat, ofNat, UInt64.toNat_ofNat', land_mask32] + omega + +end Limbs8 + +/-- A recomposed value with 32-bit limbs stays below `2 ^ 256`. -/ +theorem sum8_lt {a0 a1 a2 a3 a4 a5 a6 a7 : ℕ} (h0 : a0 < 2 ^ 32) (h1 : a1 < 2 ^ 32) + (h2 : a2 < 2 ^ 32) (h3 : a3 < 2 ^ 32) (h4 : a4 < 2 ^ 32) (h5 : a5 < 2 ^ 32) + (h6 : a6 < 2 ^ 32) (h7 : a7 < 2 ^ 32) : + a0 + 2 ^ 32 * a1 + 2 ^ 64 * a2 + 2 ^ 96 * a3 + 2 ^ 128 * a4 + 2 ^ 160 * a5 + 2 ^ 192 * a6 + + 2 ^ 224 * a7 < 2 ^ 256 := by + omega + +/-- Bounded limb vectors are determined by their value. -/ +theorem Limbs8.ext_of_toNat {x y : Limbs8} (hx : x.Bounded) (hy : y.Bounded) + (h : x.toNat = y.toNat) : x = y := by + obtain ⟨hx0, hx1, hx2, hx3, hx4, hx5, hx6, hx7⟩ := hx + obtain ⟨hy0, hy1, hy2, hy3, hy4, hy5, hy6, hy7⟩ := hy + simp only [Limbs8.toNat] at h + have e0 : x.l0.toNat = y.l0.toNat := by omega + have e1 : x.l1.toNat = y.l1.toNat := by omega + have e2 : x.l2.toNat = y.l2.toNat := by omega + have e3 : x.l3.toNat = y.l3.toNat := by omega + have e4 : x.l4.toNat = y.l4.toNat := by omega + have e5 : x.l5.toNat = y.l5.toNat := by omega + have e6 : x.l6.toNat = y.l6.toNat := by omega + have e7 : x.l7.toNat = y.l7.toNat := by omega + obtain ⟨x0, x1, x2, x3, x4, x5, x6, x7⟩ := x + obtain ⟨y0, y1, y2, y3, y4, y5, y6, y7⟩ := y + simp only [Limbs8.mk.injEq] + exact ⟨UInt64.toNat_inj.mp e0, UInt64.toNat_inj.mp e1, UInt64.toNat_inj.mp e2, + UInt64.toNat_inj.mp e3, UInt64.toNat_inj.mp e4, UInt64.toNat_inj.mp e5, + UInt64.toNat_inj.mp e6, UInt64.toNat_inj.mp e7⟩ + +/-- A bounded eight-limb value is below `2 ^ 256`. -/ +theorem Limbs8.toNat_lt {x : Limbs8} (hx : x.Bounded) : x.toNat < 2 ^ 256 := by + obtain ⟨h0, h1, h2, h3, h4, h5, h6, h7⟩ := hx + exact sum8_lt h0 h1 h2 h3 h4 h5 h6 h7 + +/-! ### Chain lemmas + +Telescoping identities for an eight-step carry or borrow chain. They are stated over plain +naturals and the `y` slots are arbitrary, so `carry_chain_sum` serves both the +add-with-carry chains and the multiply-accumulate chains of CIOS, where `yᵢ` is a limb +product. -/ + +/-- Weighted telescoping of an eight-step carry chain. -/ +theorem carry_chain_sum {x0 x1 x2 x3 x4 x5 x6 x7 y0 y1 y2 y3 y4 y5 y6 y7 : ℕ} + {s0 s1 s2 s3 s4 s5 s6 s7 cin c0 c1 c2 c3 c4 c5 c6 c7 : ℕ} + (e0 : s0 + 2 ^ 32 * c0 = x0 + y0 + cin) + (e1 : s1 + 2 ^ 32 * c1 = x1 + y1 + c0) + (e2 : s2 + 2 ^ 32 * c2 = x2 + y2 + c1) + (e3 : s3 + 2 ^ 32 * c3 = x3 + y3 + c2) + (e4 : s4 + 2 ^ 32 * c4 = x4 + y4 + c3) + (e5 : s5 + 2 ^ 32 * c5 = x5 + y5 + c4) + (e6 : s6 + 2 ^ 32 * c6 = x6 + y6 + c5) + (e7 : s7 + 2 ^ 32 * c7 = x7 + y7 + c6) : + s0 + 2 ^ 32 * s1 + 2 ^ 64 * s2 + 2 ^ 96 * s3 + 2 ^ 128 * s4 + 2 ^ 160 * s5 + 2 ^ 192 * s6 + + 2 ^ 224 * s7 + 2 ^ 256 * c7 = + (x0 + 2 ^ 32 * x1 + 2 ^ 64 * x2 + 2 ^ 96 * x3 + 2 ^ 128 * x4 + 2 ^ 160 * x5 + + 2 ^ 192 * x6 + 2 ^ 224 * x7) + + (y0 + 2 ^ 32 * y1 + 2 ^ 64 * y2 + 2 ^ 96 * y3 + 2 ^ 128 * y4 + 2 ^ 160 * y5 + + 2 ^ 192 * y6 + 2 ^ 224 * y7) + cin := by + omega + +/-- Weighted telescoping of an eight-step borrow chain. -/ +theorem borrow_chain_sum {x0 x1 x2 x3 x4 x5 x6 x7 y0 y1 y2 y3 y4 y5 y6 y7 : ℕ} + {m0 m1 m2 m3 m4 m5 m6 m7 bin b0 b1 b2 b3 b4 b5 b6 b7 : ℕ} + (e0 : m0 + y0 + bin = x0 + 2 ^ 32 * b0) + (e1 : m1 + y1 + b0 = x1 + 2 ^ 32 * b1) + (e2 : m2 + y2 + b1 = x2 + 2 ^ 32 * b2) + (e3 : m3 + y3 + b2 = x3 + 2 ^ 32 * b3) + (e4 : m4 + y4 + b3 = x4 + 2 ^ 32 * b4) + (e5 : m5 + y5 + b4 = x5 + 2 ^ 32 * b5) + (e6 : m6 + y6 + b5 = x6 + 2 ^ 32 * b6) + (e7 : m7 + y7 + b6 = x7 + 2 ^ 32 * b7) : + (m0 + 2 ^ 32 * m1 + 2 ^ 64 * m2 + 2 ^ 96 * m3 + 2 ^ 128 * m4 + 2 ^ 160 * m5 + + 2 ^ 192 * m6 + 2 ^ 224 * m7) + + (y0 + 2 ^ 32 * y1 + 2 ^ 64 * y2 + 2 ^ 96 * y3 + 2 ^ 128 * y4 + 2 ^ 160 * y5 + + 2 ^ 192 * y6 + 2 ^ 224 * y7) + bin = + x0 + 2 ^ 32 * x1 + 2 ^ 64 * x2 + 2 ^ 96 * x3 + 2 ^ 128 * x4 + 2 ^ 160 * x5 + 2 ^ 192 * x6 + + 2 ^ 224 * x7 + 2 ^ 256 * b7 := by + omega + +/-! ### Branch lemmas + +The situations that a borrow chain can produce, each as a standalone linear problem over a +handful of naturals. -/ + +private theorem cond_of_borrow_zero {M Q T : ℕ} (h : M + Q = T) (hM : M < 2 ^ 256) : + M = if T < Q then T else T - Q := by + rw [if_neg (by omega)] + omega + +private theorem cond_of_borrow_one {M Q T : ℕ} (h : M + Q = T + 2 ^ 256) (hM : M < 2 ^ 256) : + T = if T < Q then T else T - Q := by + rw [if_pos (by omega)] + +private theorem cond_eq_mod {T Q : ℕ} (h : T < 2 * Q) : + (if T < Q then T else T - Q) = T % Q := by + by_cases hc : T < Q + · rw [if_pos hc, Nat.mod_eq_of_lt hc] + · rw [if_neg hc] + conv_rhs => rw [show T = T - Q + Q by omega] + rw [Nat.add_mod_right, Nat.mod_eq_of_lt (by omega)] + +private theorem carry_top_zero {S A B Q c : ℕ} (h : S + 2 ^ 256 * c = A + B) + (hS : S < 2 ^ 256) (hA : A < Q) (hB : B < Q) (hQ : 2 * Q < 2 ^ 256) : + S = A + B ∧ S < 2 * Q := by + omega + +private theorem sub_of_borrow_zero {A B Q D : ℕ} (h : D + B = A + 2 ^ 256 * 0) + (hB : B < Q) (hA : A < Q) : D = (A + (Q - B)) % Q := by + conv_rhs => rw [show A + (Q - B) = D + Q by omega] + rw [Nat.add_mod_right, Nat.mod_eq_of_lt (by omega)] + +private theorem sub_of_borrow_one {A B Q D E c : ℕ} (hD : D + B = A + 2 ^ 256 * 1) + (hE : E + 2 ^ 256 * c = D + Q) (hDlt : D < 2 ^ 256) (hElt : E < 2 ^ 256) + (hB : B < Q) (hA : A < Q) (hQ : 2 * Q < 2 ^ 256) : E = (A + (Q - B)) % Q := by + rw [show A + (Q - B) = E by omega, Nat.mod_eq_of_lt (by omega)] + +/-! ## Limbwise addition and subtraction -/ + +theorem addLimbs_bounded (a b : Limbs8) : (addLimbs a b).Bounded := by + simp only [addLimbs, Limbs8.Bounded] + exact ⟨adcLo_lt _ _ _, adcLo_lt _ _ _, adcLo_lt _ _ _, adcLo_lt _ _ _, + adcLo_lt _ _ _, adcLo_lt _ _ _, adcLo_lt _ _ _, adcLo_lt _ _ _⟩ + +/-- The limbwise sum and the discarded top carry recompose the sum of the inputs. -/ +theorem addLimbs_toNat (a b : Limbs8) (ha : a.Bounded) (hb : b.Bounded) : + ∃ c : ℕ, c ≤ 1 ∧ (addLimbs a b).toNat + 2 ^ 256 * c = a.toNat + b.toNat := by + obtain ⟨ha0, ha1, ha2, ha3, ha4, ha5, ha6, ha7⟩ := ha + obtain ⟨hb0, hb1, hb2, hb3, hb4, hb5, hb6, hb7⟩ := hb + obtain ⟨s0, c0, ds0, dc0, e0, gc0, hs0⟩ := + adc_spec a.l0 b.l0 (0 : UInt64) ha0 hb0 (by decide) + obtain ⟨s1, c1, ds1, dc1, e1, gc1, hs1⟩ := + adc_spec a.l1 b.l1 _ ha1 hb1 (dc0 ▸ gc0) + obtain ⟨s2, c2, ds2, dc2, e2, gc2, hs2⟩ := + adc_spec a.l2 b.l2 _ ha2 hb2 (dc1 ▸ gc1) + obtain ⟨s3, c3, ds3, dc3, e3, gc3, hs3⟩ := + adc_spec a.l3 b.l3 _ ha3 hb3 (dc2 ▸ gc2) + obtain ⟨s4, c4, ds4, dc4, e4, gc4, hs4⟩ := + adc_spec a.l4 b.l4 _ ha4 hb4 (dc3 ▸ gc3) + obtain ⟨s5, c5, ds5, dc5, e5, gc5, hs5⟩ := + adc_spec a.l5 b.l5 _ ha5 hb5 (dc4 ▸ gc4) + obtain ⟨s6, c6, ds6, dc6, e6, gc6, hs6⟩ := + adc_spec a.l6 b.l6 _ ha6 hb6 (dc5 ▸ gc5) + obtain ⟨s7, c7, ds7, dc7, e7, gc7, hs7⟩ := + adc_spec a.l7 b.l7 _ ha7 hb7 (dc6 ▸ gc6) + simp only [dc0, dc1, dc2, dc3, dc4, dc5, dc6, toNat_zero] at e0 e1 e2 e3 e4 e5 e6 e7 + refine ⟨c7, gc7, ?_⟩ + simp only [addLimbs, Limbs8.toNat, ds0, ds1, ds2, ds3, ds4, ds5, ds6, ds7] + have h := carry_chain_sum e0 e1 e2 e3 e4 e5 e6 e7 + simp only [Nat.add_zero] at h + exact h + +theorem subLimbs_bounded (a b : Limbs8) : (subLimbs a b).Bounded := by + simp only [subLimbs, Limbs8.Bounded] + exact ⟨sbbLo_lt _ _ _, sbbLo_lt _ _ _, sbbLo_lt _ _ _, sbbLo_lt _ _ _, + sbbLo_lt _ _ _, sbbLo_lt _ _ _, sbbLo_lt _ _ _, sbbLo_lt _ _ _⟩ + +/-- The limbwise difference and the final borrow recompose the difference of the inputs. -/ +theorem subLimbs_spec (a b : Limbs8) (ha : a.Bounded) (hb : b.Bounded) : + (subBorrow a b).toNat ≤ 1 ∧ + (subLimbs a b).toNat + b.toNat = a.toNat + 2 ^ 256 * (subBorrow a b).toNat := by + obtain ⟨ha0, ha1, ha2, ha3, ha4, ha5, ha6, ha7⟩ := ha + obtain ⟨hb0, hb1, hb2, hb3, hb4, hb5, hb6, hb7⟩ := hb + obtain ⟨d0, c0, dd0, dc0, e0, gc0, hd0⟩ := + sbb_spec a.l0 b.l0 (0 : UInt64) ha0 hb0 (by decide) + obtain ⟨d1, c1, dd1, dc1, e1, gc1, hd1⟩ := + sbb_spec a.l1 b.l1 _ ha1 hb1 (dc0 ▸ gc0) + obtain ⟨d2, c2, dd2, dc2, e2, gc2, hd2⟩ := + sbb_spec a.l2 b.l2 _ ha2 hb2 (dc1 ▸ gc1) + obtain ⟨d3, c3, dd3, dc3, e3, gc3, hd3⟩ := + sbb_spec a.l3 b.l3 _ ha3 hb3 (dc2 ▸ gc2) + obtain ⟨d4, c4, dd4, dc4, e4, gc4, hd4⟩ := + sbb_spec a.l4 b.l4 _ ha4 hb4 (dc3 ▸ gc3) + obtain ⟨d5, c5, dd5, dc5, e5, gc5, hd5⟩ := + sbb_spec a.l5 b.l5 _ ha5 hb5 (dc4 ▸ gc4) + obtain ⟨d6, c6, dd6, dc6, e6, gc6, hd6⟩ := + sbb_spec a.l6 b.l6 _ ha6 hb6 (dc5 ▸ gc5) + obtain ⟨d7, c7, dd7, dc7, e7, gc7, hd7⟩ := + sbb_spec a.l7 b.l7 _ ha7 hb7 (dc6 ▸ gc6) + simp only [dc0, dc1, dc2, dc3, dc4, dc5, dc6, toNat_zero] at e0 e1 e2 e3 e4 e5 e6 e7 + refine ⟨?_, ?_⟩ + · simp only [subBorrow, dc7] + exact gc7 + · simp only [subLimbs, subBorrow, Limbs8.toNat, dd0, dd1, dd2, dd3, dd4, dd5, dd6, dd7, + dc7] + have h := borrow_chain_sum e0 e1 e2 e3 e4 e5 e6 e7 + simp only [Nat.add_zero] at h + exact h + +/-! ## Conditional subtraction -/ + +theorem condSub_bounded (q t : Limbs8) (ht : t.Bounded) : (condSub q t).Bounded := by + simp only [condSub] + split + · exact subLimbs_bounded t q + · exact ht + +/-- `condSub` subtracts the modulus exactly when the input is at least the modulus. -/ +theorem condSub_toNat (q t : Limbs8) (hq : q.Bounded) (ht : t.Bounded) : + (condSub q t).toNat = if t.toNat < q.toNat then t.toNat else t.toNat - q.toNat := by + obtain ⟨hbo, hchain⟩ := subLimbs_spec t q ht hq + have hD := Limbs8.toNat_lt (subLimbs_bounded t q) + simp only [condSub, beq_iff_eq, ← UInt64.toNat_inj, toNat_zero] + split + case isTrue h => + rw [h, Nat.mul_zero, Nat.add_zero] at hchain + exact cond_of_borrow_zero hchain hD + case isFalse h => + rw [show (subBorrow t q).toNat = 1 by omega, Nat.mul_one] at hchain + exact cond_of_borrow_one hchain hD + +/-- `condSub` returns a canonical representative for inputs below `2 * q`. -/ +theorem condSub_lt (q t : Limbs8) (hq : q.Bounded) (ht : t.Bounded) + (h : t.toNat < 2 * q.toNat) : (condSub q t).toNat < q.toNat := by + rw [condSub_toNat q t hq ht] + split <;> omega + +/-! ## Field operations -/ + +theorem add_bounded (q a b : Limbs8) : (add q a b).Bounded := + condSub_bounded _ _ (addLimbs_bounded a b) + +theorem sub_bounded (q a b : Limbs8) : (sub q a b).Bounded := by + simp only [sub] + split + · exact subLimbs_bounded a b + · exact addLimbs_bounded _ _ + +theorem neg_bounded (q a : Limbs8) : (neg q a).Bounded := sub_bounded _ _ _ + +/-- Modular addition is correct for canonical inputs of a modulus below `2 ^ 255`. -/ +theorem add_toNat (q a b : Limbs8) (hq : q.Bounded) (ha : a.Bounded) (hb : b.Bounded) + (hq2 : 2 * q.toNat < 2 ^ 256) (haq : a.toNat < q.toNat) (hbq : b.toNat < q.toNat) : + (add q a b).toNat = (a.toNat + b.toNat) % q.toNat := by + obtain ⟨c, hc, hchain⟩ := addLimbs_toNat a b ha hb + obtain ⟨hsum, hlt⟩ := + carry_top_zero hchain (Limbs8.toNat_lt (addLimbs_bounded a b)) haq hbq hq2 + simp only [add] + rw [condSub_toNat q _ hq (addLimbs_bounded a b), hsum] + rw [hsum] at hlt + exact cond_eq_mod hlt + +theorem add_lt (q a b : Limbs8) (hq : q.Bounded) (ha : a.Bounded) (hb : b.Bounded) + (hq2 : 2 * q.toNat < 2 ^ 256) (haq : a.toNat < q.toNat) (hbq : b.toNat < q.toNat) : + (add q a b).toNat < q.toNat := by + rw [add_toNat q a b hq ha hb hq2 haq hbq] + exact Nat.mod_lt _ (by omega) + +/-- Modular subtraction is correct for canonical inputs of a modulus below `2 ^ 255`. -/ +theorem sub_toNat (q a b : Limbs8) (hq : q.Bounded) (ha : a.Bounded) (hb : b.Bounded) + (hq2 : 2 * q.toNat < 2 ^ 256) (haq : a.toNat < q.toNat) (hbq : b.toNat < q.toNat) : + (sub q a b).toNat = (a.toNat + (q.toNat - b.toNat)) % q.toNat := by + obtain ⟨hbo, hchain⟩ := subLimbs_spec a b ha hb + have hD := Limbs8.toNat_lt (subLimbs_bounded a b) + simp only [sub, beq_iff_eq, ← UInt64.toNat_inj, toNat_zero] + split + case isTrue h => + rw [h] at hchain + exact sub_of_borrow_zero hchain hbq haq + case isFalse h => + rw [show (subBorrow a b).toNat = 1 by omega] at hchain + obtain ⟨c, _, hchainE⟩ := addLimbs_toNat _ q (subLimbs_bounded a b) hq + exact sub_of_borrow_one hchain hchainE hD + (Limbs8.toNat_lt (addLimbs_bounded _ _)) hbq haq hq2 + +theorem sub_lt (q a b : Limbs8) (hq : q.Bounded) (ha : a.Bounded) (hb : b.Bounded) + (hq2 : 2 * q.toNat < 2 ^ 256) (haq : a.toNat < q.toNat) (hbq : b.toNat < q.toNat) : + (sub q a b).toNat < q.toNat := by + rw [sub_toNat q a b hq ha hb hq2 haq hbq] + exact Nat.mod_lt _ (by omega) + +/-- Modular negation is correct for canonical inputs of a modulus below `2 ^ 255`. -/ +theorem neg_toNat (q a : Limbs8) (hq : q.Bounded) (ha : a.Bounded) + (hq2 : 2 * q.toNat < 2 ^ 256) (haq : a.toNat < q.toNat) : + (neg q a).toNat = (q.toNat - a.toNat) % q.toNat := by + simp only [neg] + rw [sub_toNat q Limbs8.zero a hq Limbs8.zero_bounded ha hq2 + (by rw [Limbs8.zero_toNat]; omega) haq, Limbs8.zero_toNat, Nat.zero_add] + +theorem neg_lt (q a : Limbs8) (hq : q.Bounded) (ha : a.Bounded) + (hq2 : 2 * q.toNat < 2 ^ 256) (haq : a.toNat < q.toNat) : (neg q a).toNat < q.toNat := + sub_lt q Limbs8.zero a hq Limbs8.zero_bounded ha hq2 + (by rw [Limbs8.zero_toNat]; omega) haq + +/-! ## CIOS multiplication + +The definitions of the multiplication round and of the full product. Their correctness — +the round invariant and the `2 * q` bound on the accumulator — is established in the sibling +module. -/ + +namespace State9 + +end State9 + +end Native64x8 + +end Montgomery diff --git a/CompElliptic/Vendor/CompPoly/Montgomery/Native64x8Defs.lean b/CompElliptic/Vendor/CompPoly/Montgomery/Native64x8Defs.lean new file mode 100644 index 0000000..3968f14 --- /dev/null +++ b/CompElliptic/Vendor/CompPoly/Montgomery/Native64x8Defs.lean @@ -0,0 +1,350 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Gregor Mitscha-Baude +-/ + +/-! +# Eight-limb Montgomery arithmetic: runtime definitions (core-only) + +Vendored from CompPoly branch `fast_multilimb_fields`; see `README.md` in this directory. + +This module contains **only** the runtime definitions of +`CompPoly/Fields/Montgomery/Native64x8.lean`, moved here verbatim, plus the Pasta constants +and monomorphic entry points. It deliberately imports nothing beyond Lean core, because it is +the single module in the `FastFieldNative` precompiled lane: `precompileModules` native-compiles +the whole import closure, and a mathlib closure OOMs the build machine. + +All correctness statements about these definitions live in the sibling modules +`CompElliptic.Vendor.CompPoly.Montgomery.Native64x8` (raw operations) and +`CompElliptic.Vendor.CompPoly.Montgomery.Native64x8Mul` (CIOS multiplication), which import this one. +-/ + +namespace Montgomery + +namespace Native64x8 + +/-- Mask selecting the low 32 bits of a `UInt64`. -/ +@[inline] def mask : UInt64 := 0xffffffff + +/-- Low limb of add-with-carry: `(x + y + c) mod 2 ^ 32`. -/ +@[inline] def adcLo (x y c : UInt64) : UInt64 := (x + y + c) &&& mask + +/-- Carry-out of add-with-carry: `(x + y + c) / 2 ^ 32`. -/ +@[inline] def adcCo (x y c : UInt64) : UInt64 := (x + y + c) >>> 32 + +/-- Low limb of subtract-with-borrow: `(x - y - b) mod 2 ^ 32`. -/ +@[inline] def sbbLo (x y b : UInt64) : UInt64 := (x - y - b) &&& mask + +/-- Borrow-out of subtract-with-borrow, read off the sign bit of the 64-bit difference. -/ +@[inline] def sbbBo (x y b : UInt64) : UInt64 := (x - y - b) >>> 63 + +/-- Low limb of multiply-accumulate: `(t + x * y + c) mod 2 ^ 32`. -/ +@[inline] def macLo (t x y c : UInt64) : UInt64 := (t + x * y + c) &&& mask + +/-- High word of multiply-accumulate: `(t + x * y + c) / 2 ^ 32`. -/ +@[inline] def macHi (t x y c : UInt64) : UInt64 := (t + x * y + c) >>> 32 + +/-- The Montgomery multiplier of a limb: `(s * negInv) mod 2 ^ 32`. -/ +@[inline] def montM (s negInv : UInt64) : UInt64 := ((s &&& mask) * negInv) &&& mask + +/-- A 256-bit value as eight little-endian 32-bit limbs, each stored in a `UInt64`. -/ +structure Limbs8 where + /-- Limb of weight `2 ^ 0`. -/ + l0 : UInt64 + /-- Limb of weight `2 ^ 32`. -/ + l1 : UInt64 + /-- Limb of weight `2 ^ 64`. -/ + l2 : UInt64 + /-- Limb of weight `2 ^ 96`. -/ + l3 : UInt64 + /-- Limb of weight `2 ^ 128`. -/ + l4 : UInt64 + /-- Limb of weight `2 ^ 160`. -/ + l5 : UInt64 + /-- Limb of weight `2 ^ 192`. -/ + l6 : UInt64 + /-- Limb of weight `2 ^ 224`. -/ + l7 : UInt64 +deriving DecidableEq, Repr, Inhabited + +namespace Limbs8 + +/-- The zero value. -/ +def zero : Limbs8 := ⟨0, 0, 0, 0, 0, 0, 0, 0⟩ + +/-- The value one. -/ +def one : Limbs8 := ⟨1, 0, 0, 0, 0, 0, 0, 0⟩ + +/-- Split a natural number into eight 32-bit limbs, discarding bits above `2 ^ 256`. -/ +@[inline] def ofNat (n : Nat) : Limbs8 := + ⟨UInt64.ofNat (n &&& 0xffffffff), UInt64.ofNat (n >>> 32 &&& 0xffffffff), + UInt64.ofNat (n >>> 64 &&& 0xffffffff), UInt64.ofNat (n >>> 96 &&& 0xffffffff), + UInt64.ofNat (n >>> 128 &&& 0xffffffff), UInt64.ofNat (n >>> 160 &&& 0xffffffff), + UInt64.ofNat (n >>> 192 &&& 0xffffffff), UInt64.ofNat (n >>> 224 &&& 0xffffffff)⟩ + +/-- The natural number represented by the limbs: `∑ lᵢ * 2 ^ (32 * i)`. -/ +def toNat (x : Limbs8) : Nat := + x.l0.toNat + 2 ^ 32 * x.l1.toNat + 2 ^ 64 * x.l2.toNat + 2 ^ 96 * x.l3.toNat + + 2 ^ 128 * x.l4.toNat + 2 ^ 160 * x.l5.toNat + 2 ^ 192 * x.l6.toNat + 2 ^ 224 * x.l7.toNat + +/-- Every limb holds at most 32 significant bits. -/ +def Bounded (x : Limbs8) : Prop := + x.l0.toNat < 2 ^ 32 ∧ x.l1.toNat < 2 ^ 32 ∧ x.l2.toNat < 2 ^ 32 ∧ x.l3.toNat < 2 ^ 32 ∧ + x.l4.toNat < 2 ^ 32 ∧ x.l5.toNat < 2 ^ 32 ∧ x.l6.toNat < 2 ^ 32 ∧ x.l7.toNat < 2 ^ 32 + +instance (x : Limbs8) : Decidable x.Bounded := by + unfold Bounded + infer_instance + +end Limbs8 + +/-- Limbwise add-with-carry, discarding the carry out of the top limb. -/ +@[inline] def addLimbs (a b : Limbs8) : Limbs8 := + let c0 := adcCo a.l0 b.l0 0 + let c1 := adcCo a.l1 b.l1 c0 + let c2 := adcCo a.l2 b.l2 c1 + let c3 := adcCo a.l3 b.l3 c2 + let c4 := adcCo a.l4 b.l4 c3 + let c5 := adcCo a.l5 b.l5 c4 + let c6 := adcCo a.l6 b.l6 c5 + ⟨adcLo a.l0 b.l0 0, adcLo a.l1 b.l1 c0, adcLo a.l2 b.l2 c1, adcLo a.l3 b.l3 c2, + adcLo a.l4 b.l4 c3, adcLo a.l5 b.l5 c4, adcLo a.l6 b.l6 c5, adcLo a.l7 b.l7 c6⟩ + +/-- Limbwise subtract-with-borrow. -/ +@[inline] def subLimbs (a b : Limbs8) : Limbs8 := + let b0 := sbbBo a.l0 b.l0 0 + let b1 := sbbBo a.l1 b.l1 b0 + let b2 := sbbBo a.l2 b.l2 b1 + let b3 := sbbBo a.l3 b.l3 b2 + let b4 := sbbBo a.l4 b.l4 b3 + let b5 := sbbBo a.l5 b.l5 b4 + let b6 := sbbBo a.l6 b.l6 b5 + ⟨sbbLo a.l0 b.l0 0, sbbLo a.l1 b.l1 b0, sbbLo a.l2 b.l2 b1, sbbLo a.l3 b.l3 b2, + sbbLo a.l4 b.l4 b3, sbbLo a.l5 b.l5 b4, sbbLo a.l6 b.l6 b5, sbbLo a.l7 b.l7 b6⟩ + +/-- Borrow out of the top limb of `subLimbs`. -/ +@[inline] def subBorrow (a b : Limbs8) : UInt64 := + let b0 := sbbBo a.l0 b.l0 0 + let b1 := sbbBo a.l1 b.l1 b0 + let b2 := sbbBo a.l2 b.l2 b1 + let b3 := sbbBo a.l3 b.l3 b2 + let b4 := sbbBo a.l4 b.l4 b3 + let b5 := sbbBo a.l5 b.l5 b4 + let b6 := sbbBo a.l6 b.l6 b5 + sbbBo a.l7 b.l7 b6 + +/-- Subtract the modulus once if the value is at least the modulus. The borrow chain +decides the branch, so no comparison is needed. -/ +@[inline] def condSub (q t : Limbs8) : Limbs8 := + if subBorrow t q == 0 then subLimbs t q else t + +/-- Modular addition. -/ +@[inline] def add (q a b : Limbs8) : Limbs8 := condSub q (addLimbs a b) + +/-- Modular subtraction: on a borrow, the modulus is added back. -/ +@[inline] def sub (q a b : Limbs8) : Limbs8 := + let d := subLimbs a b + if subBorrow a b == 0 then d else addLimbs d q + +/-- Modular negation. -/ +@[inline] def neg (q a : Limbs8) : Limbs8 := sub q Limbs8.zero a + +/-- The CIOS accumulator: eight limbs plus one head limb. -/ +structure State9 where + /-- Limb of weight `2 ^ 0`. -/ + t0 : UInt64 + /-- Limb of weight `2 ^ 32`. -/ + t1 : UInt64 + /-- Limb of weight `2 ^ 64`. -/ + t2 : UInt64 + /-- Limb of weight `2 ^ 96`. -/ + t3 : UInt64 + /-- Limb of weight `2 ^ 128`. -/ + t4 : UInt64 + /-- Limb of weight `2 ^ 160`. -/ + t5 : UInt64 + /-- Limb of weight `2 ^ 192`. -/ + t6 : UInt64 + /-- Limb of weight `2 ^ 224`. -/ + t7 : UInt64 + /-- Head limb of weight `2 ^ 256`. -/ + t8 : UInt64 +deriving DecidableEq, Repr, Inhabited + +namespace State9 + +/-- The zero accumulator. -/ +@[inline] def zero : State9 := ⟨0, 0, 0, 0, 0, 0, 0, 0, 0⟩ + +/-- The eight low limbs of the accumulator. -/ +@[inline] def toLimbs8 (t : State9) : Limbs8 := ⟨t.t0, t.t1, t.t2, t.t3, t.t4, t.t5, t.t6, t.t7⟩ + +/-- The natural number represented by the accumulator. -/ +def toNat (t : State9) : Nat := t.toLimbs8.toNat + 2 ^ 256 * t.t8.toNat + +/-- Every limb of the accumulator holds at most 32 significant bits. -/ +def Bounded (t : State9) : Prop := t.toLimbs8.Bounded ∧ t.t8.toNat < 2 ^ 32 + +end State9 + +/-- The multiply half of a CIOS round: accumulate `a * bi` into the accumulator. The carry +out of the top limb is kept in the head limb, so no information is lost. -/ +@[inline] def mulAccum (a : Limbs8) (bi : UInt64) (t : State9) : State9 := + let k0 := macHi t.t0 a.l0 bi 0 + let k1 := macHi t.t1 a.l1 bi k0 + let k2 := macHi t.t2 a.l2 bi k1 + let k3 := macHi t.t3 a.l3 bi k2 + let k4 := macHi t.t4 a.l4 bi k3 + let k5 := macHi t.t5 a.l5 bi k4 + let k6 := macHi t.t6 a.l6 bi k5 + let k7 := macHi t.t7 a.l7 bi k6 + ⟨macLo t.t0 a.l0 bi 0, macLo t.t1 a.l1 bi k0, macLo t.t2 a.l2 bi k1, + macLo t.t3 a.l3 bi k2, macLo t.t4 a.l4 bi k3, macLo t.t5 a.l5 bi k4, + macLo t.t6 a.l6 bi k5, macLo t.t7 a.l7 bi k6, t.t8 + k7⟩ + +/-- The reduce half of a CIOS round: add the multiple `montM s.t0 negInv` of the modulus that +cancels the low limb, then drop that limb. `negInv` has to be `-q⁻¹ mod 2 ^ 32`. -/ +@[inline] def mulReduce (q : Limbs8) (negInv : UInt64) (s : State9) : State9 := + let m := montM s.t0 negInv + let u0 := macHi s.t0 m q.l0 0 + let u1 := macHi s.t1 m q.l1 u0 + let u2 := macHi s.t2 m q.l2 u1 + let u3 := macHi s.t3 m q.l3 u2 + let u4 := macHi s.t4 m q.l4 u3 + let u5 := macHi s.t5 m q.l5 u4 + let u6 := macHi s.t6 m q.l6 u5 + let u7 := macHi s.t7 m q.l7 u6 + ⟨macLo s.t1 m q.l1 u0, macLo s.t2 m q.l2 u1, macLo s.t3 m q.l3 u2, + macLo s.t4 m q.l4 u3, macLo s.t5 m q.l5 u4, macLo s.t6 m q.l6 u5, + macLo s.t7 m q.l7 u6, adcLo s.t8 u7 0, adcCo s.t8 u7 0⟩ + +/-- One CIOS outer round: accumulate `a * bi` into the accumulator, then reduce away one +limb. -/ +@[inline] def mulRound (q : Limbs8) (negInv : UInt64) (a : Limbs8) (bi : UInt64) + (t : State9) : State9 := + mulReduce q negInv (mulAccum a bi t) + +/-- CIOS Montgomery multiplication: eight rounds followed by one conditional +subtraction. -/ +@[inline] def mul (q : Limbs8) (negInv : UInt64) (a b : Limbs8) : Limbs8 := + let t := mulRound q negInv a b.l0 State9.zero + let t := mulRound q negInv a b.l1 t + let t := mulRound q negInv a b.l2 t + let t := mulRound q negInv a b.l3 t + let t := mulRound q negInv a b.l4 t + let t := mulRound q negInv a b.l5 t + let t := mulRound q negInv a b.l6 t + let t := mulRound q negInv a b.l7 t + condSub q t.toLimbs8 + +/-- Montgomery squaring. -/ +@[inline] def square (q : Limbs8) (negInv : UInt64) (a : Limbs8) : Limbs8 := + mul q negInv a a + +/-! ## Pasta field constants and specialized entry points + +The two Pasta base fields, as plain data plus monomorphic wrappers. Everything in this +section is core-only so that it can live in the `FastFieldNative` precompiled lane; the +`Mont64x8Field` instances that carry the correctness side conditions live in +`CompElliptic.Vendor.CompPoly.Montgomery.Pasta`. -/ + +namespace VestaFq + +/-- The Vesta base field modulus in eight 32-bit limbs. -/ +def modulusLimbs : Limbs8 := ⟨0x1, 0x8c46eb21, 0x994a8dd, 0x224698fc, 0x0, 0x0, 0x0, 0x40000000⟩ + +/-- `2 ^ 256 mod q`, the Montgomery representation of one. -/ +def rModModulus : Limbs8 := + ⟨0xfffffffd, 0x5b2b3e9c, 0xe3420567, 0x992c350b, 0xffffffff, 0xffffffff, 0xffffffff, + 0x3fffffff⟩ + +/-- `(2 ^ 256) ^ 2 mod q`, used to enter Montgomery form. -/ +def r2ModModulus : Limbs8 := + ⟨0xf, 0xfc9678ff, 0x891a16e3, 0x67bb433d, 0x4ccf590, 0x7fae2310, 0x7ccfdaa9, 0x96d41af⟩ + +/-- `-q⁻¹ mod 2 ^ 32`. -/ +def negInv : UInt64 := 0xffffffff + +/-- Montgomery-form zero. -/ +def zero : Limbs8 := Limbs8.zero + +/-- Montgomery-form one. -/ +def one : Limbs8 := rModModulus + +/-- Modular addition in the Vesta base field. -/ +@[inline] def add (a b : Limbs8) : Limbs8 := Native64x8.add modulusLimbs a b + +/-- Modular subtraction in the Vesta base field. -/ +@[inline] def sub (a b : Limbs8) : Limbs8 := Native64x8.sub modulusLimbs a b + +/-- Modular negation in the Vesta base field. -/ +@[inline] def neg (a : Limbs8) : Limbs8 := Native64x8.neg modulusLimbs a + +/-- Montgomery multiplication in the Vesta base field. -/ +@[inline] def mul (a b : Limbs8) : Limbs8 := Native64x8.mul modulusLimbs negInv a b + +/-- Montgomery squaring in the Vesta base field. -/ +@[inline] def square (a : Limbs8) : Limbs8 := Native64x8.mul modulusLimbs negInv a a + +/-- Enter Montgomery form from a canonical natural number below `q`. -/ +@[inline] def ofNat (n : Nat) : Limbs8 := + Native64x8.mul modulusLimbs negInv (Limbs8.ofNat n) r2ModModulus + +/-- Leave Montgomery form: the canonical limb representative. -/ +@[inline] def toLimbs8 (a : Limbs8) : Limbs8 := + Native64x8.mul modulusLimbs negInv a Limbs8.one + +end VestaFq + +namespace PallasFq + +/-- The Pallas base field modulus in eight 32-bit limbs. -/ +def modulusLimbs : Limbs8 := ⟨0x1, 0x992d30ed, 0x94cf91b, 0x224698fc, 0x0, 0x0, 0x0, 0x40000000⟩ + +/-- `2 ^ 256 mod p`, the Montgomery representation of one. -/ +def rModModulus : Limbs8 := + ⟨0xfffffffd, 0x34786d38, 0xe41914ad, 0x992c350b, 0xffffffff, 0xffffffff, 0xffffffff, + 0x3fffffff⟩ + +/-- `(2 ^ 256) ^ 2 mod p`, used to enter Montgomery form. -/ +def r2ModModulus : Limbs8 := + ⟨0xf, 0x8c78ecb3, 0x8b0de0e7, 0xd7d30dbd, 0xc3c95d18, 0x7797a99b, 0x7b9cb714, 0x96d41af⟩ + +/-- `-p⁻¹ mod 2 ^ 32`. -/ +def negInv : UInt64 := 0xffffffff + +/-- Montgomery-form zero. -/ +def zero : Limbs8 := Limbs8.zero + +/-- Montgomery-form one. -/ +def one : Limbs8 := rModModulus + +/-- Modular addition in the Pallas base field. -/ +@[inline] def add (a b : Limbs8) : Limbs8 := Native64x8.add modulusLimbs a b + +/-- Modular subtraction in the Pallas base field. -/ +@[inline] def sub (a b : Limbs8) : Limbs8 := Native64x8.sub modulusLimbs a b + +/-- Modular negation in the Pallas base field. -/ +@[inline] def neg (a : Limbs8) : Limbs8 := Native64x8.neg modulusLimbs a + +/-- Montgomery multiplication in the Pallas base field. -/ +@[inline] def mul (a b : Limbs8) : Limbs8 := Native64x8.mul modulusLimbs negInv a b + +/-- Montgomery squaring in the Pallas base field. -/ +@[inline] def square (a : Limbs8) : Limbs8 := Native64x8.mul modulusLimbs negInv a a + +/-- Enter Montgomery form from a canonical natural number below `p`. -/ +@[inline] def ofNat (n : Nat) : Limbs8 := + Native64x8.mul modulusLimbs negInv (Limbs8.ofNat n) r2ModModulus + +/-- Leave Montgomery form: the canonical limb representative. -/ +@[inline] def toLimbs8 (a : Limbs8) : Limbs8 := + Native64x8.mul modulusLimbs negInv a Limbs8.one + +end PallasFq + +end Native64x8 + +end Montgomery diff --git a/CompElliptic/Vendor/CompPoly/Montgomery/Native64x8Field.lean b/CompElliptic/Vendor/CompPoly/Montgomery/Native64x8Field.lean new file mode 100644 index 0000000..d5164d7 --- /dev/null +++ b/CompElliptic/Vendor/CompPoly/Montgomery/Native64x8Field.lean @@ -0,0 +1,551 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Gregor Mitscha-Baude +-/ + +import CompPoly.Fields.Basic +import CompElliptic.Vendor.CompPoly.Montgomery.Native64x8Mul +import Mathlib.Algebra.Field.TransferInstance +import Mathlib.FieldTheory.Finite.Basic + +/-! +# Fast eight-limb Montgomery fields + +The bounded carrier, conversions, arithmetic, and field instances built on the raw eight-limb +Montgomery operations of `Montgomery/Native64x8`. This is the multi-limb analogue of +`Montgomery/Native32Field`, for prime moduli below `2 ^ 255`. + +A carrier element stores the Montgomery residue `x * 2 ^ 256 mod q` as eight 32-bit limbs; +`toField` divides that residue by `2 ^ 256` again and lands in `ZMod modulus`. All arithmetic +is transported along `toField`, which is a ring isomorphism onto `ZMod modulus`. + +## Main results + +* `FastField` — the carrier, `{ x : Limbs8 // x.Bounded ∧ x.toNat < modulus }` +* `toField_add`, `toField_mul`, … — `@[simp]` equivalences with the canonical field +* `Mont64x8Field` — the per-field constants class +* `ringEquiv`, `instField` — the ring isomorphism and the transferred `Field` instance +-/ + +namespace Montgomery +namespace Native64x8 + +/-! ## Per-field constants -/ + +/-- Per-field data for a fast eight-limb Montgomery field with radix `R = 2 ^ 256`. All +side conditions are concrete numeral facts, discharged by `decide` at instantiation. -/ +class Mont64x8Field (modulus : ℕ) where + /-- `modulus` is prime. -/ + prime : modulus.Prime + /-- The modulus in eight 32-bit limbs. -/ + modulusLimbs : Limbs8 + /-- `2 ^ 256 mod modulus`, the Montgomery representation of one. -/ + rModModulus : Limbs8 + /-- `(2 ^ 256) ^ 2 mod modulus`, used to enter Montgomery form. -/ + r2ModModulus : Limbs8 + /-- `-modulus⁻¹ mod 2 ^ 32`, used by Montgomery reduction. -/ + montgomeryNegInv : UInt64 + modulusLimbs_bounded : modulusLimbs.Bounded := by decide + modulusLimbs_toNat : modulusLimbs.toNat = modulus := by decide + two_lt_modulus : 2 < modulus := by decide + two_mul_modulus_lt : 2 * modulus < 2 ^ 256 := by decide + rModModulus_bounded : rModModulus.Bounded := by decide + rModModulus_toNat : rModModulus.toNat = 2 ^ 256 % modulus := by decide + r2ModModulus_bounded : r2ModModulus.Bounded := by decide + r2ModModulus_toNat : r2ModModulus.toNat = (2 ^ 256) ^ 2 % modulus := by decide + montgomeryNegInv_lt : montgomeryNegInv.toNat < 2 ^ 32 := by decide + montgomeryNegInv_mul_modulus_mod_two_pow_32 : + montgomeryNegInv.toNat * modulus % 2 ^ 32 = 2 ^ 32 - 1 := by decide + +/-! ## Modulus facts -/ + +namespace Mont64x8Field + +variable {modulus : ℕ} [P : Mont64x8Field modulus] + +instance : Fact (Nat.Prime modulus) := ⟨P.prime⟩ + +theorem modulus_pos : 0 < modulus := Nat.zero_lt_of_lt P.two_lt_modulus + +theorem modulus_lt : modulus < 2 ^ 256 := by + have := P.two_mul_modulus_lt + omega + +theorem q_toNat : P.modulusLimbs.toNat = modulus := P.modulusLimbs_toNat + +theorem two_mul_q_lt : 2 * P.modulusLimbs.toNat < 2 ^ 256 := by + rw [q_toNat] + exact P.two_mul_modulus_lt + +theorem negInv_mul_q : P.montgomeryNegInv.toNat * P.modulusLimbs.toNat % 2 ^ 32 = 2 ^ 32 - 1 := by + rw [q_toNat] + exact P.montgomeryNegInv_mul_modulus_mod_two_pow_32 + +instance : NeZero modulus := ⟨modulus_pos.ne'⟩ + +theorem two_ne_zero : (2 : ZMod modulus) ≠ 0 := by + intro h + have hdvd : modulus ∣ 2 := (ZMod.natCast_eq_zero_iff 2 modulus).mp h + exact (Nat.not_le_of_gt P.two_lt_modulus) (Nat.le_of_dvd (by decide) hdvd) + +theorem r_ne_zero : ((2 ^ 256 : ℕ) : ZMod modulus) ≠ 0 := by + rw [Nat.cast_pow, Nat.cast_ofNat] + exact pow_ne_zero _ two_ne_zero + +end Mont64x8Field + +/-! ## The carrier -/ + +/-- The fast carrier for a prime modulus: eight 32-bit limbs holding a value below `modulus`, +interpreted as a Montgomery residue. At runtime this erases to `Limbs8`. -/ +def FastField (modulus : ℕ) [Mont64x8Field modulus] : Type := + { x : Limbs8 // x.Bounded ∧ x.toNat < modulus } + +section + +variable {modulus : ℕ} [P : Mont64x8Field modulus] + +instance : DecidableEq (FastField modulus) := + inferInstanceAs (DecidableEq { x : Limbs8 // x.Bounded ∧ x.toNat < modulus }) + +namespace FastField + +theorem val_bounded (x : FastField modulus) : x.val.Bounded := x.property.1 + +theorem val_lt (x : FastField modulus) : x.val.toNat < P.modulusLimbs.toNat := by + rw [Mont64x8Field.q_toNat] + exact x.property.2 + +/-! ### Arithmetic -/ + +/-- The zero element. -/ +def zero (modulus : ℕ) [P : Mont64x8Field modulus] : FastField modulus := + ⟨Limbs8.zero, Limbs8.zero_bounded, by + rw [Limbs8.zero_toNat] + exact Mont64x8Field.modulus_pos⟩ + +/-- The one element, the Montgomery residue `2 ^ 256 mod modulus`. -/ +def one (modulus : ℕ) [P : Mont64x8Field modulus] : FastField modulus := + ⟨P.rModModulus, P.rModModulus_bounded, by + rw [P.rModModulus_toNat] + exact Nat.mod_lt _ Mont64x8Field.modulus_pos⟩ + +/-- Fast modular addition in Montgomery form. -/ +@[inline] def add (x y : FastField modulus) : FastField modulus := + ⟨Native64x8.add P.modulusLimbs x.val y.val, add_bounded _ _ _, + by + have h := add_lt _ _ _ P.modulusLimbs_bounded x.val_bounded y.val_bounded + Mont64x8Field.two_mul_q_lt x.val_lt y.val_lt + rwa [Mont64x8Field.q_toNat] at h⟩ + +/-- Fast modular subtraction in Montgomery form. -/ +@[inline] def sub (x y : FastField modulus) : FastField modulus := + ⟨Native64x8.sub P.modulusLimbs x.val y.val, sub_bounded _ _ _, + by + have h := sub_lt _ _ _ P.modulusLimbs_bounded x.val_bounded y.val_bounded + Mont64x8Field.two_mul_q_lt x.val_lt y.val_lt + rwa [Mont64x8Field.q_toNat] at h⟩ + +/-- Fast modular negation in Montgomery form. -/ +@[inline] def neg (x : FastField modulus) : FastField modulus := + ⟨Native64x8.neg P.modulusLimbs x.val, neg_bounded _ _, + by + have h := neg_lt _ _ P.modulusLimbs_bounded x.val_bounded Mont64x8Field.two_mul_q_lt + x.val_lt + rwa [Mont64x8Field.q_toNat] at h⟩ + +/-- Fast Montgomery multiplication. -/ +@[inline] def mul (x y : FastField modulus) : FastField modulus := + ⟨Native64x8.mul P.modulusLimbs P.montgomeryNegInv x.val y.val, + (mul_spec _ _ _ _ P.modulusLimbs_bounded x.val_bounded y.val_bounded + P.montgomeryNegInv_lt Mont64x8Field.negInv_mul_q x.val_lt + Mont64x8Field.two_mul_q_lt).1, + by + have h := (mul_spec _ _ _ _ P.modulusLimbs_bounded x.val_bounded y.val_bounded + P.montgomeryNegInv_lt Mont64x8Field.negInv_mul_q x.val_lt + Mont64x8Field.two_mul_q_lt).2.1 + rwa [Mont64x8Field.q_toNat] at h⟩ + +/-- Fast squaring. -/ +@[inline] def square (x : FastField modulus) : FastField modulus := mul x x + +/-- Exponentiation by repeated squaring. -/ +@[specialize] def pow (x : FastField modulus) (n : ℕ) : FastField modulus := + @npowBinRec (FastField modulus) ⟨one modulus⟩ ⟨mul⟩ n x + +/-- Inversion by Fermat's little theorem, `x⁻¹ = x ^ (modulus - 2)`. -/ +@[inline] def inv (x : FastField modulus) : FastField modulus := pow x (modulus - 2) + +/-- Division through inversion. -/ +@[inline] def div (x y : FastField modulus) : FastField modulus := mul x (inv y) + +/-! ### Conversions -/ + +/-- The canonical limb representative of a fast element: one Montgomery reduction. -/ +@[inline] def toLimbs8 (x : FastField modulus) : Limbs8 := + Native64x8.mul P.modulusLimbs P.montgomeryNegInv x.val Limbs8.one + +/-- The canonical natural representative of a fast element. -/ +@[inline] def toNat (x : FastField modulus) : ℕ := (toLimbs8 x).toNat + +/-- The canonical `ZMod` value of a fast element. -/ +@[inline] def toField (x : FastField modulus) : ZMod modulus := (toNat x : ZMod modulus) + +/-- Build a fast element from a canonical natural representative. -/ +@[inline] def ofCanonicalNat (n : ℕ) (h : n < modulus) : FastField modulus := + ⟨Native64x8.mul P.modulusLimbs P.montgomeryNegInv (Limbs8.ofNat n) P.r2ModModulus, + (mul_spec _ _ _ _ P.modulusLimbs_bounded (Limbs8.ofNat_bounded n) P.r2ModModulus_bounded + P.montgomeryNegInv_lt Mont64x8Field.negInv_mul_q + (by + rw [Limbs8.ofNat_toNat, Mont64x8Field.q_toNat, + Nat.mod_eq_of_lt (h.trans Mont64x8Field.modulus_lt)] + exact h) + Mont64x8Field.two_mul_q_lt).1, + by + have hlt := (mul_spec _ _ _ _ P.modulusLimbs_bounded (Limbs8.ofNat_bounded n) + P.r2ModModulus_bounded P.montgomeryNegInv_lt Mont64x8Field.negInv_mul_q + (by + rw [Limbs8.ofNat_toNat, Mont64x8Field.q_toNat, + Nat.mod_eq_of_lt (h.trans Mont64x8Field.modulus_lt)] + exact h) + Mont64x8Field.two_mul_q_lt).2.1 + rwa [Mont64x8Field.q_toNat] at hlt⟩ + +/-- Convert a natural number into fast Montgomery form. -/ +@[inline] def ofNat (modulus : ℕ) [P : Mont64x8Field modulus] (n : ℕ) : FastField modulus := + ofCanonicalNat (n % modulus) (Nat.mod_lt _ Mont64x8Field.modulus_pos) + +/-- Convert from the canonical field into fast Montgomery form. -/ +@[inline] def ofField (x : ZMod modulus) : FastField modulus := + ofCanonicalNat x.val (ZMod.val_lt x) + +/-- Convert an integer into fast Montgomery form. -/ +@[inline] private def ofInt (modulus : ℕ) [P : Mont64x8Field modulus] (n : Int) : + FastField modulus := + ofField (n : ZMod modulus) + +instance : Zero (FastField modulus) := ⟨FastField.zero modulus⟩ +instance : One (FastField modulus) := ⟨FastField.one modulus⟩ +instance : Add (FastField modulus) := ⟨FastField.add⟩ +instance : Neg (FastField modulus) := ⟨FastField.neg⟩ +instance : Sub (FastField modulus) := ⟨FastField.sub⟩ +instance : Mul (FastField modulus) := ⟨FastField.mul⟩ +instance : Inv (FastField modulus) := ⟨FastField.inv⟩ +instance : Div (FastField modulus) := ⟨FastField.div⟩ +instance : NatCast (FastField modulus) := ⟨FastField.ofNat modulus⟩ +instance : IntCast (FastField modulus) := ⟨FastField.ofInt modulus⟩ +instance : Pow (FastField modulus) ℕ := ⟨FastField.pow⟩ + +instance : SMul ℕ (FastField modulus) where + smul n x := FastField.ofNat modulus n * x + +instance : SMul Int (FastField modulus) where + smul n x := FastField.ofInt modulus n * x + +instance : Pow (FastField modulus) Int where + pow x n := + match n with + | Int.ofNat k => FastField.pow x k + | Int.negSucc k => FastField.pow (FastField.inv x) (k + 1) + +theorem zero_def : (0 : FastField modulus) = FastField.zero modulus := rfl +theorem one_def : (1 : FastField modulus) = FastField.one modulus := rfl +theorem add_def (x y : FastField modulus) : x + y = FastField.add x y := rfl +theorem neg_def (x : FastField modulus) : -x = FastField.neg x := rfl +theorem sub_def (x y : FastField modulus) : x - y = FastField.sub x y := rfl +theorem mul_def (x y : FastField modulus) : x * y = FastField.mul x y := rfl +theorem inv_def (x : FastField modulus) : x⁻¹ = FastField.inv x := rfl +theorem div_def (x y : FastField modulus) : x / y = x * y⁻¹ := rfl +theorem square_def (x : FastField modulus) : FastField.square x = x * x := rfl + +/-! ## Correctness -/ + +section Bridge + +instance : NNRatCast (FastField modulus) where + nnratCast q := FastField.ofField (q : ZMod modulus) + +instance : RatCast (FastField modulus) where + ratCast q := FastField.ofField (q : ZMod modulus) + +instance : SMul ℚ≥0 (FastField modulus) where + smul q x := FastField.ofField (q • FastField.toField x) + +instance : SMul ℚ (FastField modulus) where + smul q x := FastField.ofField (q • FastField.toField x) + +/-- The raw Montgomery product divides by `2 ^ 256` in the canonical field. -/ +private theorem mul_cast (x y : Limbs8) (hx : x.Bounded) (hy : y.Bounded) + (hxq : x.toNat < modulus) : + ((Native64x8.mul P.modulusLimbs P.montgomeryNegInv x y).toNat : ZMod modulus) = + (x.toNat : ZMod modulus) * (y.toNat : ZMod modulus) * + ((2 ^ 256 : ℕ) : ZMod modulus)⁻¹ := by + have hmod := (mul_spec _ _ _ _ P.modulusLimbs_bounded hx hy P.montgomeryNegInv_lt + Mont64x8Field.negInv_mul_q (by rw [Mont64x8Field.q_toNat]; exact hxq) + Mont64x8Field.two_mul_q_lt).2.2 + rw [Mont64x8Field.q_toNat] at hmod + have hcast := (ZMod.natCast_eq_natCast_iff _ _ _).2 hmod + rw [Nat.cast_mul, Nat.cast_mul] at hcast + rw [← hcast, mul_comm ((2 ^ 256 : ℕ) : ZMod modulus), mul_assoc, + mul_inv_cancel₀ Mont64x8Field.r_ne_zero, mul_one] + +/-- The canonical value of a fast element is its residue divided by `2 ^ 256`. -/ +theorem toField_eq (x : FastField modulus) : + toField x = (x.val.toNat : ZMod modulus) * ((2 ^ 256 : ℕ) : ZMod modulus)⁻¹ := by + rw [toField, toNat, toLimbs8, + mul_cast x.val Limbs8.one x.val_bounded Limbs8.one_bounded x.property.2, Limbs8.one_toNat] + simp + +private theorem val_cast (x : FastField modulus) : + (x.val.toNat : ZMod modulus) = toField x * ((2 ^ 256 : ℕ) : ZMod modulus) := by + rw [toField_eq, mul_assoc, inv_mul_cancel₀ Mont64x8Field.r_ne_zero, mul_one] + +theorem toNat_lt (x : FastField modulus) : toNat x < modulus := by + have := (mul_spec _ _ _ _ P.modulusLimbs_bounded x.val_bounded Limbs8.one_bounded + P.montgomeryNegInv_lt Mont64x8Field.negInv_mul_q x.val_lt Mont64x8Field.two_mul_q_lt).2.1 + rwa [Mont64x8Field.q_toNat] at this + +private theorem ofCanonicalNat_val_cast {n : ℕ} (h : n < modulus) : + ((ofCanonicalNat n h).val.toNat : ZMod modulus) = + (n : ZMod modulus) * ((2 ^ 256 : ℕ) : ZMod modulus) := by + have hR2 : ((((2 ^ 256) ^ 2 : ℕ) % modulus : ℕ) : ZMod modulus) + = ((2 ^ 256 : ℕ) : ZMod modulus) * ((2 ^ 256 : ℕ) : ZMod modulus) := by + rw [ZMod.natCast_mod, pow_two, Nat.cast_mul] + rw [ofCanonicalNat, + mul_cast _ _ (Limbs8.ofNat_bounded n) P.r2ModModulus_bounded + (by rw [Limbs8.ofNat_toNat, Nat.mod_eq_of_lt (h.trans Mont64x8Field.modulus_lt)]; exact h), + Limbs8.ofNat_toNat, Nat.mod_eq_of_lt (h.trans Mont64x8Field.modulus_lt), + P.r2ModModulus_toNat, hR2, mul_assoc, mul_assoc, + mul_inv_cancel₀ Mont64x8Field.r_ne_zero, mul_one] + +@[simp] +theorem toField_ofCanonicalNat {n : ℕ} (h : n < modulus) : + toField (ofCanonicalNat n h) = (n : ZMod modulus) := by + rw [toField_eq, ofCanonicalNat_val_cast h, mul_assoc, + mul_inv_cancel₀ Mont64x8Field.r_ne_zero, mul_one] + +@[simp] +theorem toNat_ofCanonicalNat {n : ℕ} (h : n < modulus) : toNat (ofCanonicalNat n h) = n := + Montgomery.natCast_inj_of_lt (toField_ofCanonicalNat h) (toNat_lt _) h + +/-- Converting from the canonical field to fast form and back is the identity. -/ +@[simp] +theorem toField_ofField (x : ZMod modulus) : toField (ofField x) = x := by + simp [ofField, toField_ofCanonicalNat] + +/-- Converting from fast form to the canonical field and back is the identity. -/ +@[simp] +theorem ofField_toField (x : FastField modulus) : ofField (toField x) = x := by + apply Subtype.ext + have hlt : ∀ y : FastField modulus, y.val.toNat < 2 ^ 256 := fun y => + Limbs8.toNat_lt y.val_bounded + have hval : ((ofField (toField x)).val.toNat : ZMod modulus) = (x.val.toNat : ZMod modulus) := by + rw [val_cast, val_cast, toField_ofField] + have hnat : (ofField (toField x)).val.toNat = x.val.toNat := + Montgomery.natCast_inj_of_lt hval (ofField (toField x)).property.2 x.property.2 + exact Limbs8.ext_of_toNat (ofField (toField x)).val_bounded x.val_bounded hnat + +/-- The canonical-field interpretation distinguishes fast values. -/ +theorem toField_injective : Function.Injective (toField (modulus := modulus)) := + Function.LeftInverse.injective ofField_toField + +/-! ### Field operations -/ + +@[simp] +theorem toField_zero : toField (0 : FastField modulus) = 0 := by + rw [toField_eq, zero_def, zero] + simp [Limbs8.zero_toNat] + +@[simp] +theorem toField_one : toField (1 : FastField modulus) = 1 := by + rw [toField_eq, one_def, one] + simp only [P.rModModulus_toNat, ZMod.natCast_mod] + exact mul_inv_cancel₀ Mont64x8Field.r_ne_zero + +@[simp] +theorem toField_add (x y : FastField modulus) : toField (x + y) = toField x + toField y := by + rw [toField_eq, toField_eq x, toField_eq y, add_def, add] + have h := add_toNat P.modulusLimbs x.val y.val P.modulusLimbs_bounded x.val_bounded + y.val_bounded Mont64x8Field.two_mul_q_lt x.val_lt y.val_lt + rw [h, Mont64x8Field.q_toNat, ZMod.natCast_mod, Nat.cast_add] + ring + +@[simp] +theorem toField_sub (x y : FastField modulus) : toField (x - y) = toField x - toField y := by + rw [toField_eq, toField_eq x, toField_eq y, sub_def, sub] + have h := sub_toNat P.modulusLimbs x.val y.val P.modulusLimbs_bounded x.val_bounded + y.val_bounded Mont64x8Field.two_mul_q_lt x.val_lt y.val_lt + rw [h, Mont64x8Field.q_toNat, ZMod.natCast_mod, Nat.cast_add, + Nat.cast_sub (le_of_lt y.property.2), ZMod.natCast_self] + ring + +@[simp] +theorem toField_neg (x : FastField modulus) : toField (-x) = -toField x := by + rw [toField_eq, toField_eq x, neg_def, neg] + have h := neg_toNat P.modulusLimbs x.val P.modulusLimbs_bounded x.val_bounded + Mont64x8Field.two_mul_q_lt x.val_lt + rw [h, Mont64x8Field.q_toNat, ZMod.natCast_mod, Nat.cast_sub (le_of_lt x.property.2), + ZMod.natCast_self] + ring + +@[simp] +theorem toField_mul (x y : FastField modulus) : toField (x * y) = toField x * toField y := by + rw [toField_eq, toField_eq x, toField_eq y, mul_def, mul, + mul_cast x.val y.val x.val_bounded y.val_bounded x.property.2] + ring + +private theorem mul_assoc' (x y z : FastField modulus) : x * y * z = x * (y * z) := by + apply toField_injective + rw [toField_mul, toField_mul, toField_mul, toField_mul] + ring + +private theorem pow_succ_field (x : FastField modulus) (n : ℕ) : pow x (n + 1) = pow x n * x := by + unfold pow + letI : Semigroup (FastField modulus) := { mul, mul_assoc := mul_assoc' } + exact npowBinRec_succ n x + +@[simp] +theorem toField_square (x : FastField modulus) : toField (square x) = toField x * toField x := by + simp only [square_def, toField_mul] + +@[simp] +theorem toField_pow (x : FastField modulus) (n : ℕ) : toField (pow x n) = toField x ^ n := by + induction n with + | zero => + unfold pow + rw [npowBinRec_zero, toField_one] + simp + | succ n ih => rw [pow_succ_field, toField_mul, ih, _root_.pow_succ] + +/-- Fermat-style inversion in `ZMod modulus`. -/ +private theorem inv_eq_pow {a : ZMod modulus} (ha : a ≠ 0) : a⁻¹ = a ^ (modulus - 2) := by + have hcard : Fintype.card (ZMod modulus) = modulus := ZMod.card modulus + have h1 : a ^ (modulus - 1) = 1 := by + have h := FiniteField.pow_card_sub_one_eq_one a ha + rw [hcard] at h + exact h + have hmul : a * a ^ (modulus - 2) = 1 := by + rw [← pow_succ'] + show a ^ (modulus - 2 + 1) = 1 + have hsucc : modulus - 2 + 1 = modulus - 1 := by + have := P.two_lt_modulus + omega + rw [hsucc] + exact h1 + exact (eq_inv_of_mul_eq_one_left (by rwa [mul_comm])).symm + +@[simp] +theorem toField_inv (x : FastField modulus) : toField x⁻¹ = (toField x)⁻¹ := by + simp only [inv_def, inv, toField_pow] + by_cases hx : toField x = 0 + · rw [hx, inv_zero, zero_pow] + have := P.two_lt_modulus + omega + · rw [inv_eq_pow hx] + +@[simp] +theorem toField_div (x y : FastField modulus) : toField (x / y) = toField x / toField y := by + simp only [div_def, toField_mul, toField_inv] + rfl + +@[simp] +theorem toField_natCast (n : ℕ) : toField (n : FastField modulus) = (n : ZMod modulus) := by + change toField (ofNat modulus n) = (n : ZMod modulus) + rw [ofNat, toField_ofCanonicalNat, ZMod.natCast_eq_natCast_iff] + exact Nat.mod_modEq _ _ + +@[simp] +theorem toField_intCast (n : Int) : toField (n : FastField modulus) = (n : ZMod modulus) := by + change toField (ofField n) = (n : ZMod modulus) + rw [toField_ofField] + +@[simp] +theorem toField_nsmul (n : ℕ) (x : FastField modulus) : toField (n • x) = n • toField x := by + change toField ((n : FastField modulus) * x) = n • toField x + rw [toField_mul, toField_natCast, nsmul_eq_mul] + +@[simp] +theorem toField_zsmul (n : Int) (x : FastField modulus) : toField (n • x) = n • toField x := by + change toField ((n : FastField modulus) * x) = n • toField x + rw [toField_mul, toField_intCast, zsmul_eq_mul] + +@[simp] +theorem toField_npow (x : FastField modulus) (n : ℕ) : toField (x ^ n) = toField x ^ n := by + change toField (pow x n) = toField x ^ n + rw [toField_pow] + +@[simp] +theorem toField_zpow (x : FastField modulus) (n : Int) : toField (x ^ n) = toField x ^ n := by + cases n with + | ofNat n => + change toField (pow x n) = toField x ^ (Int.ofNat n) + rw [toField_pow] + exact (zpow_natCast (toField x) n).symm + | negSucc n => + change toField (pow (inv x) (n + 1)) = toField x ^ (Int.negSucc n) + have hinv : toField (inv x) = (toField x)⁻¹ := by + change toField x⁻¹ = (toField x)⁻¹ + rw [toField_inv] + rw [toField_pow, hinv, zpow_negSucc, inv_pow] + +@[simp] +theorem toField_nnratCast (q : ℚ≥0) : toField (q : FastField modulus) = (q : ZMod modulus) := by + change toField (ofField (q : ZMod modulus)) = (q : ZMod modulus) + rw [toField_ofField] + +@[simp] +theorem toField_ratCast (q : ℚ) : toField (q : FastField modulus) = (q : ZMod modulus) := by + change toField (ofField (q : ZMod modulus)) = (q : ZMod modulus) + rw [toField_ofField] + +@[simp] +theorem toField_nnqsmul (q : ℚ≥0) (x : FastField modulus) : toField (q • x) = q • toField x := by + change toField (ofField (q • toField x)) = q • toField x + rw [toField_ofField] + +@[simp] +theorem toField_qsmul (q : ℚ) (x : FastField modulus) : toField (q • x) = q • toField x := by + change toField (ofField (q • toField x)) = q • toField x + rw [toField_ofField] + +/-! ### Algebraic structure -/ + +/-- Ring equivalence between the fast Montgomery representation and the canonical field. -/ +def ringEquiv (modulus : ℕ) [P : Mont64x8Field modulus] : + FastField modulus ≃+* ZMod modulus where + toFun := toField + invFun := ofField + left_inv := ofField_toField + right_inv := toField_ofField + map_add' := toField_add + map_mul' := toField_mul + +@[simp] +theorem ringEquiv_apply {x : FastField modulus} : ringEquiv modulus x = toField x := rfl + +@[simp] +theorem ringEquiv_symm_apply {x : ZMod modulus} : (ringEquiv modulus).symm x = ofField x := rfl + +/-- Field instance transferred from the canonical field through `toField`. -/ +instance instField : _root_.Field (FastField modulus) := by + apply toField_injective.field toField <;> simp + +/-- A fast eight-limb field is non-binary. -/ +instance instNonBinaryField : NonBinaryField (FastField modulus) where + char_neq_2 := by + intro h + apply Mont64x8Field.two_ne_zero (modulus := modulus) + calc + _ = toField ((2 : ℕ) : FastField modulus) := (toField_natCast 2).symm + _ = toField (0 : FastField modulus) := congrArg toField h + _ = 0 := toField_zero + +end Bridge + +end FastField + +end + +end Native64x8 +end Montgomery diff --git a/CompElliptic/Vendor/CompPoly/Montgomery/Native64x8Mul.lean b/CompElliptic/Vendor/CompPoly/Montgomery/Native64x8Mul.lean new file mode 100644 index 0000000..2cbcfe3 --- /dev/null +++ b/CompElliptic/Vendor/CompPoly/Montgomery/Native64x8Mul.lean @@ -0,0 +1,312 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Gregor Mitscha-Baude +-/ + +import CompElliptic.Vendor.CompPoly.Montgomery.Native64x8 +import Mathlib.Tactic.Linarith + +/-! +# Correctness of eight-limb CIOS Montgomery multiplication + +A CIOS round is the composition of an accumulation and a reduction step: + +* `mulAccum a bi t` accumulates `a * bi` into the accumulator without losing information, + `⟦mulAccum a bi t⟧ = ⟦t⟧ + ⟦a⟧ * bi`; +* `mulReduce q negInv s` adds the multiple `m * q` of the modulus that cancels the low limb + and drops that limb, `2 ^ 32 * ⟦mulReduce q negInv s⟧ = ⟦s⟧ + m * q`. + +Composing them gives the round invariant `2 ^ 32 * ⟦mulRound⟧ = ⟦t⟧ + ⟦a⟧ * bi + m * q`, and +folding eight rounds gives `2 ^ 256 * ⟦t₈⟧ = ⟦a⟧ * ⟦b⟧ + M * q`, so the accumulator is the +Montgomery product up to the final conditional subtraction. + +## Main results + +* `mulAccum_spec`, `mulReduce_spec` — the two halves of a round +* `mulRound_spec` — the round invariant, with limb bounds and the `2 * q` bound +* `mul_spec` — `mul` is canonical and satisfies `2 ^ 256 * ⟦mul a b⟧ ≡ ⟦a⟧ * ⟦b⟧ [MOD q]` +-/ + +namespace Montgomery +namespace Native64x8 + +/-! ### Arithmetic helpers -/ + +/-- Scalar multiplication distributes over a limb recomposition. -/ +theorem sum8_mul (a0 a1 a2 a3 a4 a5 a6 a7 b : ℕ) : + (a0 + 2 ^ 32 * a1 + 2 ^ 64 * a2 + 2 ^ 96 * a3 + 2 ^ 128 * a4 + 2 ^ 160 * a5 + + 2 ^ 192 * a6 + 2 ^ 224 * a7) * b = + a0 * b + 2 ^ 32 * (a1 * b) + 2 ^ 64 * (a2 * b) + 2 ^ 96 * (a3 * b) + + 2 ^ 128 * (a4 * b) + 2 ^ 160 * (a5 * b) + 2 ^ 192 * (a6 * b) + 2 ^ 224 * (a7 * b) := by + ring + +/-- Scalar multiplication distributes over a limb recomposition, from the left. -/ +theorem mul_sum8 (b a0 a1 a2 a3 a4 a5 a6 a7 : ℕ) : + b * (a0 + 2 ^ 32 * a1 + 2 ^ 64 * a2 + 2 ^ 96 * a3 + 2 ^ 128 * a4 + 2 ^ 160 * a5 + + 2 ^ 192 * a6 + 2 ^ 224 * a7) = + b * a0 + 2 ^ 32 * (b * a1) + 2 ^ 64 * (b * a2) + 2 ^ 96 * (b * a3) + + 2 ^ 128 * (b * a4) + 2 ^ 160 * (b * a5) + 2 ^ 192 * (b * a6) + 2 ^ 224 * (b * a7) := by + ring + +/-- The low limb of a bounded value is its residue modulo `2 ^ 32`. -/ +theorem Limbs8.toNat_mod (x : Limbs8) (h : x.l0.toNat < 2 ^ 32) : + x.toNat % 2 ^ 32 = x.l0.toNat := by + simp only [Limbs8.toNat] + omega + +/-- The Montgomery multiplier makes the low limb of the reduction vanish. -/ +private theorem montM_low_zero {s negInv Q q0 w u : ℕ} (hs : s < 2 ^ 32) + (hq0 : Q % 2 ^ 32 = q0) (hnq : negInv * Q % 2 ^ 32 = 2 ^ 32 - 1) + (h : w + 2 ^ 32 * u = s + s * negInv % 2 ^ 32 * q0) (hw : w < 2 ^ 32) : w = 0 := by + have hdvd : 2 ^ 32 ∣ s + s % 2 ^ 32 * negInv % 2 ^ 32 * Q := + Montgomery.dvd_add (2 ^ 32) Q negInv (by norm_num) hnq s + rw [Nat.mod_eq_of_lt hs] at hdvd + have hq0' : q0 % 2 ^ 32 = Q % 2 ^ 32 := by omega + have hcong : (s + s * negInv % 2 ^ 32 * q0) % 2 ^ 32 + = (s + s * negInv % 2 ^ 32 * Q) % 2 ^ 32 := + Nat.ModEq.add_left s (Nat.ModEq.mul_left _ hq0') + obtain ⟨c, hc⟩ := hdvd + omega + +/-- The accumulator stays below `2 * q` from one round to the next. -/ +private theorem round_bound {T A bi m Q Tp : ℕ} (h : 2 ^ 32 * Tp = T + A * bi + m * Q) + (hT : T < 2 * Q) (hA : A < Q) (hbi : bi < 2 ^ 32) (hm : m < 2 ^ 32) : Tp < 2 * Q := by + nlinarith [h, hT, hA, hbi, hm] + +/-- Assembling the limb chain of the accumulation step. -/ +private theorem accum_assemble {Slow Tlow T Abi k7 t8 s8 : ℕ} (hT : T = Tlow + 2 ^ 256 * t8) + (hchain : Slow + 2 ^ 256 * k7 = Tlow + Abi) (hs8 : s8 = t8 + k7) : + Slow + 2 ^ 256 * s8 = T + Abi := by + omega + +/-- Assembling the limb chain of the reduction step. -/ +private theorem reduce_assemble {V W Slow mQ S u7 v8 u8 s8 : ℕ} + (hS : S = Slow + 2 ^ 256 * s8) (hW : W = 2 ^ 32 * V) + (hchain : W + 2 ^ 256 * u7 = Slow + mQ) (hv8 : v8 + 2 ^ 32 * u8 = s8 + u7) : + 2 ^ 32 * (V + 2 ^ 224 * v8 + 2 ^ 256 * u8) = S + mQ := by + omega + +/-! ### The accumulator -/ + +theorem State9.zero_bounded : State9.zero.Bounded := by + simp only [State9.Bounded, State9.toLimbs8, State9.zero, Limbs8.Bounded, UInt64.toNat_zero] + norm_num + +theorem State9.zero_toNat : State9.zero.toNat = 0 := by + simp only [State9.toNat, State9.toLimbs8, State9.zero, Limbs8.toNat, UInt64.toNat_zero] + +/-- The value of an accumulator in terms of its limbs. -/ +theorem State9.toNat_eq (t : State9) : + t.toNat = t.t0.toNat + 2 ^ 32 * t.t1.toNat + 2 ^ 64 * t.t2.toNat + 2 ^ 96 * t.t3.toNat + + 2 ^ 128 * t.t4.toNat + 2 ^ 160 * t.t5.toNat + 2 ^ 192 * t.t6.toNat + + 2 ^ 224 * t.t7.toNat + 2 ^ 256 * t.t8.toNat := by + simp only [State9.toNat, State9.toLimbs8, Limbs8.toNat] + +/-! ### The accumulation step -/ + +/-- `mulAccum` adds `a * bi` to the accumulator exactly: the carry out of the top limb is +retained in the head limb. -/ +theorem mulAccum_spec (a : Limbs8) (bi : UInt64) (t : State9) (ha : a.Bounded) + (ht : t.Bounded) (hbi : bi.toNat < 2 ^ 32) : + (mulAccum a bi t).toLimbs8.Bounded ∧ (mulAccum a bi t).t8.toNat < 2 ^ 33 ∧ + (mulAccum a bi t).toNat = t.toNat + a.toNat * bi.toNat := by + obtain ⟨⟨ht0, ht1, ht2, ht3, ht4, ht5, ht6, ht7⟩, ht8⟩ := ht + obtain ⟨ha0, ha1, ha2, ha3, ha4, ha5, ha6, ha7⟩ := ha + obtain ⟨s0, k0, ds0, dk0, e0, gk0, hs0⟩ := + mac_spec t.t0 a.l0 bi (0 : UInt64) ht0 ha0 hbi (by decide) + obtain ⟨s1, k1, ds1, dk1, e1, gk1, hs1⟩ := + mac_spec t.t1 a.l1 bi _ ht1 ha1 hbi (dk0 ▸ gk0) + obtain ⟨s2, k2, ds2, dk2, e2, gk2, hs2⟩ := + mac_spec t.t2 a.l2 bi _ ht2 ha2 hbi (dk1 ▸ gk1) + obtain ⟨s3, k3, ds3, dk3, e3, gk3, hs3⟩ := + mac_spec t.t3 a.l3 bi _ ht3 ha3 hbi (dk2 ▸ gk2) + obtain ⟨s4, k4, ds4, dk4, e4, gk4, hs4⟩ := + mac_spec t.t4 a.l4 bi _ ht4 ha4 hbi (dk3 ▸ gk3) + obtain ⟨s5, k5, ds5, dk5, e5, gk5, hs5⟩ := + mac_spec t.t5 a.l5 bi _ ht5 ha5 hbi (dk4 ▸ gk4) + obtain ⟨s6, k6, ds6, dk6, e6, gk6, hs6⟩ := + mac_spec t.t6 a.l6 bi _ ht6 ha6 hbi (dk5 ▸ gk5) + obtain ⟨s7, k7, ds7, dk7, e7, gk7, hs7⟩ := + mac_spec t.t7 a.l7 bi _ ht7 ha7 hbi (dk6 ▸ gk6) + simp only [dk0, dk1, dk2, dk3, dk4, dk5, dk6, UInt64.toNat_zero] at e0 e1 e2 e3 e4 e5 e6 e7 + have hchain := carry_chain_sum e0 e1 e2 e3 e4 e5 e6 e7 + rw [← sum8_mul, ← Limbs8.toNat, Nat.add_zero] at hchain + have hk7 : ∀ x : UInt64, x.toNat = k7 → (t.t8 + x).toNat = t.t8.toNat + k7 := by + intro x hx + rw [UInt64.toNat_add, hx, Nat.mod_eq_of_lt (by omega)] + refine ⟨⟨?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_⟩, ?_, ?_⟩ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · simp only [mulAccum] + rw [hk7 _ dk7] + omega + · rw [State9.toNat_eq, State9.toNat_eq t] + simp only [mulAccum, ds0, ds1, ds2, ds3, ds4, ds5, ds6, ds7] + rw [hk7 _ dk7] + exact accum_assemble rfl hchain rfl + +/-! ### The reduction step -/ + +/-- `mulReduce` adds the multiple of the modulus that cancels the low limb, and the division +by `2 ^ 32` performed by dropping that limb is exact. -/ +theorem mulReduce_spec (q : Limbs8) (negInv : UInt64) (s : State9) (hq : q.Bounded) + (hs : s.toLimbs8.Bounded) (hs8 : s.t8.toNat < 2 ^ 33) (hn : negInv.toNat < 2 ^ 32) + (hnq : negInv.toNat * q.toNat % 2 ^ 32 = 2 ^ 32 - 1) : + (mulReduce q negInv s).Bounded ∧ + 2 ^ 32 * (mulReduce q negInv s).toNat = + s.toNat + (montM s.t0 negInv).toNat * q.toNat := by + obtain ⟨hs0, hs1, hs2, hs3, hs4, hs5, hs6, hs7⟩ := hs + obtain ⟨hq0, hq1, hq2, hq3, hq4, hq5, hq6, hq7⟩ := hq + have hmv : (montM s.t0 negInv).toNat = s.t0.toNat * negInv.toNat % 2 ^ 32 := + montM_toNat _ _ hs0 hn + have hmv_lt : (montM s.t0 negInv).toNat < 2 ^ 32 := montM_lt _ _ + obtain ⟨w0, u0, dw0, du0, f0, gu0, hw0⟩ := + mac_spec s.t0 (montM s.t0 negInv) q.l0 (0 : UInt64) hs0 hmv_lt hq0 (by decide) + obtain ⟨w1, u1, dw1, du1, f1, gu1, hw1⟩ := + mac_spec s.t1 (montM s.t0 negInv) q.l1 _ hs1 hmv_lt hq1 (du0 ▸ gu0) + obtain ⟨w2, u2, dw2, du2, f2, gu2, hw2⟩ := + mac_spec s.t2 (montM s.t0 negInv) q.l2 _ hs2 hmv_lt hq2 (du1 ▸ gu1) + obtain ⟨w3, u3, dw3, du3, f3, gu3, hw3⟩ := + mac_spec s.t3 (montM s.t0 negInv) q.l3 _ hs3 hmv_lt hq3 (du2 ▸ gu2) + obtain ⟨w4, u4, dw4, du4, f4, gu4, hw4⟩ := + mac_spec s.t4 (montM s.t0 negInv) q.l4 _ hs4 hmv_lt hq4 (du3 ▸ gu3) + obtain ⟨w5, u5, dw5, du5, f5, gu5, hw5⟩ := + mac_spec s.t5 (montM s.t0 negInv) q.l5 _ hs5 hmv_lt hq5 (du4 ▸ gu4) + obtain ⟨w6, u6, dw6, du6, f6, gu6, hw6⟩ := + mac_spec s.t6 (montM s.t0 negInv) q.l6 _ hs6 hmv_lt hq6 (du5 ▸ gu5) + obtain ⟨w7, u7, dw7, du7, f7, gu7, hw7⟩ := + mac_spec s.t7 (montM s.t0 negInv) q.l7 _ hs7 hmv_lt hq7 (du6 ▸ gu6) + simp only [du0, du1, du2, du3, du4, du5, du6, UInt64.toNat_zero, hmv] at f0 f1 f2 f3 f4 f5 f6 f7 + have hw0zero : w0 = 0 := + montM_low_zero hs0 (Limbs8.toNat_mod q hq0) hnq (by rw [Nat.add_zero] at f0; exact f0) hw0 + rw [hw0zero] at f0 + have hchain := carry_chain_sum f0 f1 f2 f3 f4 f5 f6 f7 + rw [← mul_sum8, ← Limbs8.toNat, Nat.add_zero] at hchain + obtain ⟨v8, u8, dv8, du8, g8, gu8, hv8⟩ := + adc_spec_wide s.t8 _ (0 : UInt64) hs8 (du7 ▸ gu7) (by decide) + rw [UInt64.toNat_zero, Nat.add_zero, du7] at g8 + refine ⟨⟨⟨?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_⟩, ?_⟩, ?_⟩ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · exact macLo_lt _ _ _ _ + · exact adcLo_lt _ _ _ + · simp only [mulReduce] + rw [du8] + omega + · rw [State9.toNat_eq, State9.toNat_eq s, hmv] + simp only [mulReduce, dw1, dw2, dw3, dw4, dw5, dw6, dw7, dv8, du8] + exact reduce_assemble rfl (by ring) hchain g8 + +/-! ### The round invariant -/ + +/-- The CIOS round invariant: one round accumulates `a * bi` and cancels the low limb against +a multiple `m` of the modulus; the accumulator stays limbwise bounded and below `2 * q`. -/ +theorem mulRound_spec (q a : Limbs8) (negInv bi : UInt64) (t : State9) + (hq : q.Bounded) (ha : a.Bounded) (ht : t.Bounded) (hbi : bi.toNat < 2 ^ 32) + (hn : negInv.toNat < 2 ^ 32) (hnq : negInv.toNat * q.toNat % 2 ^ 32 = 2 ^ 32 - 1) + (haq : a.toNat < q.toNat) (htq : t.toNat < 2 * q.toNat) : + (mulRound q negInv a bi t).Bounded ∧ (mulRound q negInv a bi t).toNat < 2 * q.toNat ∧ + ∃ m : ℕ, m < 2 ^ 32 ∧ + 2 ^ 32 * (mulRound q negInv a bi t).toNat = + t.toNat + a.toNat * bi.toNat + m * q.toNat := by + obtain ⟨hab, ha8, hav⟩ := mulAccum_spec a bi t ha ht hbi + obtain ⟨hrb, hrv⟩ := mulReduce_spec q negInv (mulAccum a bi t) hq hab ha8 hn hnq + rw [hav] at hrv + have hinv : 2 ^ 32 * (mulRound q negInv a bi t).toNat = + t.toNat + a.toNat * bi.toNat + (montM (mulAccum a bi t).t0 negInv).toNat * q.toNat := by + rw [mulRound] + exact hrv + exact ⟨hrb, round_bound hinv htq haq hbi (montM_lt _ _), _, montM_lt _ _, hinv⟩ + +/-! ### The eight-round fold -/ + +/-- Folding the eight round invariants into a single Montgomery identity. -/ +private theorem fold8 {T1 T2 T3 T4 T5 T6 T7 T8 A Q : ℕ} + {B0 B1 B2 B3 B4 B5 B6 B7 m0 m1 m2 m3 m4 m5 m6 m7 : ℕ} + (h0 : 2 ^ 32 * T1 = 0 + A * B0 + m0 * Q) + (h1 : 2 ^ 32 * T2 = T1 + A * B1 + m1 * Q) + (h2 : 2 ^ 32 * T3 = T2 + A * B2 + m2 * Q) + (h3 : 2 ^ 32 * T4 = T3 + A * B3 + m3 * Q) + (h4 : 2 ^ 32 * T5 = T4 + A * B4 + m4 * Q) + (h5 : 2 ^ 32 * T6 = T5 + A * B5 + m5 * Q) + (h6 : 2 ^ 32 * T7 = T6 + A * B6 + m6 * Q) + (h7 : 2 ^ 32 * T8 = T7 + A * B7 + m7 * Q) : + 2 ^ 256 * T8 = + (A * B0 + 2 ^ 32 * (A * B1) + 2 ^ 64 * (A * B2) + 2 ^ 96 * (A * B3) + + 2 ^ 128 * (A * B4) + 2 ^ 160 * (A * B5) + 2 ^ 192 * (A * B6) + 2 ^ 224 * (A * B7)) + + (m0 * Q + 2 ^ 32 * (m1 * Q) + 2 ^ 64 * (m2 * Q) + 2 ^ 96 * (m3 * Q) + + 2 ^ 128 * (m4 * Q) + 2 ^ 160 * (m5 * Q) + 2 ^ 192 * (m6 * Q) + + 2 ^ 224 * (m7 * Q)) := by + omega + +/-- The final conditional subtraction of `mul`, given the folded Montgomery identity. -/ +private theorem mul_finish (q : Limbs8) (t : State9) {A M : ℕ} (hq : q.Bounded) + (hb : t.Bounded) (hlt : t.toNat < 2 * q.toNat) (hq2 : 2 * q.toNat < 2 ^ 256) + (hfold : 2 ^ 256 * t.toNat = A + M * q.toNat) : + (condSub q t.toLimbs8).Bounded ∧ (condSub q t.toLimbs8).toNat < q.toNat ∧ + 2 ^ 256 * (condSub q t.toLimbs8).toNat ≡ A [MOD q.toNat] := by + have hlimbs : t.toLimbs8.toNat = t.toNat := by + have h1 := Limbs8.toNat_lt hb.1 + have h2 : t.toNat = t.toLimbs8.toNat + 2 ^ 256 * t.t8.toNat := rfl + omega + have hmod : 2 ^ 256 * t.toNat ≡ A [MOD q.toNat] := by + unfold Nat.ModEq + rw [hfold, Nat.add_mul_mod_self_right] + refine ⟨condSub_bounded _ _ hb.1, condSub_lt q _ hq hb.1 (by omega), ?_⟩ + rw [condSub_toNat q _ hq hb.1, hlimbs] + split + · exact hmod + · refine (Nat.ModEq.mul_left _ ?_).trans hmod + exact (Nat.modEq_iff_dvd' (by omega)).2 ⟨1, by omega⟩ + +/-- Montgomery multiplication is canonical and computes `a * b * (2 ^ 256)⁻¹ mod q`. -/ +theorem mul_spec (q : Limbs8) (negInv : UInt64) (a b : Limbs8) + (hq : q.Bounded) (ha : a.Bounded) (hb : b.Bounded) (hn : negInv.toNat < 2 ^ 32) + (hnq : negInv.toNat * q.toNat % 2 ^ 32 = 2 ^ 32 - 1) + (haq : a.toNat < q.toNat) (hq2 : 2 * q.toNat < 2 ^ 256) : + (mul q negInv a b).Bounded ∧ (mul q negInv a b).toNat < q.toNat ∧ + 2 ^ 256 * (mul q negInv a b).toNat ≡ a.toNat * b.toNat [MOD q.toNat] := by + obtain ⟨hb0, hb1, hb2, hb3, hb4, hb5, hb6, hb7⟩ := hb + obtain ⟨B1, L1, m0, hm0, r0⟩ := + mulRound_spec q a negInv b.l0 State9.zero hq ha State9.zero_bounded hb0 hn hnq haq + (by rw [State9.zero_toNat]; omega) + obtain ⟨B2, L2, m1, hm1, r1⟩ := + mulRound_spec q a negInv b.l1 _ hq ha B1 hb1 hn hnq haq + L1 + obtain ⟨B3, L3, m2, hm2, r2⟩ := + mulRound_spec q a negInv b.l2 _ hq ha B2 hb2 hn hnq haq + L2 + obtain ⟨B4, L4, m3, hm3, r3⟩ := + mulRound_spec q a negInv b.l3 _ hq ha B3 hb3 hn hnq haq + L3 + obtain ⟨B5, L5, m4, hm4, r4⟩ := + mulRound_spec q a negInv b.l4 _ hq ha B4 hb4 hn hnq haq + L4 + obtain ⟨B6, L6, m5, hm5, r5⟩ := + mulRound_spec q a negInv b.l5 _ hq ha B5 hb5 hn hnq haq + L5 + obtain ⟨B7, L7, m6, hm6, r6⟩ := + mulRound_spec q a negInv b.l6 _ hq ha B6 hb6 hn hnq haq + L6 + obtain ⟨B8, L8, m7, hm7, r7⟩ := + mulRound_spec q a negInv b.l7 _ hq ha B7 hb7 hn hnq haq + L7 + rw [State9.zero_toNat] at r0 + have hfold := fold8 r0 r1 r2 r3 r4 r5 r6 r7 + rw [← mul_sum8, ← sum8_mul, ← Limbs8.toNat] at hfold + simp only [mul] + exact mul_finish q _ hq B8 L8 hq2 hfold + +end Native64x8 +end Montgomery diff --git a/CompElliptic/Vendor/CompPoly/Montgomery/Pasta.lean b/CompElliptic/Vendor/CompPoly/Montgomery/Pasta.lean new file mode 100644 index 0000000..81d9980 --- /dev/null +++ b/CompElliptic/Vendor/CompPoly/Montgomery/Pasta.lean @@ -0,0 +1,52 @@ +/- +Copyright (c) 2026 CompElliptic Contributors. All rights reserved. +Released under the Apache License, Version 2.0, or the MIT license, at your option, +as described in the files LICENSE-APACHE and LICENSE-MIT. +Authors: Gregor Mitscha-Baude +-/ +import CompElliptic.Vendor.CompPoly.Montgomery.Native64x8Field +import CompElliptic.Fields.Pasta + +/-! +# The Pasta base fields as fast eight-limb Montgomery fields + +The per-field data realizing the two Pasta base fields as `Mont64x8Field`s. Unlike the +upstream CompPoly branch, which is self-contained and therefore carries its own Pratt +certificates, this vendored copy reuses the primes and certificates of +`CompElliptic.Fields.Pasta`, so that there is exactly one +`Field (ZMod PALLAS_SCALAR_CARD)` instance in the repository and the Montgomery carrier +bridges directly into the `Fq` of `CompElliptic.Curves.Pasta.Fast.Projective`. + +Note the Pasta naming: `PALLAS_SCALAR_CARD` is the **Vesta base** field size and +`PALLAS_BASE_CARD` is the **Vesta scalar** field size. +-/ + +namespace Montgomery +namespace Native64x8 + +open CompElliptic.Fields.Pasta + +/-- The Vesta base field (`= PALLAS_SCALAR_CARD`) as a fast eight-limb Montgomery field. -/ +instance instVestaBase : Mont64x8Field PALLAS_SCALAR_CARD where + prime := PALLAS_SCALAR_is_prime + modulusLimbs := VestaFq.modulusLimbs + rModModulus := VestaFq.rModModulus + r2ModModulus := VestaFq.r2ModModulus + montgomeryNegInv := VestaFq.negInv + +/-- The Pallas base field (`= PALLAS_BASE_CARD`) as a fast eight-limb Montgomery field. -/ +instance instPallasBase : Mont64x8Field PALLAS_BASE_CARD where + prime := PALLAS_BASE_is_prime + modulusLimbs := PallasFq.modulusLimbs + rModModulus := PallasFq.rModModulus + r2ModModulus := PallasFq.r2ModModulus + montgomeryNegInv := PallasFq.negInv + +/-- The fast Vesta base field carrier, stored as a Montgomery residue. -/ +abbrev VestaBaseFast : Type := FastField PALLAS_SCALAR_CARD + +/-- The fast Pallas base field carrier, stored as a Montgomery residue. -/ +abbrev PallasBaseFast : Type := FastField PALLAS_BASE_CARD + +end Native64x8 +end Montgomery diff --git a/CompElliptic/Vendor/README.md b/CompElliptic/Vendor/README.md new file mode 100644 index 0000000..75d8dbe --- /dev/null +++ b/CompElliptic/Vendor/README.md @@ -0,0 +1,48 @@ +# Vendored material (temporary) + +This directory holds code that belongs to an **upstream dependency** and lives here only until +CompElliptic's pin of that dependency can provide it. Nothing else should be added to it, and +nothing in here may import anything from the rest of CompElliptic *except* what is explicitly +noted below. + +Each subdirectory carries its own deletion criterion. When that criterion is met the +subdirectory is deleted and the import paths are rewritten — the code itself does not change. + +## `CompPoly/Montgomery/` — eight-limb Montgomery field + +A **temporary vendoring** of the [CompPoly](https://github.com/Verified-zkEVM/CompPoly) branch +`fast_multilimb_fields` (commit `b3850f0`), which adds eight-limb (8 × 32-bit packed in +`UInt64`) Montgomery arithmetic for 255-bit prime moduli next to the existing single-word +`Montgomery/Native32*`. It is what makes `CompElliptic/Curves/Pasta/Fast/ProjectiveMontDefs.lean` +fast enough to be worth precompiling. + +**Delete this directory** when CompElliptic's CompPoly pin moves past +[CompPoly#258](https://github.com/Verified-zkEVM/CompPoly/pull/258) (the eight-limb Montgomery +field) and [CompPoly#274](https://github.com/Verified-zkEVM/CompPoly/pull/274). On landing, the +import paths become `CompElliptic.Vendor.CompPoly.Montgomery.X` → +`CompPoly.Fields.Montgomery.X`, and `Pasta.lean` is dropped (see below). + +| File | Upstream original | Change | +|---|---|---| +| `Basic.lean` | `CompPoly/Fields/Montgomery/Basic.lean` | none | +| `Native64x8Defs.lean` | `CompPoly/Fields/Montgomery/Native64x8.lean` (definitions) | split out, `ℕ` → `Nat`, plus the Pasta constants and monomorphic entry points | +| `Native64x8.lean` | `CompPoly/Fields/Montgomery/Native64x8.lean` (theorems) | imports the split-out definitions; two `norm_num` calls dropped (Lean 4.30 vs 4.31 `simp` drift) | +| `Native64x8Mul.lean` | same name | import path; one `norm_num` dropped | +| `Native64x8Field.lean` | same name | import path | +| `Pasta.lean` | `CompPoly/Fields/Pasta/{Basic,Fast}.lean` | **rewritten**: reuses `CompElliptic.Fields.Pasta`'s primes and Pratt certificates instead of vendoring a second copy, so the repo keeps exactly one `Field (ZMod PALLAS_SCALAR_CARD)` instance and the Montgomery carrier bridges directly into `CompElliptic.Curves.Pasta.Fast.Projective.Fq` | + +`Pasta.lean` is the one file here that imports CompElliptic (`CompElliptic.Fields.Pasta`); it is +also the one file that is dropped rather than re-pointed when the vendoring ends, precisely +because upstream carries its own copy of that data. + +The definition/proof split exists for the precompiled lane: `Native64x8Defs.lean` imports +nothing beyond Lean core, so it can be native-compiled (the `FastFieldNative` library, see the +lakefile) without dragging a mathlib import closure through codegen. Keep it that way. +Upstream has since **also** adopted the same definitions/proofs split, so the vendored copy and +upstream have converged on the same shape; the eventual migration is an import-path rewrite +rather than a restructuring. + +Namespaces are deliberately **unchanged** from upstream (`Montgomery.Native64x8`) — that is the +whole point of vendoring these files unmodified, and it is what makes the migration a pure +import-path rewrite. The pinned CompPoly predates the `Montgomery` directory, so nothing +clashes. diff --git a/FastFieldNative.lean b/FastFieldNative.lean new file mode 100644 index 0000000..a3599fc --- /dev/null +++ b/FastFieldNative.lean @@ -0,0 +1,17 @@ +/- +Copyright (c) 2026 CompElliptic Contributors. All rights reserved. +Released under the Apache License, Version 2.0, or the MIT license, at your option, +as described in the files LICENSE-APACHE and LICENSE-MIT. +Authors: Gregor Mitscha-Baude +-/ +import CompElliptic.Vendor.CompPoly.Montgomery.Native64x8Defs +import CompElliptic.Curves.Pasta.Fast.ProjectiveMontDefs + +/-! +# The precompiled native lane + +Root module of the `FastFieldNative` library (see the lakefile): it exists only so that the +emitted `.so` carries an initializer named after the library. Its import closure — the Montgomery +field definitions and the Vesta kernel — is core-only by construction; the proofs about them live +in sibling modules outside the lane. +-/ diff --git a/README.md b/README.md index 7ddb894..1be1fdf 100644 --- a/README.md +++ b/README.md @@ -110,7 +110,14 @@ Early work in progress. Present so far: (`CompElliptic/Encoding.lean`); - a computable Tonelli–Shanks square root for prime fields, soundness and completeness proved, with `pallasBase` / `vestaBase` instances (`CompElliptic/Fields/Sqrt.lean`); -- the compressed Pasta point encoding (`toBytes`) for Pallas and Vesta (`CompElliptic/Encodings/`). +- the compressed Pasta point encoding (`toBytes`) for Pallas and Vesta (`CompElliptic/Encodings/`); +- fast Vesta group arithmetic for computing with points rather than only reasoning about them — + complete Renes–Costello–Batina projective addition, a windowed Pippenger multi-scalar + multiplication and a scalar ladder, each proven to compute the affine group operation it replaces, + together with a natively compilable Montgomery-limb kernel proven against them + (`CompElliptic/Curves/Pasta/Fast/`, not imported by `CompElliptic.lean`). These interfaces are + provisional: they are not guaranteed to remain public, and may be folded into the existing API + or otherwise changed incompatibly. Uses of `sorry` are kept minimal and limited to work-in-progress. The library's general theorems depend only on the standard `propext` / `Classical.choice` / `Quot.sound` axioms. Facts specific to diff --git a/lakefile.lean b/lakefile.lean index e63414f..789703c 100644 --- a/lakefile.lean +++ b/lakefile.lean @@ -9,5 +9,22 @@ package CompElliptic where require CompPoly from git "https://github.com/Verified-zkEVM/CompPoly.git" @ "d458e5cebd364b15660ff8de20ec964dcd52c120" +-- The glob covers modules the root `CompElliptic.lean` deliberately does not import: the fast +-- arithmetic and its vendored field are opt-in, so that `import CompElliptic` stays free of the +-- precompiled lane below. @[default_target] -lean_lib CompElliptic +lean_lib CompElliptic where + globs := #[.andSubmodules `CompElliptic] + +-- Native-compiles the Montgomery arithmetic, which is meant to be run (`#eval`, `native_decide`), +-- not only proven about. Every module here must import nothing outside Lean core: codegen runs over +-- the whole import closure, so one mathlib-side import silently makes the build enormous -- hence +-- the definitions/proofs split, and `scripts/check_native_lane.sh`. `FastFieldNative.lean` exists +-- only because Lean derives a dynlib's `initialize_...` symbol from the library name. +lean_lib FastFieldNative where + precompileModules := true + globs := #[ + .one `FastFieldNative, + .one `CompElliptic.Vendor.CompPoly.Montgomery.Native64x8Defs, + .one `CompElliptic.Curves.Pasta.Fast.ProjectiveMontDefs + ] diff --git a/scripts/check_native_lane.sh b/scripts/check_native_lane.sh new file mode 100755 index 0000000..9e6d2dd --- /dev/null +++ b/scripts/check_native_lane.sh @@ -0,0 +1,40 @@ +#!/usr/bin/env bash +# Check the core-only invariant of the `FastFieldNative` precompiled lane: codegen runs over the +# lane's entire import closure, so a single mathlib-side import makes the build enormous without +# failing loudly. Run from the repository root; exits non-zero on violation. +set -euo pipefail + +cd "$(dirname "$0")/.." + +# The lane, mirroring the `FastFieldNative` globs in lakefile.lean. +LANE=( + "FastFieldNative.lean" + "CompElliptic/Vendor/CompPoly/Montgomery/Native64x8Defs.lean" + "CompElliptic/Curves/Pasta/Fast/ProjectiveMontDefs.lean" +) + +# `FastFieldNative.lean` is the root module; it may import the rest of the lane and nothing else. +status=0 +for f in "${LANE[@]}"; do + if [[ ! -f "$f" ]]; then + echo "MISSING: $f is listed in the lane but does not exist" >&2 + status=1 + continue + fi + while read -r _ mod _; do + [[ -z "${mod:-}" ]] && continue + path="${mod//.//}.lean" + for allowed in "${LANE[@]}"; do + if [[ "$path" == "$allowed" ]]; then + continue 2 + fi + done + echo "VIOLATION: $f imports $mod, which is outside the core-only lane" >&2 + status=1 + done < <(grep -E '^import ' "$f" || true) +done + +if [[ $status -eq 0 ]]; then + echo "FastFieldNative lane is core-only (${#LANE[@]} modules)." +fi +exit $status