From 9348931699fc24c39d1b5aadad71113515f3c2a0 Mon Sep 17 00:00:00 2001 From: Varun Thakore Date: Mon, 29 Jun 2026 20:40:11 +0530 Subject: [PATCH 1/3] feat(fields): add fast Goldilocks arithmetic --- CompPoly/Fields/Goldilocks.lean | 24 +- CompPoly/Fields/Goldilocks/Basic.lean | 45 ++ CompPoly/Fields/Goldilocks/Fast.lean | 17 + .../Fields/Goldilocks/Fast/Arithmetic.lean | 243 +++++++++ CompPoly/Fields/Goldilocks/Fast/Field.lean | 71 +++ CompPoly/Fields/Goldilocks/Fast/Internal.lean | 488 ++++++++++++++++++ .../Fields/Goldilocks/Fast/Reduction.lean | 478 +++++++++++++++++ CompPoly/Fields/Goldilocks/Fast/Theorems.lean | 321 ++++++++++++ 8 files changed, 1670 insertions(+), 17 deletions(-) create mode 100644 CompPoly/Fields/Goldilocks/Basic.lean create mode 100644 CompPoly/Fields/Goldilocks/Fast.lean create mode 100644 CompPoly/Fields/Goldilocks/Fast/Arithmetic.lean create mode 100644 CompPoly/Fields/Goldilocks/Fast/Field.lean create mode 100644 CompPoly/Fields/Goldilocks/Fast/Internal.lean create mode 100644 CompPoly/Fields/Goldilocks/Fast/Reduction.lean create mode 100644 CompPoly/Fields/Goldilocks/Fast/Theorems.lean diff --git a/CompPoly/Fields/Goldilocks.lean b/CompPoly/Fields/Goldilocks.lean index c4d4e55a..2c30a1f6 100644 --- a/CompPoly/Fields/Goldilocks.lean +++ b/CompPoly/Fields/Goldilocks.lean @@ -1,26 +1,16 @@ /- Copyright (c) 2024 ArkLib Contributors. All rights reserved. Released under Apache 2.0 license as described in the file LICENSE. -Authors: Quang Dao +Authors: Quang Dao, Varun Thakore -/ -import CompPoly.Fields.PrattCertificate +import CompPoly.Fields.Goldilocks.Basic +import CompPoly.Fields.Goldilocks.Fast /-! - # Goldilocks prime field `2^{64} - 2^{32} + 1` +# Goldilocks Field - This is the field used in Plonky2/3. +Facade module for the Goldilocks field. It re-exports the canonical `ZMod` model +from `CompPoly.Fields.Goldilocks.Basic` and the native-word implementation from +`CompPoly.Fields.Goldilocks.Fast`. -/ - -namespace Goldilocks - -@[reducible] -def fieldSize : Nat := 2 ^ 64 - 2 ^ 32 + 1 - -abbrev Field := ZMod fieldSize - -theorem is_prime : Nat.Prime fieldSize := by - unfold fieldSize - pratt - -end Goldilocks diff --git a/CompPoly/Fields/Goldilocks/Basic.lean b/CompPoly/Fields/Goldilocks/Basic.lean new file mode 100644 index 00000000..0228e5f1 --- /dev/null +++ b/CompPoly/Fields/Goldilocks/Basic.lean @@ -0,0 +1,45 @@ +/- +Copyright (c) 2024 ArkLib Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Quang Dao, Varun Thakore +-/ + +import CompPoly.Fields.Basic +import CompPoly.Fields.PrattCertificate + +/-! +# Goldilocks Prime Field `2^{64} - 2^{32} + 1` + +This module defines the canonical `ZMod` model of the Goldilocks prime field. +-/ + +namespace Goldilocks +namespace Basic + +/-- The Goldilocks prime modulus, `2^64 - 2^32 + 1`. -/ +@[reducible] +def fieldSize : Nat := 2 ^ 64 - 2 ^ 32 + 1 + +/-- The canonical mathematical Goldilocks field, implemented as integers modulo +`fieldSize`. -/ +abbrev Field := ZMod fieldSize + +/-- The Goldilocks modulus is prime, verified by a Pratt certificate. -/ +theorem is_prime : Nat.Prime fieldSize := by + unfold fieldSize + pratt + +/-- Register primality of `fieldSize` for Mathlib instances such as `ZMod.instField`. -/ +instance : Fact (Nat.Prime fieldSize) := ⟨is_prime⟩ + +/-- The canonical Goldilocks carrier is a field because its modulus is prime. -/ +instance : _root_.Field Field := ZMod.instField fieldSize + +/-- Goldilocks has characteristic different from two. -/ +instance : NonBinaryField Field where + char_neq_2 := by + simpa [Field, fieldSize] using + (by decide : (2 : ZMod (2 ^ 64 - 2 ^ 32 + 1)) ≠ 0) + +end Basic +end Goldilocks diff --git a/CompPoly/Fields/Goldilocks/Fast.lean b/CompPoly/Fields/Goldilocks/Fast.lean new file mode 100644 index 00000000..f5588ee9 --- /dev/null +++ b/CompPoly/Fields/Goldilocks/Fast.lean @@ -0,0 +1,17 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Varun Thakore +-/ + +import CompPoly.Fields.Goldilocks.Fast.Field + +/-! +# Fast Goldilocks Field + +Public entry point for the native-word Goldilocks implementation. + +Importing this module provides the fast `UInt64`-backed `Goldilocks.Fast.Field`, +its arithmetic operations, correctness theorems relating it to +`Goldilocks.Basic.Field` and the transferred field instances. +-/ diff --git a/CompPoly/Fields/Goldilocks/Fast/Arithmetic.lean b/CompPoly/Fields/Goldilocks/Fast/Arithmetic.lean new file mode 100644 index 00000000..60c37ffc --- /dev/null +++ b/CompPoly/Fields/Goldilocks/Fast/Arithmetic.lean @@ -0,0 +1,243 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Varun Thakore +-/ + +import CompPoly.Fields.Goldilocks.Fast.Reduction + +/-! +# Fast Goldilocks Arithmetic + +This module defines the fast native-word Goldilocks carrier, conversions, +arithmetic operations, and notation instances. + +Fast field values are stored as canonical `UInt64` representatives below the +Goldilocks prime `p = 2^64 - 2^32 + 1`. Multiplication uses a 128-bit +intermediate represented as a pair of `UInt64` words `(lo, hi)` and reduces it +with the Goldilocks identity + +`2^64 ≡ 2^32 - 1 (mod p)`. +-/ + +namespace Goldilocks +namespace Fast + +/-- The fast native-word Goldilocks field carrier, stored as a canonical residue. -/ +abbrev Field : Type := { x : UInt64 // x.toNat < Goldilocks.Basic.fieldSize } + +/-- Fast representatives have decidable equality through their `UInt64` value. -/ +instance : DecidableEq Field := inferInstance + +/-- The raw canonical word backing a fast Goldilocks element. -/ +@[inline] +def raw (x : Field) : UInt64 := x.val + +/-- Reduce a native `UInt64` modulo Goldilocks. -/ +@[inline] +def reduceUInt64 (x : UInt64) : Field := + ⟨Reduction.reduceUInt64Raw x, Reduction.reduceUInt64Raw_lt x⟩ + +/-- One-word reduction preserves the represented canonical field element. -/ +@[simp] +theorem reduceUInt64_cast (x : UInt64) : + ((reduceUInt64 x).val.toNat : Goldilocks.Basic.Field) = + (x.toNat : Goldilocks.Basic.Field) := by + exact Reduction.reduceUInt64Raw_cast x + +/-- The zero fast Goldilocks element. -/ +@[inline] +def zero : Field := ⟨0, by decide⟩ + +/-- The one fast Goldilocks element. -/ +@[inline] +def one : Field := ⟨1, by decide⟩ + +/-- Build a fast element from a canonical natural representative. -/ +@[inline] +def ofCanonicalNat (n : Nat) (h : n < Goldilocks.Basic.fieldSize) : Field := + ⟨UInt64.ofNat n, by + have hn : n < UInt64.size := Nat.lt_trans h fieldSize_lt_uint64Size + rw [UInt64.toNat_ofNat'] + rw [Nat.mod_eq_of_lt] + · exact h + · simpa [UInt64.size] using hn⟩ + +/-- Convert a natural number into fast canonical representation. -/ +@[inline] +def ofNat (n : Nat) : Field := + ofCanonicalNat (n % Goldilocks.Basic.fieldSize) (Nat.mod_lt _ fieldSize_pos) + +/-- Convert a 64-bit word into fast canonical representation. -/ +@[inline] +def ofUInt64 (x : UInt64) : Field := + reduceUInt64 x + +/-- Convert from the canonical `ZMod` Goldilocks field into fast canonical form. -/ +@[inline] +def ofField (x : Goldilocks.Basic.Field) : Field := + ofCanonicalNat x.val (ZMod.val_lt x) + +/-- Convert an integer into fast canonical representation. -/ +@[inline] +def ofInt (z : Int) : Field := + ofField (z : Goldilocks.Basic.Field) + +/-- Convert a fast Goldilocks element to its canonical natural representative. -/ +@[inline] +def toNat (x : Field) : Nat := + x.val.toNat + +/-- Convert a fast Goldilocks element to the canonical `ZMod` Goldilocks field. -/ +@[inline] +def toField (x : Field) : Goldilocks.Basic.Field := + (toNat x : Goldilocks.Basic.Field) + +/-- Fast modular addition in canonical form. -/ +@[inline] +def add (x y : Field) : Field := + let lo := x.val + y.val + let carry := decide (lo < x.val) + ⟨Reduction.reduceAddWithCarryRaw lo carry, + Reduction.reduceAddWithCarryRaw_lt lo carry + (Reduction.addWithCarry_bound x.val y.val x.property y.property)⟩ + +/-- Fast modular negation in canonical form. -/ +@[inline] +def neg (x : Field) : Field := + ⟨Reduction.negRaw x.val, Reduction.negRaw_lt x.val x.property⟩ + +/-- Fast modular subtraction in canonical form. -/ +@[inline] +def sub (x y : Field) : Field := + ⟨Reduction.subRaw x.val y.val, Reduction.subRaw_lt x.val y.val x.property y.property⟩ + +/-- Fast modular multiplication in canonical form. -/ +@[inline] +def mul (x y : Field) : Field := + ⟨Reduction.reduceMulRaw x.val y.val, Reduction.reduceMulRaw_lt x.val y.val⟩ + +/-- Fast squaring. -/ +@[inline] +def square (x : Field) : Field := + mul x x + +/-- Repeated squaring: `squareN x n` computes `x^(2^n)`. -/ +@[inline] +def squareN (x : Field) : Nat → Field + | 0 => x + | n + 1 => square (squareN x n) + +/-- Exponentiation over the fast representation using binary exponentiation. -/ +@[inline] +def pow (x : Field) (n : Nat) : Field := + @npowBinRec Field ⟨one⟩ ⟨mul⟩ n x + +/-- Fermat exponent used for inversion in the Goldilocks prime field. -/ +@[inline] +def invExponent : Nat := Goldilocks.Basic.fieldSize - 2 + +/-- Fast modular inversion using an addition chain for `p - 2`. + +For Goldilocks, `p - 2 = 0xFFFFFFFEFFFFFFFF`. The chain builds +`x^(2^31 - 1)`, derives `x^(2^32 - 2)` and `x^(2^32 - 1)`, then combines them as + +`(2^32 - 2) * 2^32 + (2^32 - 1) = p - 2`. +-/ +@[noinline] +def inv (x : Field) : Field := + let t2 := mul (square x) x + let t4 := mul (squareN t2 2) t2 + let t8 := mul (squareN t4 4) t4 + let t16 := mul (squareN t8 8) t8 + let t31 := + mul (squareN t16 15) + (mul (squareN t8 7) + (mul (squareN t4 3) + (mul (square t2) x))) + let t32m2 := square t31 + let t32m1 := mul t32m2 x + mul (squareN t32m2 32) t32m1 + +/-- Division through inversion and fast multiplication. -/ +@[inline] +def div (x y : Field) : Field := + mul x (inv y) + +/-- Use fast zero for standard `0` notation. -/ +instance instZeroField : Zero Field where + zero := zero + +/-- Use fast one for standard `1` notation. -/ +instance instOneField : One Field where + one := one + +/-- Use fast addition for standard `+` notation. -/ +instance instAddField : Add Field where + add := add + +/-- Use fast negation for standard unary `-` notation. -/ +instance instNegField : Neg Field where + neg := neg + +/-- Use fast subtraction for standard `-` notation. -/ +instance instSubField : Sub Field where + sub := sub + +/-- Use fast multiplication for standard `*` notation. -/ +instance instMulField : Mul Field where + mul := mul + +/-- Use fast inversion for standard inverse notation. -/ +instance instInvField : Inv Field where + inv := inv + +/-- Use fast division for standard `/` notation. -/ +instance instDivField : Div Field where + div := div + +/-- Use `ofNat` for natural-number casts into fast Goldilocks. -/ +instance instNatCastField : NatCast Field where + natCast := ofNat + +/-- Interpret integer casts through the canonical Goldilocks field. -/ +instance instIntCastField : IntCast Field where + intCast := ofInt + +/-- Natural scalar multiplication is multiplication by the corresponding fast natural cast. -/ +instance instNatSMulField : SMul Nat Field where + smul n x := (n : Field) * x + +/-- Integer scalar multiplication is multiplication by the corresponding fast integer cast. -/ +instance instIntSMulField : SMul Int Field where + smul n x := (n : Field) * x + +/-- Use fast binary exponentiation for natural powers. -/ +instance instPowFieldNat : Pow Field Nat where + pow := pow + +/-- Use fast natural powers and inversion for integer powers. -/ +instance instPowFieldInt : Pow Field Int where + pow x n := + match n with + | Int.ofNat k => pow x k + | Int.negSucc k => pow (inv x) (k + 1) + +/-- Interpret nonnegative rational casts through the canonical Goldilocks field. -/ +instance instNNRatCastField : NNRatCast Field where + nnratCast q := ofField (q : Goldilocks.Basic.Field) + +/-- Interpret rational casts through the canonical Goldilocks field. -/ +instance instRatCastField : RatCast Field where + ratCast q := ofField (q : Goldilocks.Basic.Field) + +/-- Nonnegative rational scalar multiplication is transported through the canonical field. -/ +instance instNNRatSMulField : SMul ℚ≥0 Field where + smul q x := ofField (q • toField x) + +/-- Rational scalar multiplication is transported through the canonical field. -/ +instance instRatSMulField : SMul ℚ Field where + smul q x := ofField (q • toField x) + +end Fast +end Goldilocks diff --git a/CompPoly/Fields/Goldilocks/Fast/Field.lean b/CompPoly/Fields/Goldilocks/Fast/Field.lean new file mode 100644 index 00000000..9cf1b802 --- /dev/null +++ b/CompPoly/Fields/Goldilocks/Fast/Field.lean @@ -0,0 +1,71 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Varun Thakore +-/ + +import CompPoly.Fields.Goldilocks.Fast.Theorems +import Mathlib.Algebra.Field.TransferInstance + +/-! +# Field Structure for Fast Goldilocks Arithmetic + +This module transfers the canonical Goldilocks field structure to the fast +`UInt64` representation. +-/ + +namespace Goldilocks +namespace Fast + +/-- Ring equivalence between the fast representation and canonical Goldilocks. -/ +def ringEquiv : Field ≃+* Goldilocks.Basic.Field where + toFun := toField + invFun := ofField + left_inv := ofField_toField + right_inv := toField_ofField + map_add' := toField_add + map_mul' := toField_mul + +/-- Applying `ringEquiv` is the same as interpreting a fast value canonically. -/ +@[simp] +theorem ringEquiv_apply (x : Field) : ringEquiv x = toField x := rfl + +/-- Applying the inverse `ringEquiv` converts a canonical value into fast form. -/ +@[simp] +theorem ringEquiv_symm_apply (x : Goldilocks.Basic.Field) : ringEquiv.symm x = ofField x := rfl + +/-- Field instance transferred from canonical Goldilocks through `toField`. -/ +instance (priority := low) instField : _root_.Field Field := + toField_injective.field toField + toField_zero + toField_one + toField_add + toField_mul + toField_neg + toField_sub + toField_inv + toField_div + toField_nsmul + toField_zsmul + toField_nnqsmul + toField_qsmul + toField_npow + toField_zpow + toField_natCast + toField_intCast + toField_nnratCast + toField_ratCast + +/-- Commutative-ring instance inherited from the transferred field structure. -/ +instance (priority := low) instCommRing : CommRing Field := by + infer_instance + +/-- Fast Goldilocks is a non-binary field. -/ +instance (priority := low) instNonBinaryField : NonBinaryField Field where + char_neq_2 := by + intro h + have hv := congrArg Subtype.val h + exact (by decide : (2 : UInt64) ≠ 0) hv + +end Fast +end Goldilocks diff --git a/CompPoly/Fields/Goldilocks/Fast/Internal.lean b/CompPoly/Fields/Goldilocks/Fast/Internal.lean new file mode 100644 index 00000000..ca197fe8 --- /dev/null +++ b/CompPoly/Fields/Goldilocks/Fast/Internal.lean @@ -0,0 +1,488 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Varun Thakore +-/ + +import CompPoly.Fields.Goldilocks.Basic + +/-! +# Internal Lemmas for Fast Goldilocks Arithmetic + +This module contains low-level constants and `UInt64` facts used by the native-word +Goldilocks implementation. +-/ + +namespace Goldilocks +namespace Fast + +/-- Goldilocks modulus `2^64 - 2^32 + 1` as a native word. -/ +@[inline] +def modulus : UInt64 := 0xffffffff00000001 + +/-- Two's complement of the modulus: `2^64 - modulus = 2^32 - 1 = 0xFFFFFFFF`. -/ +@[inline] +def neg_modulus : UInt64 := 0xffffffff + +/-- The native `UInt64` modulus agrees with the mathematical Goldilocks modulus. -/ +@[simp] +theorem modulus_toNat : modulus.toNat = Goldilocks.Basic.fieldSize := by + decide + +/-- The native negated-modulus constant agrees with `2^32 - 1`. -/ +@[simp] +theorem neg_modulus_toNat : neg_modulus.toNat = 2 ^ 32 - 1 := by + decide + +/-- The Goldilocks modulus is positive. -/ +theorem fieldSize_pos : 0 < Goldilocks.Basic.fieldSize := by + decide + +/-- The Goldilocks modulus fits in a `UInt64`. -/ +theorem fieldSize_lt_uint64Size : Goldilocks.Basic.fieldSize < UInt64.size := by + decide + +/-- Every `UInt64` value is below twice the Goldilocks modulus. -/ +theorem uint64_toNat_lt_two_fieldSize (x : UInt64) : + x.toNat < 2 * Goldilocks.Basic.fieldSize := by + exact Nat.lt_trans (UInt64.toNat_lt_size x) (by decide) + +/-- The folding congruence used by the Goldilocks reducer. + +`2^64 ≡ 2^32 - 1 (mod p)`. +-/ +theorem uint64_cast_eq_neg_modulus : + (UInt64.size : Goldilocks.Basic.Field) = (neg_modulus.toNat : Goldilocks.Basic.Field) := by + decide + +/-- Multiplying by `2^32 - 1` after shifting by `2^32` is negation modulo Goldilocks. -/ +theorem pow32_mul_neg_modulus_cast : + ((2 ^ 32 : Nat) : Goldilocks.Basic.Field) * + (neg_modulus.toNat : Goldilocks.Basic.Field) = + -1 := by + decide + +/-- Right shifting a `UInt64` by 32 gives division by `2^32` on naturals. -/ +theorem shiftRight32_toNat (x : UInt64) : + (x >>> 32).toNat = x.toNat / 2 ^ 32 := by + rw [UInt64.toNat_shiftRight] + have h : (32 : UInt64).toNat % 64 = 32 := by + decide + rw [h, Nat.shiftRight_eq_div_pow] + +/-- Masking with `2^32 - 1` gives the low 32 bits on naturals. -/ +theorem and_neg_modulus_toNat (x : UInt64) : + (x &&& neg_modulus).toNat = x.toNat % 2 ^ 32 := by + rw [← UInt64.toNat_toBitVec (x &&& neg_modulus)] + rw [UInt64.toBitVec_and] + rw [BitVec.toNat_and] + rw [UInt64.toNat_toBitVec, UInt64.toNat_toBitVec] + rw [neg_modulus_toNat] + rw [Nat.and_two_pow_sub_one_eq_mod] + +/-- A `UInt64` subtraction with the Goldilocks borrow correction represents subtraction modulo `p`. + +The assumption says the subtrahend is a 32-bit limb, which is the case for +`hi >>> 32` in the 128-bit reducer. +-/ +theorem subBorrow_cast (a b : UInt64) (hb : b.toNat < 2 ^ 32) : + (((if a < b then a - b - neg_modulus else a - b).toNat) : Goldilocks.Basic.Field) = + (a.toNat : Goldilocks.Basic.Field) - (b.toNat : Goldilocks.Basic.Field) := by + by_cases h : a < b + · rw [if_pos h] + have hlt : a.toNat < b.toNat := by + simpa [UInt64.lt_iff_toNat_lt] using h + have hb_le_size : b.toNat ≤ UInt64.size := Nat.le_of_lt (UInt64.toNat_lt_size b) + have hsub_lt_size : UInt64.size - b.toNat + a.toNat < 2 ^ 64 := by + have hsize : UInt64.size = 2 ^ 64 := rfl + rw [hsize] at hb_le_size ⊢ + omega + have hsub_raw : (a - b).toNat = UInt64.size - b.toNat + a.toNat := by + rw [UInt64.toNat_sub] + exact Nat.mod_eq_of_lt hsub_lt_size + have hneg_le : neg_modulus ≤ a - b := by + rw [UInt64.le_iff_toNat_le] + rw [hsub_raw, neg_modulus_toNat] + have hsize : UInt64.size = 2 ^ 64 := rfl + rw [hsize] + omega + rw [UInt64.toNat_sub_of_le _ _ hneg_le] + rw [hsub_raw] + have hneg_le_nat : neg_modulus.toNat ≤ UInt64.size - b.toNat + a.toNat := by + rw [← hsub_raw] + rwa [UInt64.le_iff_toNat_le] at hneg_le + rw [Nat.cast_sub hneg_le_nat] + rw [Nat.cast_add] + rw [Nat.cast_sub hb_le_size] + rw [uint64_cast_eq_neg_modulus] + ring + · rw [if_neg h] + have hle : b ≤ a := by + rw [UInt64.le_iff_toNat_le] + rw [UInt64.lt_iff_toNat_lt] at h + exact Nat.le_of_not_gt h + have hle_nat : b.toNat ≤ a.toNat := by + rwa [UInt64.le_iff_toNat_le] at hle + rw [UInt64.toNat_sub_of_le _ _ hle] + rw [Nat.cast_sub hle_nat] + +/-- A bounded `UInt64` addition with the Goldilocks overflow correction represents +addition modulo `p`. -/ +theorem addOverflowBounded_cast (a b : UInt64) + (hbound : a.toNat + b.toNat < 2 * UInt64.size - neg_modulus.toNat) : + (((if a + b < a then a + b + neg_modulus else a + b).toNat) : + Goldilocks.Basic.Field) = + (a.toNat : Goldilocks.Basic.Field) + (b.toNat : Goldilocks.Basic.Field) := by + by_cases hsum : a.toNat + b.toNat < UInt64.size + · have hnot : ¬a + b < a := by + intro hlt + have hlt_nat : (a + b).toNat < a.toNat := by + simpa [UInt64.lt_iff_toNat_lt] using hlt + rw [UInt64.toNat_add, Nat.mod_eq_of_lt hsum] at hlt_nat + omega + rw [if_neg hnot] + rw [UInt64.toNat_add, Nat.mod_eq_of_lt hsum, Nat.cast_add] + · have hsize_le : UInt64.size ≤ a.toNat + b.toNat := Nat.le_of_not_gt hsum + have hlt : a + b < a := by + rw [UInt64.lt_iff_toNat_lt] + rw [UInt64.toNat_add] + rw [Nat.mod_eq_sub_mod hsize_le] + have hdiff_lt : a.toNat + b.toNat - UInt64.size < UInt64.size := by + have ha := UInt64.toNat_lt_size a + have hb := UInt64.toNat_lt_size b + omega + rw [Nat.mod_eq_of_lt hdiff_lt] + have hb := UInt64.toNat_lt_size b + omega + rw [if_pos hlt] + have hsum_mod : (a + b).toNat = a.toNat + b.toNat - UInt64.size := by + rw [UInt64.toNat_add] + rw [Nat.mod_eq_sub_mod hsize_le] + have hdiff_lt : a.toNat + b.toNat - UInt64.size < UInt64.size := by + have ha := UInt64.toNat_lt_size a + have hb := UInt64.toNat_lt_size b + omega + rw [Nat.mod_eq_of_lt hdiff_lt] + rw [UInt64.toNat_add] + rw [hsum_mod] + have hno_second : a.toNat + b.toNat - UInt64.size + neg_modulus.toNat < UInt64.size := by + omega + rw [Nat.mod_eq_of_lt hno_second] + rw [Nat.cast_add] + rw [Nat.cast_sub hsize_le] + rw [uint64_cast_eq_neg_modulus] + rw [Nat.cast_add] + ring + +/-- Multiplication by `2^32 - 1` does not overflow for a 32-bit limb. -/ +theorem mul_neg_modulus_toNat_of_lt (x : UInt64) (hx : x.toNat < 2 ^ 32) : + (x * neg_modulus).toNat = x.toNat * neg_modulus.toNat := by + rw [UInt64.toNat_mul] + rw [Nat.mod_eq_of_lt] + rw [neg_modulus_toNat] + omega + +/-- The product of a 32-bit limb by `2^32 - 1` leaves enough headroom for correction. -/ +theorem mul_neg_modulus_toNat_le (x : UInt64) (hx : x.toNat < 2 ^ 32) : + (x * neg_modulus).toNat ≤ UInt64.size - 2 * neg_modulus.toNat := by + rw [mul_neg_modulus_toNat_of_lt x hx] + rw [neg_modulus_toNat] + have hx_le : x.toNat ≤ 2 ^ 32 - 1 := by + omega + have hmul : x.toNat * (2 ^ 32 - 1) ≤ (2 ^ 32 - 1) * (2 ^ 32 - 1) := + Nat.mul_le_mul_right _ hx_le + have hconst : (2 ^ 32 - 1) * (2 ^ 32 - 1) ≤ UInt64.size - 2 * (2 ^ 32 - 1) := by + decide + exact Nat.le_trans hmul hconst + +/-- Splitting a high word into 32-bit limbs matches the Goldilocks folding congruence. -/ +theorem hi_split_cast (hi : UInt64) : + (hi.toNat : Goldilocks.Basic.Field) * (UInt64.size : Goldilocks.Basic.Field) = + ((hi &&& neg_modulus).toNat : Goldilocks.Basic.Field) * + (neg_modulus.toNat : Goldilocks.Basic.Field) - + ((hi >>> 32).toNat : Goldilocks.Basic.Field) := by + have hsplit_nat : hi.toNat = hi.toNat % 2 ^ 32 + 2 ^ 32 * (hi.toNat / 2 ^ 32) := by + rw [Nat.mod_add_div] + have hcast_split : + (hi.toNat : Goldilocks.Basic.Field) = + ((hi.toNat % 2 ^ 32 : Nat) : Goldilocks.Basic.Field) + + (((2 ^ 32 : Nat) : Goldilocks.Basic.Field) * + ((hi.toNat / 2 ^ 32 : Nat) : Goldilocks.Basic.Field)) := by + simpa [Nat.cast_add, Nat.cast_mul] using + congrArg (fun n : Nat => (n : Goldilocks.Basic.Field)) hsplit_nat + rw [hcast_split] + rw [and_neg_modulus_toNat, shiftRight32_toNat, uint64_cast_eq_neg_modulus] + rw [add_mul] + conv_lhs => + enter [2] + rw [mul_assoc] + rw [mul_comm ((hi.toNat / 2 ^ 32 : Nat) : Goldilocks.Basic.Field)] + rw [← mul_assoc] + rw [pow32_mul_neg_modulus_cast] + ring + +/-- A `UInt64` value decomposes into its low and high 32-bit limbs. -/ +theorem uint64_split32 (x : UInt64) : + x.toNat = (x &&& neg_modulus).toNat + 2 ^ 32 * (x >>> 32).toNat := by + rw [and_neg_modulus_toNat, shiftRight32_toNat] + rw [Nat.mod_add_div] + +/-- The low 32-bit limb of a `UInt64` is below `2^32`. -/ +theorem uint64_low32_lt (x : UInt64) : + (x &&& neg_modulus).toNat < 2 ^ 32 := by + rw [and_neg_modulus_toNat] + exact Nat.mod_lt _ (by decide) + +/-- The high 32-bit limb of a `UInt64` is below `2^32`. -/ +theorem uint64_high32_lt (x : UInt64) : + (x >>> 32).toNat < 2 ^ 32 := by + rw [shiftRight32_toNat] + have hx := UInt64.toNat_lt_size x + change x.toNat < 2 ^ 64 at hx + exact Nat.div_lt_of_lt_mul hx + +/-- Multiplying two 32-bit limbs does not overflow `UInt64`. -/ +theorem mul32_toNat (a b : UInt64) (ha : a.toNat < 2 ^ 32) (hb : b.toNat < 2 ^ 32) : + (a * b).toNat = a.toNat * b.toNat := by + rw [UInt64.toNat_mul] + rw [Nat.mod_eq_of_lt] + nlinarith + +/-- Algebraic decomposition of a product after splitting both factors into 32-bit limbs. -/ +theorem product_split32 (x y : UInt64) : + x.toNat * y.toNat = + (x &&& neg_modulus).toNat * (y &&& neg_modulus).toNat + + 2 ^ 32 * + ((x &&& neg_modulus).toNat * (y >>> 32).toNat + + (x >>> 32).toNat * (y &&& neg_modulus).toNat) + + 2 ^ 64 * ((x >>> 32).toNat * (y >>> 32).toNat) := by + rw [uint64_split32 x, uint64_split32 y] + ring_nf + +/-- Low word returned by a 64-by-64 product implementation. -/ +theorem wideMul_low_toNat (lo : UInt64) (x y : UInt64) (hlo : lo = x * y) : + lo.toNat = x.toNat * y.toNat % UInt64.size := by + rw [hlo, UInt64.toNat_mul] + +/-- Pure Nat carry formula for the high word of a 32-bit-limb 64-by-64 product. -/ +theorem wideMul_high_nat + (p00 p01 p10 p11 : Nat) + (_hp00 : p00 < 2 ^ 64) + (_hp01 : p01 < 2 ^ 64) + (_hp10 : p10 < 2 ^ 64) + (_hp11 : p11 < 2 ^ 64) : + let B := 2 ^ 32 + let carry := p00 / B + p01 % B + p10 % B + p11 + p01 / B + p10 / B + carry / B = + (p00 + B * (p01 + p10) + B ^ 2 * p11) / B ^ 2 := by + dsimp + let carry := p00 / 4294967296 + p01 % 4294967296 + p10 % 4294967296 + let q := p11 + p01 / 4294967296 + p10 / 4294967296 + carry / 4294967296 + have hN : + p00 + 4294967296 * (p01 + p10) + 18446744073709551616 * p11 = + p00 % 4294967296 + 4294967296 * (carry % 4294967296) + + 18446744073709551616 * q := by + have hp00d : p00 % 4294967296 + 4294967296 * (p00 / 4294967296) = p00 := + Nat.mod_add_div p00 4294967296 + have hp01d : p01 % 4294967296 + 4294967296 * (p01 / 4294967296) = p01 := + Nat.mod_add_div p01 4294967296 + have hp10d : p10 % 4294967296 + 4294967296 * (p10 / 4294967296) = p10 := + Nat.mod_add_div p10 4294967296 + have hcd : carry % 4294967296 + 4294967296 * (carry / 4294967296) = carry := + Nat.mod_add_div carry 4294967296 + subst q + subst carry + omega + rw [hN] + change q = + (p00 % 4294967296 + 4294967296 * (carry % 4294967296) + + 18446744073709551616 * q) / + 18446744073709551616 + rw [Nat.mul_comm 18446744073709551616 q] + rw [Nat.add_mul_div_right _ _ (show 0 < 18446744073709551616 by decide)] + rw [Nat.div_eq_of_lt] + · rw [Nat.zero_add] + · subst carry + have hmod0 : p00 % 4294967296 < 4294967296 := Nat.mod_lt _ (by decide) + have hmod1 : + (p00 / 4294967296 + p01 % 4294967296 + p10 % 4294967296) % + 4294967296 < + 4294967296 := Nat.mod_lt _ (by decide) + omega + +/-- High word returned by the 32-bit-limb `UInt64` multiplication algorithm. -/ +theorem wideMul_high_toNat + (x y hi : UInt64) + (hhi : + hi = + let xLo := x &&& neg_modulus + let xHi := x >>> 32 + let yLo := y &&& neg_modulus + let yHi := y >>> 32 + let p00 := xLo * yLo + let p01 := xLo * yHi + let p10 := xHi * yLo + let p11 := xHi * yHi + let carry := (p00 >>> 32) + (p01 &&& neg_modulus) + (p10 &&& neg_modulus) + p11 + (p01 >>> 32) + (p10 >>> 32) + (carry >>> 32)) : + hi.toNat = x.toNat * y.toNat / UInt64.size := by + let xLo := x &&& neg_modulus + let xHi := x >>> 32 + let yLo := y &&& neg_modulus + let yHi := y >>> 32 + let p00 := xLo * yLo + let p01 := xLo * yHi + let p10 := xHi * yLo + let p11 := xHi * yHi + let carry := (p00 >>> 32) + (p01 &&& neg_modulus) + (p10 &&& neg_modulus) + have hxLo_lt : xLo.toNat < 2 ^ 32 := by + subst xLo + exact uint64_low32_lt x + have hxHi_lt : xHi.toNat < 2 ^ 32 := by + subst xHi + exact uint64_high32_lt x + have hyLo_lt : yLo.toNat < 2 ^ 32 := by + subst yLo + exact uint64_low32_lt y + have hyHi_lt : yHi.toNat < 2 ^ 32 := by + subst yHi + exact uint64_high32_lt y + have hp00_nat : p00.toNat = xLo.toNat * yLo.toNat := by + subst p00 + exact mul32_toNat xLo yLo hxLo_lt hyLo_lt + have hp01_nat : p01.toNat = xLo.toNat * yHi.toNat := by + subst p01 + exact mul32_toNat xLo yHi hxLo_lt hyHi_lt + have hp10_nat : p10.toNat = xHi.toNat * yLo.toNat := by + subst p10 + exact mul32_toNat xHi yLo hxHi_lt hyLo_lt + have hp11_nat : p11.toNat = xHi.toNat * yHi.toNat := by + subst p11 + exact mul32_toNat xHi yHi hxHi_lt hyHi_lt + have hp00_lt : p00.toNat < 2 ^ 64 := by + rw [hp00_nat] + nlinarith [hxLo_lt, hyLo_lt] + have hp01_lt : p01.toNat < 2 ^ 64 := by + rw [hp01_nat] + nlinarith [hxLo_lt, hyHi_lt] + have hp10_lt : p10.toNat < 2 ^ 64 := by + rw [hp10_nat] + nlinarith [hxHi_lt, hyLo_lt] + have hp11_lt : p11.toNat < 2 ^ 64 := by + rw [hp11_nat] + nlinarith [hxHi_lt, hyHi_lt] + have hcarry_nat : + carry.toNat = p00.toNat / 2 ^ 32 + p01.toNat % 2 ^ 32 + p10.toNat % 2 ^ 32 := by + subst carry + rw [UInt64.toNat_add, UInt64.toNat_add] + rw [shiftRight32_toNat, and_neg_modulus_toNat, and_neg_modulus_toNat] + have hp00_hi_lt : p00.toNat / 2 ^ 32 < 2 ^ 32 := by + rw [Nat.div_lt_iff_lt_mul (by decide : 0 < 2 ^ 32)] + simpa [pow_add] using hp00_lt + have hp01_lo_lt : p01.toNat % 2 ^ 32 < 2 ^ 32 := Nat.mod_lt _ (by decide) + have hp10_lo_lt : p10.toNat % 2 ^ 32 < 2 ^ 32 := Nat.mod_lt _ (by decide) + have hsum01 : + p00.toNat / 2 ^ 32 + p01.toNat % 2 ^ 32 < 2 ^ 64 := by + omega + have hsum012 : + p00.toNat / 2 ^ 32 + p01.toNat % 2 ^ 32 + p10.toNat % 2 ^ 32 < + 2 ^ 64 := by + omega + rw [Nat.mod_eq_of_lt hsum01, Nat.mod_eq_of_lt hsum012] + have hwide := + wideMul_high_nat p00.toNat p01.toNat p10.toNat p11.toNat hp00_lt hp01_lt hp10_lt hp11_lt + have hwide' : + p11.toNat + p01.toNat / 2 ^ 32 + p10.toNat / 2 ^ 32 + + (p00.toNat / 2 ^ 32 + p01.toNat % 2 ^ 32 + p10.toNat % 2 ^ 32) / 2 ^ 32 = + (p00.toNat + 2 ^ 32 * (p01.toNat + p10.toNat) + (2 ^ 32) ^ 2 * p11.toNat) / + (2 ^ 32) ^ 2 := by + simpa using hwide + have hquot_lt : + p11.toNat + p01.toNat / 2 ^ 32 + p10.toNat / 2 ^ 32 + carry.toNat / 2 ^ 32 < + UInt64.size := by + rw [hcarry_nat] + rw [hwide'] + have hprod_bound : x.toNat * y.toNat < UInt64.size * UInt64.size := by + exact + mul_lt_mul'' (UInt64.toNat_lt_size x) (UInt64.toNat_lt_size y) (Nat.zero_le _) + (Nat.zero_le _) + have hsplit : + p00.toNat + 2 ^ 32 * (p01.toNat + p10.toNat) + (2 ^ 32) ^ 2 * p11.toNat = + x.toNat * y.toNat := by + rw [hp00_nat, hp01_nat, hp10_nat, hp11_nat] + subst p00 + subst p01 + subst p10 + subst p11 + subst xLo + subst xHi + subst yLo + subst yHi + simpa [pow_add, pow_mul] using (product_split32 x y).symm + rw [hsplit] + change x.toNat * y.toNat / UInt64.size < UInt64.size + rw [Nat.div_lt_iff_lt_mul (by decide : 0 < UInt64.size)] + exact hprod_bound + have hhi_nat : + hi.toNat = p11.toNat + p01.toNat / 2 ^ 32 + p10.toNat / 2 ^ 32 + carry.toNat / 2 ^ 32 := by + rw [hhi] + dsimp only + change (p11 + (p01 >>> 32) + (p10 >>> 32) + (carry >>> 32)).toNat = + p11.toNat + p01.toNat / 2 ^ 32 + p10.toNat / 2 ^ 32 + carry.toNat / 2 ^ 32 + rw [UInt64.toNat_add, UInt64.toNat_add, UInt64.toNat_add] + rw [shiftRight32_toNat, shiftRight32_toNat, shiftRight32_toNat] + have hsum01 : p11.toNat + p01.toNat / 2 ^ 32 < UInt64.size := by + have hle : p11.toNat + p01.toNat / 2 ^ 32 ≤ + p11.toNat + p01.toNat / 2 ^ 32 + p10.toNat / 2 ^ 32 + carry.toNat / 2 ^ 32 := by + omega + exact lt_of_le_of_lt hle hquot_lt + have hsum012 : + p11.toNat + p01.toNat / 2 ^ 32 + p10.toNat / 2 ^ 32 < UInt64.size := by + have hle : p11.toNat + p01.toNat / 2 ^ 32 + p10.toNat / 2 ^ 32 ≤ + p11.toNat + p01.toNat / 2 ^ 32 + p10.toNat / 2 ^ 32 + carry.toNat / 2 ^ 32 := by + omega + exact lt_of_le_of_lt hle hquot_lt + have hquot_lt_pow : + p11.toNat + p01.toNat / 2 ^ 32 + p10.toNat / 2 ^ 32 + carry.toNat / 2 ^ 32 < + 2 ^ 64 := by + simpa [UInt64.size] using hquot_lt + have hsum01_pow : p11.toNat + p01.toNat / 2 ^ 32 < 2 ^ 64 := by + simpa [UInt64.size] using hsum01 + have hsum012_pow : + p11.toNat + p01.toNat / 2 ^ 32 + p10.toNat / 2 ^ 32 < 2 ^ 64 := by + simpa [UInt64.size] using hsum012 + rw [Nat.mod_eq_of_lt hsum01_pow, Nat.mod_eq_of_lt hsum012_pow, + Nat.mod_eq_of_lt hquot_lt_pow] + rw [hhi_nat] + rw [hcarry_nat] + rw [hwide'] + have hsplit : + p00.toNat + 2 ^ 32 * (p01.toNat + p10.toNat) + (2 ^ 32) ^ 2 * p11.toNat = + x.toNat * y.toNat := by + rw [hp00_nat, hp01_nat, hp10_nat, hp11_nat] + subst p00 + subst p01 + subst p10 + subst p11 + subst xLo + subst xHi + subst yLo + subst yHi + simpa [pow_add, pow_mul] using (product_split32 x y).symm + rw [hsplit] + rfl + +/-- Combined semantic correctness of a 64-by-64 product represented by low and high words. -/ +theorem wideMul_cast + (x y lo hi : UInt64) + (hlo : lo = x * y) + (hhi : hi.toNat = x.toNat * y.toNat / UInt64.size) : + (lo.toNat : Goldilocks.Basic.Field) + + (hi.toNat : Goldilocks.Basic.Field) * (UInt64.size : Goldilocks.Basic.Field) = + (x.toNat : Goldilocks.Basic.Field) * (y.toNat : Goldilocks.Basic.Field) := by + rw [wideMul_low_toNat lo x y hlo, hhi] + rw [← Nat.cast_mul, ← Nat.cast_add, Nat.mul_comm (x.toNat * y.toNat / UInt64.size), + Nat.mod_add_div, Nat.cast_mul] + +end Fast +end Goldilocks diff --git a/CompPoly/Fields/Goldilocks/Fast/Reduction.lean b/CompPoly/Fields/Goldilocks/Fast/Reduction.lean new file mode 100644 index 00000000..95d0448d --- /dev/null +++ b/CompPoly/Fields/Goldilocks/Fast/Reduction.lean @@ -0,0 +1,478 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Varun Thakore +-/ + +import CompPoly.Fields.Goldilocks.Fast.Internal + +/-! +# Raw Reducers for Fast Goldilocks Arithmetic + +This module contains the raw `UInt64` arithmetic kernels and their correctness +lemmas. The public fast field API wraps these kernels in +`CompPoly.Fields.Goldilocks.Fast`. +-/ + +namespace Goldilocks +namespace Fast +namespace Reduction + +/-- Full 64-by-64 product as `(lo, hi)` words, computed from 32-bit limbs. -/ +@[inline] +def wideMul (x y : UInt64) : UInt64 × UInt64 := + let xLo := x &&& neg_modulus + let xHi := x >>> 32 + let yLo := y &&& neg_modulus + let yHi := y >>> 32 + let p00 := xLo * yLo + let p01 := xLo * yHi + let p10 := xHi * yLo + let p11 := xHi * yHi + let carry := (p00 >>> 32) + (p01 &&& neg_modulus) + (p10 &&& neg_modulus) + let hi := p11 + (p01 >>> 32) + (p10 >>> 32) + (carry >>> 32) + (x * y, hi) + +/-- Raw one-word reduction for a `UInt64` value. + +Since every `UInt64` is below `2^64 = p + 2^32 - 1`, one subtraction by `p` +is enough to canonicalize a native word. +-/ +@[inline] +def reduceUInt64Raw (x : UInt64) : UInt64 := + if x < modulus then x else x - modulus + +/-- The raw one-word reducer returns a canonical representative. -/ +theorem reduceUInt64Raw_lt (x : UInt64) : + (reduceUInt64Raw x).toNat < Goldilocks.Basic.fieldSize := by + unfold reduceUInt64Raw + by_cases hx : x < modulus + · rw [if_pos hx] + rw [UInt64.lt_iff_toNat_lt, modulus_toNat] at hx + exact hx + · rw [if_neg hx] + have hmod_le_x_nat : Goldilocks.Basic.fieldSize ≤ x.toNat := by + rw [UInt64.lt_iff_toNat_lt, modulus_toNat] at hx + exact Nat.le_of_not_gt hx + have hmod_le_x : modulus ≤ x := by + rw [UInt64.le_iff_toNat_le, modulus_toNat] + exact hmod_le_x_nat + rw [UInt64.toNat_sub_of_le _ _ hmod_le_x, modulus_toNat] + have hx_lt_two := uint64_toNat_lt_two_fieldSize x + omega + +/-- One-word reduction preserves the represented canonical field element. -/ +theorem reduceUInt64Raw_cast (x : UInt64) : + ((reduceUInt64Raw x).toNat : Goldilocks.Basic.Field) = + (x.toNat : Goldilocks.Basic.Field) := by + unfold reduceUInt64Raw + by_cases hx : x < modulus + · rw [if_pos hx] + · rw [if_neg hx] + have hmod_le_x : modulus ≤ x := by + rw [UInt64.le_iff_toNat_le, modulus_toNat] + rw [UInt64.lt_iff_toNat_lt, modulus_toNat] at hx + exact Nat.le_of_not_gt hx + rw [UInt64.toNat_sub_of_le _ _ hmod_le_x, modulus_toNat] + rw [Nat.cast_sub (by + rw [UInt64.le_iff_toNat_le, modulus_toNat] at hmod_le_x + exact hmod_le_x)] + simp + +/-- Raw reduction of a 128-bit integer represented by low and high words modulo Goldilocks. -/ +@[inline] +def reduceUInt128Raw (lo hi : UInt64) : UInt64 := + let hi_hi := hi >>> 32 + let hi_lo := hi &&& neg_modulus + + let borrow := lo < hi_hi + let t0 := lo - hi_hi + let t0 := if borrow then t0 - neg_modulus else t0 + + let t1 := hi_lo * neg_modulus + + let t2 := t0 + t1 + let overflow := t2 < t0 + let t2 := if overflow then t2 + neg_modulus else t2 + + reduceUInt64Raw t2 + +/-- The raw 128-bit reducer returns a canonical representative below the modulus. -/ +theorem reduceUInt128Raw_lt (lo hi : UInt64) : + (reduceUInt128Raw lo hi).toNat < Goldilocks.Basic.fieldSize := by + unfold reduceUInt128Raw + apply reduceUInt64Raw_lt + +/-- Semantic correctness of raw 128-bit Goldilocks reduction. -/ +theorem reduceUInt128Raw_cast (lo hi : UInt64) : + ((reduceUInt128Raw lo hi).toNat : Goldilocks.Basic.Field) = + (lo.toNat : Goldilocks.Basic.Field) + + (hi.toNat : Goldilocks.Basic.Field) * (UInt64.size : Goldilocks.Basic.Field) := by + let hi_hi := hi >>> 32 + let hi_lo := hi &&& neg_modulus + let t0 := if lo < hi_hi then lo - hi_hi - neg_modulus else lo - hi_hi + let t1 := hi_lo * neg_modulus + let t2 := if t0 + t1 < t0 then t0 + t1 + neg_modulus else t0 + t1 + change ((reduceUInt64Raw t2).toNat : Goldilocks.Basic.Field) = + (lo.toNat : Goldilocks.Basic.Field) + + (hi.toNat : Goldilocks.Basic.Field) * (UInt64.size : Goldilocks.Basic.Field) + have hred := reduceUInt64Raw_cast t2 + change ((reduceUInt64Raw t2).toNat : Goldilocks.Basic.Field) = + (t2.toNat : Goldilocks.Basic.Field) at hred + rw [hred] + have hhi_hi_lt : hi_hi.toNat < 2 ^ 32 := by + rw [show hi_hi = hi >>> 32 by rfl, shiftRight32_toNat] + have hhi := UInt64.toNat_lt_size hi + change hi.toNat < 2 ^ 64 at hhi + omega + have hhi_lo_lt : hi_lo.toNat < 2 ^ 32 := by + rw [show hi_lo = hi &&& neg_modulus by rfl, and_neg_modulus_toNat] + exact Nat.mod_lt _ (by decide) + have ht0_cast : + (t0.toNat : Goldilocks.Basic.Field) = + (lo.toNat : Goldilocks.Basic.Field) - (hi_hi.toNat : Goldilocks.Basic.Field) := by + rw [show t0 = if lo < hi_hi then lo - hi_hi - neg_modulus else lo - hi_hi by rfl] + exact subBorrow_cast lo hi_hi hhi_hi_lt + have ht1_cast : + (t1.toNat : Goldilocks.Basic.Field) = + (hi_lo.toNat : Goldilocks.Basic.Field) * + (neg_modulus.toNat : Goldilocks.Basic.Field) := by + rw [show t1 = hi_lo * neg_modulus by rfl] + rw [mul_neg_modulus_toNat_of_lt hi_lo hhi_lo_lt] + rw [Nat.cast_mul] + have ht2_bound : t0.toNat + t1.toNat < 2 * UInt64.size - neg_modulus.toNat := by + have ht0_lt := UInt64.toNat_lt_size t0 + have ht1_le : t1.toNat ≤ UInt64.size - 2 * neg_modulus.toNat := by + simpa [t1] using mul_neg_modulus_toNat_le hi_lo hhi_lo_lt + have htwice_neg_le_size : 2 * neg_modulus.toNat ≤ UInt64.size := by + decide + omega + have ht2_cast : + (t2.toNat : Goldilocks.Basic.Field) = + (t0.toNat : Goldilocks.Basic.Field) + (t1.toNat : Goldilocks.Basic.Field) := by + rw [show t2 = if t0 + t1 < t0 then t0 + t1 + neg_modulus else t0 + t1 by rfl] + exact addOverflowBounded_cast t0 t1 ht2_bound + rw [ht2_cast, ht0_cast, ht1_cast] + rw [hi_split_cast hi] + ring + +/-- Raw reduction of a 64-by-64 product modulo Goldilocks. -/ +@[inline] +def reduceMulRaw (x y : UInt64) : UInt64 := + let product := wideMul x y + reduceUInt128Raw product.1 product.2 + +/-- Product reduction returns a canonical representative below the modulus. -/ +theorem reduceMulRaw_lt (x y : UInt64) : + (reduceMulRaw x y).toNat < Goldilocks.Basic.fieldSize := by + unfold reduceMulRaw + apply reduceUInt128Raw_lt + +/-- Semantic correctness of native 64-by-64 product reduction. -/ +theorem reduceMulRaw_cast (x y : UInt64) : + ((reduceMulRaw x y).toNat : Goldilocks.Basic.Field) = + (x.toNat : Goldilocks.Basic.Field) * (y.toNat : Goldilocks.Basic.Field) := by + unfold reduceMulRaw + rw [reduceUInt128Raw_cast] + exact + wideMul_cast x y (wideMul x y).1 (wideMul x y).2 + (by unfold wideMul; rfl) + (by + unfold wideMul + exact wideMul_high_toNat x y _ rfl) + +/-- Raw one-step reduction for a 65-bit addition represented by low word and carry. -/ +@[inline] +def reduceAddWithCarryRaw (lo : UInt64) (carry : Bool) : UInt64 := + if carry then + lo + neg_modulus + else + reduceUInt64Raw lo + +/-- Addition reduction returns a canonical representative. -/ +theorem reduceAddWithCarryRaw_lt (lo : UInt64) (carry : Bool) + (h : + lo.toNat + (if carry then UInt64.size else 0) < 2 * Goldilocks.Basic.fieldSize) : + (reduceAddWithCarryRaw lo carry).toNat < Goldilocks.Basic.fieldSize := by + unfold reduceAddWithCarryRaw + cases carry + · simp only [Bool.false_eq_true, if_false] + exact reduceUInt64Raw_lt lo + · simp only [↓reduceIte] at h ⊢ + rw [UInt64.toNat_add] + have hsum_lt_field : lo.toNat + neg_modulus.toNat < Goldilocks.Basic.fieldSize := by + rw [neg_modulus_toNat, Goldilocks.Basic.fieldSize] + rw [Goldilocks.Basic.fieldSize, UInt64.size] at h + omega + have hsum_lt_size : lo.toNat + neg_modulus.toNat < UInt64.size := + Nat.lt_trans hsum_lt_field fieldSize_lt_uint64Size + rw [Nat.mod_eq_of_lt hsum_lt_size] + exact hsum_lt_field + +/-- Semantic correctness of addition reduction with carry. -/ +theorem reduceAddWithCarryRaw_cast (lo : UInt64) (carry : Bool) + (h : + lo.toNat + (if carry then UInt64.size else 0) < 2 * Goldilocks.Basic.fieldSize) : + ((reduceAddWithCarryRaw lo carry).toNat : Goldilocks.Basic.Field) = + (lo.toNat : Goldilocks.Basic.Field) + + (if carry then (UInt64.size : Goldilocks.Basic.Field) else 0) := by + unfold reduceAddWithCarryRaw + cases carry + · simp only [Bool.false_eq_true, if_false, add_zero] + exact reduceUInt64Raw_cast lo + · simp only [↓reduceIte] + change lo.toNat + UInt64.size < 2 * Goldilocks.Basic.fieldSize at h + rw [UInt64.toNat_add] + have hsum_lt_field : lo.toNat + neg_modulus.toNat < Goldilocks.Basic.fieldSize := by + rw [neg_modulus_toNat, Goldilocks.Basic.fieldSize] + rw [Goldilocks.Basic.fieldSize, UInt64.size] at h + omega + have hsum_lt_size : lo.toNat + neg_modulus.toNat < UInt64.size := + Nat.lt_trans hsum_lt_field fieldSize_lt_uint64Size + rw [Nat.mod_eq_of_lt hsum_lt_size] + rw [Nat.cast_add] + rw [uint64_cast_eq_neg_modulus] + +/-- The wrapped word and carry produced by adding two canonical representatives is bounded. -/ +theorem addWithCarry_bound (x y : UInt64) + (hx : x.toNat < Goldilocks.Basic.fieldSize) + (hy : y.toNat < Goldilocks.Basic.fieldSize) : + let lo := x + y + let carry := decide (lo < x) + lo.toNat + (if carry then UInt64.size else 0) < 2 * Goldilocks.Basic.fieldSize := by + intro lo carry + by_cases hcarry : carry + · simp only [hcarry] + have hlo_lt_x : lo.toNat < x.toNat := by + simpa [carry, UInt64.lt_iff_toNat_lt] using hcarry + have hsum_ge_size : UInt64.size ≤ x.toNat + y.toNat := by + by_contra hnot + have hsum_lt_size : x.toNat + y.toNat < UInt64.size := + Nat.lt_of_not_ge hnot + have hlo_eq : lo.toNat = x.toNat + y.toNat := by + rw [show lo = x + y by rfl, UInt64.toNat_add] + exact Nat.mod_eq_of_lt hsum_lt_size + omega + have hsum_lt_2size : x.toNat + y.toNat < 2 * UInt64.size := by + nlinarith [UInt64.toNat_lt_size x, UInt64.toNat_lt_size y] + have hlo_eq : lo.toNat = x.toNat + y.toNat - UInt64.size := by + rw [show lo = x + y by rfl, UInt64.toNat_add] + rw [Nat.mod_eq_sub_mod (show x.toNat + y.toNat ≥ UInt64.size by + exact hsum_ge_size)] + rw [Nat.mod_eq_of_lt] + omega + rw [hlo_eq] + change x.toNat + y.toNat - UInt64.size + UInt64.size < + 2 * Goldilocks.Basic.fieldSize + have hsum_lt_field : x.toNat + y.toNat < 2 * Goldilocks.Basic.fieldSize := by + omega + omega + · simp only [hcarry] + exact uint64_toNat_lt_two_fieldSize lo + +/-- The wrapped word and carry produced by native addition reconstruct the exact Nat sum. -/ +theorem addWithCarry_value (x y : UInt64) : + let lo := x + y + let carry := decide (lo < x) + lo.toNat + (if carry then UInt64.size else 0) = x.toNat + y.toNat := by + intro lo carry + by_cases hcarry : carry + · simp only [hcarry] + have hsum_ge_size : UInt64.size ≤ x.toNat + y.toNat := by + by_contra hnot + have hsum_lt_size : x.toNat + y.toNat < UInt64.size := + Nat.lt_of_not_ge hnot + have hlo_eq : lo.toNat = x.toNat + y.toNat := by + rw [show lo = x + y by rfl, UInt64.toNat_add] + exact Nat.mod_eq_of_lt hsum_lt_size + have hlo_lt_x : lo.toNat < x.toNat := by + simpa [carry, UInt64.lt_iff_toNat_lt] using hcarry + omega + have hsum_lt_2size : x.toNat + y.toNat < 2 * UInt64.size := by + nlinarith [UInt64.toNat_lt_size x, UInt64.toNat_lt_size y] + have hlo_eq : lo.toNat = x.toNat + y.toNat - UInt64.size := by + rw [show lo = x + y by rfl, UInt64.toNat_add] + rw [Nat.mod_eq_sub_mod (show x.toNat + y.toNat ≥ UInt64.size by + exact hsum_ge_size)] + rw [Nat.mod_eq_of_lt] + omega + rw [hlo_eq] + change x.toNat + y.toNat - UInt64.size + UInt64.size = + x.toNat + y.toNat + omega + · simp only [hcarry] + have hnot_lo_lt_x : ¬lo.toNat < x.toNat := by + intro hlo_lt_x + apply hcarry + simpa [carry, UInt64.lt_iff_toNat_lt] using hlo_lt_x + have hsum_lt_size : x.toNat + y.toNat < UInt64.size := by + by_contra hnot + have hsum_ge_size : UInt64.size ≤ x.toNat + y.toNat := + Nat.le_of_not_gt hnot + have hsum_lt_2size : x.toNat + y.toNat < 2 * UInt64.size := by + nlinarith [UInt64.toNat_lt_size x, UInt64.toNat_lt_size y] + have hlo_eq : lo.toNat = x.toNat + y.toNat - UInt64.size := by + rw [show lo = x + y by rfl, UInt64.toNat_add] + rw [Nat.mod_eq_sub_mod (show x.toNat + y.toNat ≥ UInt64.size by + exact hsum_ge_size)] + rw [Nat.mod_eq_of_lt] + omega + have hy_lt_size : y.toNat < UInt64.size := UInt64.toNat_lt_size y + omega + rw [show lo = x + y by rfl, UInt64.toNat_add] + exact Nat.mod_eq_of_lt hsum_lt_size + +/-- Raw modular negation in canonical form. -/ +@[inline] +def negRaw (x : UInt64) : UInt64 := + if x = 0 then 0 else modulus - x + +/-- Raw negation returns a canonical representative when given one. -/ +theorem negRaw_lt (x : UInt64) (hx : x.toNat < Goldilocks.Basic.fieldSize) : + (negRaw x).toNat < Goldilocks.Basic.fieldSize := by + unfold negRaw + by_cases hzero : x = 0 + · rw [if_pos hzero] + decide + · rw [if_neg hzero] + have hx_ne_nat : x.toNat ≠ 0 := by + intro hz + apply hzero + apply UInt64.toNat_inj.mp + rw [hz] + decide + have hx_pos : 0 < x.toNat := Nat.pos_of_ne_zero hx_ne_nat + have hx_le_mod : x ≤ modulus := by + rw [UInt64.le_iff_toNat_le, modulus_toNat] + exact Nat.le_of_lt hx + rw [UInt64.toNat_sub_of_le _ _ hx_le_mod, modulus_toNat] + omega + +/-- Raw negation agrees with canonical-field negation. -/ +theorem negRaw_cast (x : UInt64) (hx : x.toNat < Goldilocks.Basic.fieldSize) : + ((negRaw x).toNat : Goldilocks.Basic.Field) = + -((x.toNat : Goldilocks.Basic.Field)) := by + unfold negRaw + by_cases hzero : x = 0 + · rw [if_pos hzero] + have hxNat : x.toNat = 0 := by + simpa using congrArg UInt64.toNat hzero + rw [hxNat] + simp + · rw [if_neg hzero] + have hle : x ≤ modulus := by + rw [UInt64.le_iff_toNat_le, modulus_toNat] + exact Nat.le_of_lt hx + rw [UInt64.toNat_sub_of_le _ _ hle, modulus_toNat] + rw [Nat.cast_sub (by + rw [UInt64.le_iff_toNat_le, modulus_toNat] at hle + exact hle)] + rw [ZMod.natCast_self] + ring + +/-- Raw modular subtraction in canonical form. -/ +@[inline] +def subRaw (x y : UInt64) : UInt64 := + if y ≤ x then x - y else x - y - neg_modulus + +/-- Raw subtraction returns a canonical representative when given canonical operands. -/ +theorem subRaw_lt (x y : UInt64) + (hx : x.toNat < Goldilocks.Basic.fieldSize) + (hy : y.toNat < Goldilocks.Basic.fieldSize) : + (subRaw x y).toNat < Goldilocks.Basic.fieldSize := by + unfold subRaw + by_cases hxy : y ≤ x + · rw [if_pos hxy] + rw [UInt64.toNat_sub_of_le _ _ hxy] + have hy_le_x : y.toNat ≤ x.toNat := by + simpa [UInt64.le_iff_toNat_le] using hxy + omega + · rw [if_neg hxy] + have hx_lt_y : x.toNat < y.toNat := by + have hnot : ¬y.toNat ≤ x.toNat := by + intro hle + apply hxy + rw [UInt64.le_iff_toNat_le] + exact hle + exact Nat.lt_of_not_ge hnot + have hraw_lt_size : 2 ^ 64 - y.toNat + x.toNat < 2 ^ 64 := by + have hx_lt_size : x.toNat < 2 ^ 64 := by + simpa [UInt64.size] using UInt64.toNat_lt_size x + have hy_lt_size : y.toNat < 2 ^ 64 := by + simpa [UInt64.size] using UInt64.toNat_lt_size y + omega + have hraw_toNat : + (x - y).toNat = 2 ^ 64 - y.toNat + x.toNat := by + rw [UInt64.toNat_sub] + exact Nat.mod_eq_of_lt hraw_lt_size + have hneg_le_raw_nat : neg_modulus.toNat ≤ (x - y).toNat := by + rw [hraw_toNat, neg_modulus_toNat] + have hy_lt_field : y.toNat < 2 ^ 64 - 2 ^ 32 + 1 := by + simpa [Goldilocks.Basic.fieldSize] using hy + omega + have hneg_le_raw : neg_modulus ≤ x - y := by + rw [UInt64.le_iff_toNat_le] + exact hneg_le_raw_nat + rw [UInt64.toNat_sub_of_le _ _ hneg_le_raw, hraw_toNat, neg_modulus_toNat] + change 2 ^ 64 - y.toNat + x.toNat - (2 ^ 32 - 1) < + 2 ^ 64 - 2 ^ 32 + 1 + have hx_lt_field : x.toNat < 2 ^ 64 - 2 ^ 32 + 1 := by + simpa [Goldilocks.Basic.fieldSize] using hx + have hy_lt_field : y.toNat < 2 ^ 64 - 2 ^ 32 + 1 := by + simpa [Goldilocks.Basic.fieldSize] using hy + omega + +/-- Raw subtraction agrees with canonical-field subtraction for canonical operands. -/ +theorem subRaw_cast (x y : UInt64) + (_hx : x.toNat < Goldilocks.Basic.fieldSize) + (hy : y.toNat < Goldilocks.Basic.fieldSize) : + ((subRaw x y).toNat : Goldilocks.Basic.Field) = + (x.toNat : Goldilocks.Basic.Field) - (y.toNat : Goldilocks.Basic.Field) := by + unfold subRaw + by_cases hxy : y ≤ x + · rw [if_pos hxy] + rw [UInt64.toNat_sub_of_le _ _ hxy] + rw [Nat.cast_sub (by + rw [UInt64.le_iff_toNat_le] at hxy + exact hxy)] + · rw [if_neg hxy] + have hx_lt_y : x.toNat < y.toNat := by + have hnot : ¬y.toNat ≤ x.toNat := by + intro hle + apply hxy + rw [UInt64.le_iff_toNat_le] + exact hle + exact Nat.lt_of_not_ge hnot + have hraw_lt_size : 2 ^ 64 - y.toNat + x.toNat < 2 ^ 64 := by + have hx_lt_size : x.toNat < 2 ^ 64 := by + simpa [UInt64.size] using UInt64.toNat_lt_size x + have hy_lt_size : y.toNat < 2 ^ 64 := by + simpa [UInt64.size] using UInt64.toNat_lt_size y + omega + have hraw_toNat : + (x - y).toNat = 2 ^ 64 - y.toNat + x.toNat := by + rw [UInt64.toNat_sub] + exact Nat.mod_eq_of_lt hraw_lt_size + have hneg_le_raw_nat : neg_modulus.toNat ≤ (x - y).toNat := by + rw [hraw_toNat, neg_modulus_toNat] + have hy_lt_field : y.toNat < 2 ^ 64 - 2 ^ 32 + 1 := by + simpa [Goldilocks.Basic.fieldSize] using hy + omega + have hneg_le_raw : neg_modulus ≤ x - y := by + rw [UInt64.le_iff_toNat_le] + exact hneg_le_raw_nat + rw [UInt64.toNat_sub_of_le _ _ hneg_le_raw, hraw_toNat, neg_modulus_toNat] + change (((UInt64.size - y.toNat + x.toNat - neg_modulus.toNat : Nat) : + Goldilocks.Basic.Field) = + (x.toNat : Goldilocks.Basic.Field) - (y.toNat : Goldilocks.Basic.Field)) + have hneg_le_concrete : neg_modulus.toNat ≤ UInt64.size - y.toNat + x.toNat := by + simpa [UInt64.size, hraw_toNat] using hneg_le_raw_nat + rw [Nat.cast_sub hneg_le_concrete] + rw [Nat.cast_add] + rw [Nat.cast_sub (Nat.le_of_lt (UInt64.toNat_lt_size y))] + rw [uint64_cast_eq_neg_modulus] + ring + +end Reduction +end Fast +end Goldilocks diff --git a/CompPoly/Fields/Goldilocks/Fast/Theorems.lean b/CompPoly/Fields/Goldilocks/Fast/Theorems.lean new file mode 100644 index 00000000..f9c3a212 --- /dev/null +++ b/CompPoly/Fields/Goldilocks/Fast/Theorems.lean @@ -0,0 +1,321 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Varun Thakore +-/ + +import CompPoly.Fields.Goldilocks.Fast.Arithmetic +import Mathlib.FieldTheory.Finite.Basic + +/-! +# Correctness Theorems for Fast Goldilocks Arithmetic + +This module proves that the fast `UInt64` operations agree with the canonical +`ZMod` Goldilocks field. +-/ + +namespace Goldilocks +namespace Fast + +/-- Converting a canonical natural representative to fast form preserves its value. -/ +@[simp] +private theorem toField_ofCanonicalNat (n : Nat) (h : n < Goldilocks.Basic.fieldSize) : + toField (ofCanonicalNat n h) = (n : Goldilocks.Basic.Field) := by + unfold toField toNat ofCanonicalNat + have hn : n < UInt64.size := Nat.lt_trans h fieldSize_lt_uint64Size + rw [UInt64.toNat_ofNat'] + rw [Nat.mod_eq_of_lt (by simpa [UInt64.size] using hn)] + +/-- Converting a canonical natural representative to fast form and reading it back is +the identity. -/ +@[simp] +private theorem toNat_ofCanonicalNat (n : Nat) (h : n < Goldilocks.Basic.fieldSize) : + toNat (ofCanonicalNat n h) = n := by + unfold toNat ofCanonicalNat + have hn : n < UInt64.size := Nat.lt_trans h fieldSize_lt_uint64Size + rw [UInt64.toNat_ofNat'] + exact Nat.mod_eq_of_lt (by simpa [UInt64.size] using hn) + +/-- Converting a natural number to fast form agrees with the same natural cast in the +canonical field. -/ +@[simp] +theorem toField_ofNat (n : Nat) : + toField (ofNat n) = (n : Goldilocks.Basic.Field) := by + unfold ofNat + rw [toField_ofCanonicalNat] + rw [← ZMod.natCast_zmod_val (n : Goldilocks.Basic.Field)] + rw [ZMod.val_natCast] + +/-- Converting a `UInt64` to fast form agrees with casting its natural value into the +canonical field. -/ +@[simp] +theorem toField_ofUInt64 (x : UInt64) : + toField (ofUInt64 x) = (x.toNat : Goldilocks.Basic.Field) := by + unfold toField toNat ofUInt64 + exact reduceUInt64_cast x + +/-- Converting an integer to fast form agrees with casting it into the canonical field. -/ +@[simp] +theorem toField_ofInt (z : Int) : + toField (ofInt z) = (z : Goldilocks.Basic.Field) := by + unfold ofInt ofField + rw [toField_ofCanonicalNat] + exact ZMod.natCast_zmod_val (z : Goldilocks.Basic.Field) + +/-- Converting from the canonical field to fast form and back is the identity. -/ +@[simp] +theorem toField_ofField (x : Goldilocks.Basic.Field) : toField (ofField x) = x := by + unfold ofField + rw [toField_ofCanonicalNat] + exact ZMod.natCast_zmod_val x + +/-- Converting from fast form to the canonical field and back is the identity. -/ +@[simp] +theorem ofField_toField (x : Field) : ofField (toField x) = x := by + apply Subtype.ext + apply UInt64.toNat_inj.mp + change toNat (ofField (toField x)) = toNat x + unfold ofField toField + rw [toNat_ofCanonicalNat] + exact ZMod.val_natCast_of_lt x.property + +/-- The canonical-field interpretation distinguishes fast Goldilocks values. -/ +theorem toField_injective : Function.Injective toField := + Function.LeftInverse.injective ofField_toField + +/-- Fermat-style inversion in the canonical Goldilocks field. -/ +private lemma canonical_inv_eq_pow (a : Goldilocks.Basic.Field) (ha : a ≠ 0) : + a⁻¹ = a ^ (Goldilocks.Basic.fieldSize - 2) := by + have hcard : Fintype.card Goldilocks.Basic.Field = Goldilocks.Basic.fieldSize := + ZMod.card Goldilocks.Basic.fieldSize + have h1 : a ^ (Goldilocks.Basic.fieldSize - 1) = 1 := by + have h := FiniteField.pow_card_sub_one_eq_one a ha + rw [hcard] at h + exact h + have hmul : a * a ^ (Goldilocks.Basic.fieldSize - 2) = 1 := by + rw [← pow_succ'] + show a ^ (Goldilocks.Basic.fieldSize - 2 + 1) = 1 + have : Goldilocks.Basic.fieldSize - 2 + 1 = Goldilocks.Basic.fieldSize - 1 := by + unfold Goldilocks.Basic.fieldSize + omega + rw [this] + exact h1 + exact (eq_inv_of_mul_eq_one_left (by rwa [mul_comm])).symm + +/-- Fast zero maps to canonical zero. -/ +@[simp] +theorem toField_zero : toField (0 : Field) = 0 := by + decide + +/-- Fast one maps to canonical one. -/ +@[simp] +theorem toField_one : toField (1 : Field) = 1 := by + decide + +/-- Fast addition agrees with canonical-field addition. -/ +@[simp] +theorem toField_add (x y : Field) : toField (x + y) = toField x + toField y := by + change + (((add x y).val.toNat : Goldilocks.Basic.Field) = + (x.val.toNat : Goldilocks.Basic.Field) + (y.val.toNat : Goldilocks.Basic.Field)) + unfold add + rw [Reduction.reduceAddWithCarryRaw_cast _ _ + (Reduction.addWithCarry_bound x.val y.val x.property y.property)] + let lo := x.val + y.val + let carry := decide (lo < x.val) + have hvalue := Reduction.addWithCarry_value x.val y.val + change lo.toNat + (if carry then UInt64.size else 0) = x.val.toNat + y.val.toNat at hvalue + change + ((lo.toNat : Goldilocks.Basic.Field) + + (if carry then (UInt64.size : Goldilocks.Basic.Field) else 0) = + (x.val.toNat : Goldilocks.Basic.Field) + (y.val.toNat : Goldilocks.Basic.Field)) + by_cases hcarry : carry = true + · simp only [hcarry, if_true] at hvalue ⊢ + rw [← Nat.cast_add, hvalue, Nat.cast_add] + · simp only [hcarry, Bool.false_eq_true, if_false, add_zero] at hvalue ⊢ + rw [hvalue, Nat.cast_add] + +/-- Fast negation agrees with canonical-field negation. -/ +@[simp] +theorem toField_neg (x : Field) : toField (-x) = -toField x := by + change toField (neg x) = -(toField x) + unfold neg toField toNat + exact Reduction.negRaw_cast x.val x.property + +/-- Fast subtraction agrees with canonical-field subtraction. -/ +@[simp] +theorem toField_sub (x y : Field) : toField (x - y) = toField x - toField y := by + change toField (sub x y) = toField x - toField y + unfold sub toField toNat + exact Reduction.subRaw_cast x.val y.val x.property y.property + +/-- Fast multiplication agrees with canonical-field multiplication. -/ +@[simp] +theorem toField_mul (x y : Field) : toField (x * y) = toField x * toField y := by + change toField (mul x y) = toField x * toField y + unfold mul toField toNat + exact Reduction.reduceMulRaw_cast x.val y.val + +/-- The named fast multiplication function agrees with canonical-field multiplication. -/ +@[simp] +theorem toField_mul_def (x y : Field) : toField (mul x y) = toField x * toField y := + toField_mul x y + +/-- Fast squaring agrees with multiplying the canonical field value by itself. -/ +@[simp] +theorem toField_square (x : Field) : toField (square x) = toField x * toField x := by + change toField (x * x) = toField x * toField x + rw [toField_mul] + +/-- Repeated fast squaring agrees with raising to `2^n` in the canonical field. -/ +@[simp] +theorem toField_squareN (x : Field) (n : Nat) : + toField (squareN x n) = toField x ^ (2 ^ n) := by + induction n generalizing x with + | zero => + unfold squareN + simp + | succ n ih => + unfold squareN + rw [toField_square, ih] + rw [← pow_add] + congr 1 + rw [Nat.pow_succ] + omega + +/-- Fast multiplication is associative, proved by transporting to the canonical field. -/ +private theorem mul_assoc_field (x y z : Field) : (x * y) * z = x * (y * z) := by + apply toField_injective + rw [toField_mul, toField_mul, toField_mul, toField_mul] + ring + +/-- Binary exponentiation satisfies the expected successor equation. -/ +private theorem pow_succ (x : Field) (n : Nat) : pow x (n + 1) = pow x n * x := by + unfold pow + letI : Semigroup Field := { + mul := (· * ·) + mul_assoc := mul_assoc_field + } + exact npowBinRec_succ n x + +/-- Fast natural-power computation agrees with powers in the canonical field. -/ +@[simp] +theorem toField_pow (x : Field) (n : Nat) : toField (pow x n) = toField x ^ n := by + induction n with + | zero => + unfold pow + rw [npowBinRec_zero] + rw [toField_one] + simp + | succ n ih => + rw [pow_succ, toField_mul, ih, _root_.pow_succ] + +/-- The optimized inversion chain computes the Fermat inverse exponent. -/ +private theorem toField_inv_chain (x : Field) : + toField (inv x) = toField x ^ invExponent := by + unfold inv + simp only [toField_mul_def, toField_square, toField_squareN] + ring_nf + simp [invExponent, Goldilocks.Basic.fieldSize] + +/-- Fast inversion agrees with canonical inversion before notation is unfolded. -/ +private theorem toField_inv_raw (x : Field) : toField (inv x) = (toField x)⁻¹ := by + rw [toField_inv_chain] + by_cases hx : toField x = 0 + · rw [hx] + simp [invExponent, Goldilocks.Basic.fieldSize] + · simpa [invExponent] using (canonical_inv_eq_pow (toField x) hx).symm + +/-- Fast inversion agrees with inversion in the canonical field. -/ +@[simp] +theorem toField_inv (x : Field) : toField x⁻¹ = (toField x)⁻¹ := by + change toField (inv x) = (toField x)⁻¹ + exact toField_inv_raw x + +/-- Division is multiplication by inverse at the level of canonical interpretation. -/ +private theorem toField_div_mul_inv (x y : Field) : + toField (div x y) = toField x * toField (inv y) := by + unfold div + change toField (x * inv y) = toField x * toField (inv y) + exact toField_mul x (inv y) + +/-- Fast division agrees with division in the canonical field. -/ +@[simp] +theorem toField_div (x y : Field) : toField (x / y) = toField x / toField y := by + change toField (div x y) = toField x / toField y + rw [toField_div_mul_inv, toField_inv_raw y] + rfl + +/-- Natural casts in the fast field agree with natural casts in the canonical field. -/ +@[simp] +theorem toField_natCast (n : Nat) : toField (n : Field) = (n : Goldilocks.Basic.Field) := by + change toField (ofNat n) = (n : Goldilocks.Basic.Field) + rw [toField_ofNat] + +/-- Integer casts in the fast field agree with integer casts in the canonical field. -/ +@[simp] +theorem toField_intCast (n : Int) : toField (n : Field) = (n : Goldilocks.Basic.Field) := by + change toField (ofInt n) = (n : Goldilocks.Basic.Field) + rw [toField_ofInt] + +/-- Fast natural scalar multiplication agrees with canonical-field scalar multiplication. -/ +@[simp] +theorem toField_nsmul (n : Nat) (x : Field) : toField (n • x) = n • toField x := by + change toField ((n : Field) * x) = n • toField x + rw [toField_mul, toField_natCast] + rw [nsmul_eq_mul] + +/-- Fast integer scalar multiplication agrees with canonical-field scalar multiplication. -/ +@[simp] +theorem toField_zsmul (n : Int) (x : Field) : toField (n • x) = n • toField x := by + change toField ((n : Field) * x) = n • toField x + rw [toField_mul, toField_intCast] + rw [zsmul_eq_mul] + +/-- Standard fast natural powers agree with powers in the canonical field. -/ +@[simp] +theorem toField_npow (x : Field) (n : Nat) : toField (x ^ n) = toField x ^ n := by + change toField (pow x n) = toField x ^ n + rw [toField_pow] + +/-- Standard fast integer powers agree with integer powers in the canonical field. -/ +@[simp] +theorem toField_zpow (x : Field) (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)⁻¹ := toField_inv_raw x + rw [toField_pow, hinv, zpow_negSucc, inv_pow] + +/-- Nonnegative rational casts in the fast field agree with canonical-field casts. -/ +@[simp] +theorem toField_nnratCast (q : ℚ≥0) : toField (q : Field) = (q : Goldilocks.Basic.Field) := by + change toField (ofField (q : Goldilocks.Basic.Field)) = (q : Goldilocks.Basic.Field) + rw [toField_ofField] + +/-- Rational casts in the fast field agree with canonical-field casts. -/ +@[simp] +theorem toField_ratCast (q : ℚ) : toField (q : Field) = (q : Goldilocks.Basic.Field) := by + change toField (ofField (q : Goldilocks.Basic.Field)) = (q : Goldilocks.Basic.Field) + rw [toField_ofField] + +/-- Fast nonnegative rational scalar multiplication agrees with canonical-field scalar +multiplication. -/ +@[simp] +theorem toField_nnqsmul (q : ℚ≥0) (x : Field) : toField (q • x) = q • toField x := by + change toField (ofField (q • toField x)) = q • toField x + rw [toField_ofField] + +/-- Fast rational scalar multiplication agrees with canonical-field scalar multiplication. -/ +@[simp] +theorem toField_qsmul (q : ℚ) (x : Field) : toField (q • x) = q • toField x := by + change toField (ofField (q • toField x)) = q • toField x + rw [toField_ofField] + + +end Fast +end Goldilocks From 8f5df6271592c5b90952f0c5ee258b2e942ed278 Mon Sep 17 00:00:00 2001 From: Varun Thakore Date: Tue, 30 Jun 2026 17:43:12 +0530 Subject: [PATCH 2/3] feat(fields): add external mul for goldilocks --- CompPoly.lean | 8 ++ CompPoly/Fields/Goldilocks/FastExt.lean | 106 ++++++++++++++++++++++++ lakefile.lean | 26 +++++- native/goldilocks_native.c | 44 ++++++++++ 4 files changed, 183 insertions(+), 1 deletion(-) create mode 100644 CompPoly/Fields/Goldilocks/FastExt.lean create mode 100644 native/goldilocks_native.c diff --git a/CompPoly.lean b/CompPoly.lean index 58d6f008..b8c3f9e6 100644 --- a/CompPoly.lean +++ b/CompPoly.lean @@ -99,6 +99,14 @@ import CompPoly.Fields.Binary.Tower.Support.LinearIndependentFin2 import CompPoly.Fields.Binary.Tower.Support.Preliminaries import CompPoly.Fields.Binary.Tower.TensorAlgebra import CompPoly.Fields.Goldilocks +import CompPoly.Fields.Goldilocks.Basic +import CompPoly.Fields.Goldilocks.Fast +import CompPoly.Fields.Goldilocks.Fast.Arithmetic +import CompPoly.Fields.Goldilocks.Fast.Field +import CompPoly.Fields.Goldilocks.Fast.Internal +import CompPoly.Fields.Goldilocks.Fast.Reduction +import CompPoly.Fields.Goldilocks.Fast.Theorems +import CompPoly.Fields.Goldilocks.FastExt import CompPoly.Fields.KoalaBear import CompPoly.Fields.KoalaBear.Basic import CompPoly.Fields.KoalaBear.Fast diff --git a/CompPoly/Fields/Goldilocks/FastExt.lean b/CompPoly/Fields/Goldilocks/FastExt.lean new file mode 100644 index 00000000..ac60d57f --- /dev/null +++ b/CompPoly/Fields/Goldilocks/FastExt.lean @@ -0,0 +1,106 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Varun Thakore +-/ + +import CompPoly.Fields.Goldilocks.Fast.Arithmetic + +/-! +# Extern-Backed Fast Goldilocks Operations + +This module provides opt-in extern-backed multiplication and inversion routines +for the fast Goldilocks representation. + +The verified implementation in `CompPoly.Fields.Goldilocks.Fast` remains the +default field implementation. The functions here move selected low-level +operations to external code for performance, so their arithmetic correctness +depends on the linked extern implementations. + +These declarations are a trusted native boundary: Lean checks that the wrappers +have the stated types, but it does not verify the C implementations behind the +`@[extern]` symbols. Use this module when native performance is needed and the +linked C code is accepted as part of the trusted computing base. +-/ + +namespace Goldilocks +namespace Fast +namespace Ext + +/-- High 64 bits of a 64-by-64 unsigned multiplication. + +Together with Lean's native wrapping multiplication, this gives the two words of +the full 128-bit product. + +This is an extern performance primitive. Its correctness is trusted from the +linked native implementation in `native/goldilocks_native.c`. +-/ +@[extern "lean_uint64_mul_hi"] +opaque mulHi (a b : UInt64) : UInt64 + +/-- Extern Goldilocks multiplication returning a canonical `UInt64` representative. + +This is an extern performance primitive. Lean does not verify that the native +routine implements Goldilocks multiplication correctly. The native symbol is +implemented in `native/goldilocks_native.c`. +-/ +@[extern "lean_goldilocks_mul"] +opaque goldilocksMul (a b : UInt64) : UInt64 + +/-- Multiplication using extern `UInt64.mulHi` for the high word and Lean reduction. -/ +@[inline] +def mulWithMulHi (x y : Field) : Field := + let lo := x.val * y.val + let hi := mulHi x.val y.val + ⟨Reduction.reduceUInt128Raw lo hi, Reduction.reduceUInt128Raw_lt lo hi⟩ + +/-- Multiplication using an extern Goldilocks multiplication primitive. + +The external function is expected to return a canonical representative; this +wrapper still passes the result through `reduceUInt64` to inhabit the fast field +type safely. +-/ +@[inline] +def mulNative (x y : Field) : Field := + reduceUInt64 (goldilocksMul x.val y.val) + +/-- Squaring using extern `UInt64.mulHi` for multiplication. -/ +@[inline] +def squareWithMulHi (x : Field) : Field := + mulWithMulHi x x + +/-- Squaring using the extern Goldilocks multiplication primitive. -/ +@[inline] +def squareNative (x : Field) : Field := + mulNative x x + +/-- Repeated squaring using extern-backed multiplication. -/ +@[inline] +def squareNNative (x : Field) : Nat → Field + | 0 => x + | n + 1 => squareNNative (squareNative x) n + +/-- Inversion using the Goldilocks `p - 2` addition chain and extern multiplication. -/ +@[noinline] +def invNative (x : Field) : Field := + let t2 := mulNative (squareNative x) x + let t4 := mulNative (squareNNative t2 2) t2 + let t8 := mulNative (squareNNative t4 4) t4 + let t16 := mulNative (squareNNative t8 8) t8 + let t31 := + mulNative (squareNNative t16 15) + (mulNative (squareNNative t8 7) + (mulNative (squareNNative t4 3) + (mulNative (squareNative t2) x))) + let t32m2 := squareNative t31 + let t32m1 := mulNative t32m2 x + mulNative (squareNNative t32m2 32) t32m1 + +/-- Division using extern-backed inversion and multiplication. -/ +@[inline] +def divNative (x y : Field) : Field := + mulNative x (invNative y) + +end Ext +end Fast +end Goldilocks diff --git a/lakefile.lean b/lakefile.lean index 16d958cd..a4d573be 100644 --- a/lakefile.lean +++ b/lakefile.lean @@ -8,8 +8,32 @@ package CompPoly where require "leanprover-community" / mathlib @ git "v4.30.0" +def nativeDir : FilePath := __dir__ / "native" + +def nativeBuildDir : FilePath := __dir__ / ".lake" / "build" / "native" + +def nativeLib (name src : String) (extraCcArgs : Array String := #[]) : FetchM (Job FilePath) := do + let srcFile := nativeDir / src + let oFile := nativeBuildDir / s!"{name}.o" + let libFile := nativeBuildDir / s!"lib{name}.a" + let srcJob ← inputTextFile srcFile + buildFileAfterDep libFile srcJob fun _srcTrace => do + compileO oFile srcFile (#["-O3", "-I", (← getLeanIncludeDir).toString] ++ extraCcArgs) + createParentDirs libFile + removeFileIfExists libFile + proc { + cmd := (← getLeanAr).toString + args := #["rcs", libFile.toString, oFile.toString] + } + +extern_lib liblean_goldilocks_native _pkg := nativeLib "goldilocks_native" "goldilocks_native.c" + +def nativeLinkArgs (name : String) : Array String := + #["-L", nativeBuildDir.toString, s!"-l{name}"] + @[default_target] -lean_lib CompPoly +lean_lib CompPoly where + moreLinkArgs := nativeLinkArgs "goldilocks_native" lean_lib CompPolyTests where srcDir := "tests" diff --git a/native/goldilocks_native.c b/native/goldilocks_native.c new file mode 100644 index 00000000..e899d1c5 --- /dev/null +++ b/native/goldilocks_native.c @@ -0,0 +1,44 @@ +#include +#include + +/* + * Native performance primitives for Goldilocks fast arithmetic. + * + * These functions are called from Lean through @[extern] declarations in + * CompPoly/Fields/Goldilocks/FastExt.lean. Lean verifies only the types of the + * extern declarations, not the arithmetic performed here. This file is + * therefore part of the trusted native boundary for the extern-backed API. + * + * The verified Lean implementation remains the default field implementation; + * this native code is used only by the opt-in FastExt module. + */ + +// Return the high 64 bits of a 64-by-64 unsigned multiplication. +LEAN_EXPORT uint64_t lean_uint64_mul_hi(uint64_t a, uint64_t b) { + return (uint64_t)(((unsigned __int128)a * b) >> 64); +} + +static const uint64_t GOLDILOCKS_NEG_MODULUS = 0x00000000FFFFFFFFULL; + +static inline uint64_t goldilocks_add_no_canonicalize(uint64_t x, uint64_t y) { + uint64_t res = x + y; + return res + (res < x ? GOLDILOCKS_NEG_MODULUS : 0); +} + +// Goldilocks multiplication modulo p = 2^64 - 2^32 + 1. +LEAN_EXPORT uint64_t lean_goldilocks_mul(uint64_t a, uint64_t b) { + unsigned __int128 prod = (unsigned __int128)a * b; + uint64_t lo = (uint64_t)prod; + uint64_t hi = (uint64_t)(prod >> 64); + + uint64_t hi_hi = hi >> 32; + uint64_t hi_lo = hi & GOLDILOCKS_NEG_MODULUS; + + uint64_t t0 = lo - hi_hi; + if (lo < hi_hi) { + t0 -= GOLDILOCKS_NEG_MODULUS; + } + + uint64_t t1 = hi_lo * GOLDILOCKS_NEG_MODULUS; + return goldilocks_add_no_canonicalize(t0, t1); +} From 51c678b2a02e9c73937d80045338d9f1862bb807 Mon Sep 17 00:00:00 2001 From: Varun Thakore Date: Tue, 30 Jun 2026 18:20:08 +0530 Subject: [PATCH 3/3] tests(fields): add tests for goldilocks --- lakefile.lean | 6 ++ tests/CompPolyTests.lean | 1 + .../CompPolyTests/Fields/Goldilocks/Fast.lean | 46 ++++++++++++ .../Fields/Goldilocks/FastExt.lean | 71 +++++++++++++++++++ 4 files changed, 124 insertions(+) create mode 100644 tests/CompPolyTests/Fields/Goldilocks/Fast.lean create mode 100644 tests/CompPolyTests/Fields/Goldilocks/FastExt.lean diff --git a/lakefile.lean b/lakefile.lean index a4d573be..270eaf26 100644 --- a/lakefile.lean +++ b/lakefile.lean @@ -37,6 +37,7 @@ lean_lib CompPoly where lean_lib CompPolyTests where srcDir := "tests" + moreLinkArgs := nativeLinkArgs "goldilocks_native" lean_lib CompPolyBenchLib where srcDir := "bench" @@ -44,3 +45,8 @@ lean_lib CompPolyBenchLib where lean_exe CompPolyBench where srcDir := "bench" + +lean_exe CompPolyGoldilocksFastExtTests where + srcDir := "tests" + root := `CompPolyTests.Fields.Goldilocks.FastExt + moreLinkArgs := nativeLinkArgs "goldilocks_native" diff --git a/tests/CompPolyTests.lean b/tests/CompPolyTests.lean index 0321cad5..0b6c0036 100644 --- a/tests/CompPolyTests.lean +++ b/tests/CompPolyTests.lean @@ -23,6 +23,7 @@ import CompPolyTests.Bivariate.WeightedDegree import CompPolyTests.Data.MvPolynomial.Notation import CompPolyTests.Fields.Binary.AdditiveNTT.NovelPolynomialBasis import CompPolyTests.Fields.Binary.BF128Ghash.Prelude +import CompPolyTests.Fields.Goldilocks.Fast import CompPolyTests.Fields.KoalaBear.Fast import CompPolyTests.Fields.PrattCertificate import CompPolyTests.LinearAlgebra.Dense diff --git a/tests/CompPolyTests/Fields/Goldilocks/Fast.lean b/tests/CompPolyTests/Fields/Goldilocks/Fast.lean new file mode 100644 index 00000000..79e976ca --- /dev/null +++ b/tests/CompPolyTests/Fields/Goldilocks/Fast.lean @@ -0,0 +1,46 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Varun Thakore +-/ + +import CompPoly.Fields.Goldilocks.Fast + +/-! +# Fast Goldilocks Field Tests + +Regression checks for the verified `UInt64` fast representation. + +Run this file with: +`lake build CompPolyTests.Fields.Goldilocks.Fast` +-/ + +namespace Goldilocks.Fast + +private def p : Nat := Goldilocks.Basic.fieldSize + +#guard raw (0 : Field) = 0 +#guard raw (1 : Field) = 1 +#guard toNat (ofNat 73) = 73 +#guard toNat (ofNat p) = 0 +#guard toNat (ofNat (p + 73)) = 73 +#guard toNat (ofUInt64 (UInt64.ofNat (UInt64.size - 1))) = 4294967294 +#guard toNat ((ofNat (p - 1)) + (4 : Field)) = 3 +#guard toNat ((ofNat (p - 1)) + (ofNat (p - 1))) = p - 2 +#guard toNat ((17 : Field) - (6 : Field)) = 11 +#guard toNat ((6 : Field) - (17 : Field)) = p - 11 +#guard toNat (-(0 : Field)) = 0 +#guard toNat (-(1 : Field)) = p - 1 +#guard toNat ((ofNat (p - 1)) * (ofNat (p - 1))) = 1 +#guard toField (square (54321 : Field)) = ((54321 : Goldilocks.Basic.Field) ^ 2) +#guard toNat ((73 : Field) ^ 0) = 1 +#guard toNat ((73 : Field) ^ 1) = 73 +#guard toField ((987654321 : Field) ^ 19) = ((987654321 : Goldilocks.Basic.Field) ^ 19) +#guard toField ((987654321 : Field) ^ 511) = ((987654321 : Goldilocks.Basic.Field) ^ 511) +#guard toNat ((0 : Field)⁻¹) = 0 +#guard toNat ((73 : Field)⁻¹ * (73 : Field)) = 1 +#guard toNat ((73 : Field) / (73 : Field)) = 1 +#guard toField ((73 : Field)⁻¹) = ((73 : Goldilocks.Basic.Field)⁻¹) +#guard toField ((73 : Field) ^ (-5 : Int)) = ((73 : Goldilocks.Basic.Field) ^ (-5 : Int)) + +end Goldilocks.Fast diff --git a/tests/CompPolyTests/Fields/Goldilocks/FastExt.lean b/tests/CompPolyTests/Fields/Goldilocks/FastExt.lean new file mode 100644 index 00000000..6ec24580 --- /dev/null +++ b/tests/CompPolyTests/Fields/Goldilocks/FastExt.lean @@ -0,0 +1,71 @@ +/- +Copyright (c) 2026 CompPoly Contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Varun Thakore +-/ + +import CompPoly.Fields.Goldilocks.Fast +import CompPoly.Fields.Goldilocks.FastExt + +/-! +# Extern-Backed Fast Goldilocks Tests + +Runtime regression checks for `Goldilocks.Fast.Ext`. These tests compare the +extern-backed operations against the verified `Goldilocks.Fast` implementation. + +The checks live in an executable rather than `#guard`s because Lean's module +interpreter cannot call project-local C externs during elaboration. + +Run this file with: +`lake exe CompPolyGoldilocksFastExtTests` +-/ + +namespace CompPolyTests.Fields.Goldilocks.FastExt + +private def check (name : String) (ok : Bool) : IO Bool := do + if ok then + return true + else + IO.eprintln s!"failed: {name}" + return false + +private def a : Goldilocks.Fast.Field := 987654321 + +private def b : Goldilocks.Fast.Field := 54321 + +private def c : Goldilocks.Fast.Field := 73 + +private def nearTop : Goldilocks.Fast.Field := + Goldilocks.Fast.ofNat (Goldilocks.Basic.fieldSize - 1) + +private def maxUInt64 : Goldilocks.Fast.Field := + Goldilocks.Fast.ofUInt64 (UInt64.ofNat (UInt64.size - 1)) + +private def runChecks : IO Bool := do + let ok1 ← check "mulWithMulHi" (Goldilocks.Fast.Ext.mulWithMulHi a b = a * b) + let ok2 ← check "mulNative" (Goldilocks.Fast.Ext.mulNative a b = a * b) + let ok3 ← + check "squareWithMulHi" (Goldilocks.Fast.Ext.squareWithMulHi a = Goldilocks.Fast.square a) + let ok4 ← + check "squareNative" (Goldilocks.Fast.Ext.squareNative a = Goldilocks.Fast.square a) + let ok5 ← + check "squareNNative" (Goldilocks.Fast.Ext.squareNNative a 8 = Goldilocks.Fast.squareN a 8) + let ok6 ← check "invNative" (Goldilocks.Fast.Ext.invNative c = c⁻¹) + let ok7 ← check "divNative" (Goldilocks.Fast.Ext.divNative a c = a / c) + let ok8 ← + check "mulNative near modulus" + (Goldilocks.Fast.Ext.mulNative nearTop nearTop = nearTop * nearTop) + let ok9 ← + check "mulNative max UInt64" (Goldilocks.Fast.Ext.mulNative maxUInt64 a = maxUInt64 * a) + let ok10 ← + check "divNative nontrivial" (Goldilocks.Fast.Ext.divNative b c = b / c) + return ok1 && ok2 && ok3 && ok4 && ok5 && ok6 && ok7 && ok8 && ok9 && ok10 + +end CompPolyTests.Fields.Goldilocks.FastExt + +/-- Run the extern-backed Goldilocks regression checks. -/ +def main : IO UInt32 := do + if ← CompPolyTests.Fields.Goldilocks.FastExt.runChecks then + return 0 + else + return 1