Fast inversion for eight-limb Montgomery fields

Correctness of the checked inversion of Montgomery/Native64x8InvDefs: invGcdRaw computes the field inverse and FastField.invGcd is its proof-carrying wrapper. Also proves the divstep coefficient bound and the mac-width safety of the candidate.

Divstep coefficient bounds

theorem Montgomery.Native64x8.gcdInner_zero (a b : UInt64) (f0 g0 f1 g1 : ℤ) : gcdInner 0 a b f0 g0 f1 g1 = (a, b, f0, g0, f1, g1)

Unfolding of the empty divstep run.

theorem Montgomery.Native64x8.gcdInner_succ (rounds : ℕ) (a b : UInt64) (f0 g0 f1 g1 : ℤ) : gcdInner (rounds + 1) a b f0 g0 f1 g1 = if (a &&& 1 == 0) = true then gcdInner rounds (a >>> 1) b f0 g0 (f1 * 2) (g1 * 2) else if a < b then gcdInner rounds ((b - a) >>> 1) a (f1 - f0) (g1 - g0) (f0 * 2) (g0 * 2) else gcdInner rounds ((a - b) >>> 1) b (f0 - f1) (g0 - g1) (f1 * 2) (g1 * 2)

Unfolding of one divstep.

theorem Montgomery.Native64x8.gcdInner_natAbs_le {rounds : ℕ} {a b a' b' : UInt64} {f0 g0 f1 g1 f0' g0' f1' g1' : ℤ} {n : ℕ} (h0 : f0.natAbs + g0.natAbs ≤ n) (h1 : f1.natAbs + g1.natAbs ≤ n) (heq : gcdInner rounds a b f0 g0 f1 g1 = (a', b', f0', g0', f1', g1')) : f0'.natAbs + g0'.natAbs ≤ 2 ^ rounds * n ∧ f1'.natAbs + g1'.natAbs ≤ 2 ^ rounds * n

Transition entries at most double per divstep.

theorem Montgomery.Native64x8.gcdInner_natAbs_le_31 {rounds : ℕ} (h31 : rounds ≤ 31) {a b a' b' : UInt64} {f0' g0' f1' g1' : ℤ} (heq : gcdInner rounds a b 1 0 0 1 = (a', b', f0', g0', f1', g1')) : f0'.natAbs + g0'.natAbs ≤ 2 ^ 31 ∧ f1'.natAbs + g1'.natAbs ≤ 2 ^ 31

Chunks of at most 31 divsteps satisfy the linear-combination bound.

Order from the borrow chain

theorem Montgomery.Native64x8.subBorrow_eq_one_iff {a b : Limbs8} (ha : a.Bounded) (hb : b.Bounded) : subBorrow a b = 1 ↔ a.toNat < b.toNat

The final borrow decides the numeric order.

Mac-width safety of the candidate

theorem Montgomery.Native64x8.gcdLinearCombDiv_bounded (a b : Limbs8) (f g : ℤ) : (gcdLinearCombDiv a b f g).1.Bounded

The division lincomb masks every output limb.

theorem Montgomery.Native64x8.gcdLinearCombMontyRed_lt {q : Limbs8} {negInv : UInt64} {a b : Limbs8} {f g : ℤ} (hq : q.Bounded) (hq0 : 0 < q.toNat) (hq2 : 2 * q.toNat < 2 ^ 256) (hn : negInv.toNat < 2 ^ 32) (hnq : negInv.toNat * q.toNat % 2 ^ 32 = 2 ^ 32 - 1) (ha : a.Bounded) (hb : b.Bounded) (haq : a.toNat < q.toNat) (hbq : b.toNat < q.toNat) (hfg : f.natAbs + g.natAbs ≤ 2 ^ 31) : (gcdLinearCombMontyRed q negInv a b f g).Bounded ∧ (gcdLinearCombMontyRed q negInv a b f g).toNat < q.toNat

The Montgomery lincomb stays canonical for canonical inputs.

theorem Montgomery.Native64x8.gcdMainLoop_bounded {q : Limbs8} {negInv : UInt64} {rounds : ℕ} {a u b v A U B V : Limbs8} (hq : q.Bounded) (hq0 : 0 < q.toNat) (hq2 : 2 * q.toNat < 2 ^ 256) (hn : negInv.toNat < 2 ^ 32) (hnq : negInv.toNat * q.toNat % 2 ^ 32 = 2 ^ 32 - 1) (ha : a.Bounded) (hb : b.Bounded) (hu : u.Bounded) (huq : u.toNat < q.toNat) (hv : v.Bounded) (hvq : v.toNat < q.toNat) (heq : gcdMainLoop q negInv rounds a u b v = (A, U, B, V)) : A.Bounded ∧ B.Bounded ∧ U.Bounded ∧ U.toNat < q.toNat ∧ V.Bounded ∧ V.toNat < q.toNat

The main loop keeps both tracks at mac width and the Montgomery pair canonical.

theorem Montgomery.Native64x8.gcdFinalChunks_lt {q : Limbs8} {negInv : UInt64} {finalRounds : ℕ} {a u b v : Limbs8} (hfr : finalRounds ≤ 62) (hq : q.Bounded) (hq0 : 0 < q.toNat) (hq2 : 2 * q.toNat < 2 ^ 256) (hn : negInv.toNat < 2 ^ 32) (hnq : negInv.toNat * q.toNat % 2 ^ 32 = 2 ^ 32 - 1) (hu : u.Bounded) (huq : u.toNat < q.toNat) (hv : v.Bounded) (hvq : v.toNat < q.toNat) : (gcdFinalChunks q negInv finalRounds a u b v).Bounded ∧ (gcdFinalChunks q negInv finalRounds a u b v).toNat < q.toNat

The final chunks stay canonical for canonical inputs.

theorem Montgomery.Native64x8.gcdInvCandidate_lt {modulus : ℕ} [P : GcdData modulus] {q x : Limbs8} {negInv : UInt64} (hq : q.Bounded) (hqm : q.toNat = modulus) (hq0 : 0 < q.toNat) (hq2 : 2 * q.toNat < 2 ^ 256) (hn : negInv.toNat < 2 ^ 32) (hnq : negInv.toNat * q.toNat % 2 ^ 32 = 2 ^ 32 - 1) (hx : x.Bounded) : (gcdInvCandidate modulus q negInv x).Bounded ∧ (gcdInvCandidate modulus q negInv x).toNat < q.toNat

The candidate stays at mac width and canonical.

theorem Montgomery.Native64x8.gcdInvCandidate_lt_modulus {modulus : ℕ} [P : Mont64x8Field modulus] [GcdData modulus] {x : Limbs8} (hx : x.Bounded) : (gcdInvCandidate modulus (Mont64x8Field.modulusLimbs modulus) (Mont64x8Field.montgomeryNegInv modulus) x).Bounded ∧ (gcdInvCandidate modulus (Mont64x8Field.modulusLimbs modulus) (Mont64x8Field.montgomeryNegInv modulus) x).toNat < modulus

With the class data, the candidate is bounded and canonical for any bounded input.

The Fermat fallback and the checked raw inversion

theorem Montgomery.Native64x8.montPow_eq_mul_pow {modulus : ℕ} [P : Mont64x8Field modulus] (acc x : FastField modulus) (n : ℕ) : montPow (Mont64x8Field.modulusLimbs modulus) (Mont64x8Field.montgomeryNegInv modulus) (↑acc) (↑x) n = ↑(acc * x ^ n)

montPow computes acc · xⁿ.

theorem Montgomery.Native64x8.montPow_eq_inv {modulus : ℕ} [P : Mont64x8Field modulus] (x : FastField modulus) : montPow (Mont64x8Field.modulusLimbs modulus) (Mont64x8Field.montgomeryNegInv modulus) (Mont64x8Field.rModModulus modulus) (↑x) (modulus - 2) = ↑x⁻¹

The montPow fallback computes the field inverse.

theorem Montgomery.Native64x8.invGcdRaw_eq_inv {modulus : ℕ} [P : Mont64x8Field modulus] [GcdData modulus] (x : FastField modulus) : invGcdRaw modulus (Mont64x8Field.modulusLimbs modulus) (Mont64x8Field.montgomeryNegInv modulus) (Mont64x8Field.rModModulus modulus) ↑x = ↑x⁻¹

invGcdRaw computes the field inverse.

The proof-carrying wrapper

@[inline]

def Montgomery.Native64x8.FastField.invGcd {modulus : ℕ} [P : Mont64x8Field modulus] [GcdData modulus] (x : FastField modulus) : FastField modulus

The proof-carrying wrapper of the checked inversion invGcdRaw.

@[simp]

theorem Montgomery.Native64x8.FastField.invGcd_eq_inv {modulus : ℕ} [P : Mont64x8Field modulus] [GcdData modulus] (x : FastField modulus) : x.invGcd = x⁻¹

invGcd agrees with the Field inverse.