From 7a7aa2604c77496a9b3927bb11c67ff290dfe8c4 Mon Sep 17 00:00:00 2001 From: japabu Date: Thu, 22 Oct 2026 12:52:28 +0200 Subject: [PATCH 1/2] crypto: a P-256 base-point comb, odd representatives, and a faster field P-256 key generation ran the variable-base Booth windows from G, with 255 doublings; every addition in the window loop also computed a masked doubling; squaring was a multiplication; inversion was square-and-multiply per bit; ECDSA verify built G's odd multiples on every call; and the scalar multiplications ran on identity values for scalars whose top windows are zero. - Comb: P-256's base point multiplication reads row i of a compile-time table j 2^(5i) G (52 rows of 16 affine points, 53 KiB, built by const Montgomery arithmetic with one simultaneous inversion), one masked 16-entry scan and one mixed addition (8M + 3S) per window, with no doublings. - Odd representatives: both secret multiplications run on k + j n, the j of m and m + 1 that makes it odd, m = 1 for P-256 and 2 for P-384. Its top Booth digit is nonzero, so neither the comb nor the doubling chain ever holds the identity, and a zero digit adds the row's first entry and discards the sum by mask instead of adding an identity of zeros. - Exceptional additions: bounds on the running multiple, stated at both sites, show no addition but the last can meet equal or opposite operands. The last one doubles by mask; every other addition is the incomplete formula, which removes 51 masked doublings per ECDH. - Field: P-256's multiplication, squaring, addition and subtraction are written out word by word (the shipped profile left the generic loops rolled); squaring computes the cross products once and reduces the 8-word square with the special step; inversion is an addition chain of 255 squarings and 12 multiplications for P-256's prime and 4-bit windows for the other moduli. Const forms of the arithmetic build the tables. - Verify: u1 is recoded at width 7 over a static table of G, 3G, ..., 63G, added in affine form; Q's table builds with public additions. Exploratory, on the development Mac in a shipped-like profile: P-256 ECDH from 162 to 56 us (ring 33), P-256 verify from 82 to 48 us (ring 30). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01PXEk8BBMDT8tXu4AnENB9Q --- lib/crypto/src/ecc.rs | 4 +- lib/crypto/src/montgomery.rs | 442 ++++++++++++++++++++++- lib/crypto/src/nist.rs | 567 +++++++++++++++++++++++++----- lib/crypto/src/test_arithmetic.rs | 44 +++ 4 files changed, 951 insertions(+), 106 deletions(-) diff --git a/lib/crypto/src/ecc.rs b/lib/crypto/src/ecc.rs index 3c52ff878ff..db225e6ffdb 100644 --- a/lib/crypto/src/ecc.rs +++ b/lib/crypto/src/ecc.rs @@ -9,7 +9,7 @@ macro_rules! curve_api { use crate::{ Error, Random, hash::$hash, - nist::{$curve, Point, Scalar}, + nist::{$curve, Curve, Point, Scalar}, secret::Secret, }; pub struct PublicKey { @@ -56,7 +56,7 @@ macro_rules! curve_api { pub fn public_key(&self) -> PublicKey { // The base point's multiple by a nonzero scalar below n // is finite and on the curve. - let point = Point::base().mul(&self.0); + let point = <$curve as Curve<$limbs>>::base_mul(&self.0); let mut encoded = [0; $encoded]; point.encode(&mut encoded).expect("nonzero prime-order scalar"); PublicKey { point, encoded } diff --git a/lib/crypto/src/montgomery.rs b/lib/crypto/src/montgomery.rs index d8450a434ec..ee94471d20a 100644 --- a/lib/crypto/src/montgomery.rs +++ b/lib/crypto/src/montgomery.rs @@ -87,6 +87,9 @@ impl, const N: usize> Residue { self.mul(&one).0 } pub(crate) fn add(&self, b: &Self) -> Self { + if let Form::P256 = M::FORM { + return Self(p256::add(&self.0, &b.0), PhantomData); + } let mut v = [0; N]; let mut carry = false; for i in 0..N { @@ -98,6 +101,9 @@ impl, const N: usize> Residue { self.add(self) } pub(crate) fn sub(&self, b: &Self) -> Self { + if let Form::P256 = M::FORM { + return Self(p256::sub(&self.0, &b.0), PhantomData); + } let mut v = [0; N]; let mut borrow = false; for i in 0..N { @@ -114,11 +120,55 @@ impl, const N: usize> Residue { /// reduction step (CIOS; Koc, Acar and Kaliski, "Analyzing and comparing /// Montgomery multiplication algorithms", IEEE Micro 1996). pub(crate) fn mul(&self, b: &Self) -> Self { + if let Form::P256 = M::FORM { + return Self(p256::mul(&self.0, &b.0), PhantomData); + } let (t, high) = product::(&self.0, &b.0); Self(reduce(&t, high, &M::P), PhantomData) } + /// self^2 / R mod p: the N(N - 1)/2 cross products once, doubled, the N + /// squares added, then N reduction steps over the 2N-word square. pub(crate) fn square(&self) -> Self { - self.mul(self) + if let Form::P256 = M::FORM { + return Self(p256::square(&self.0), PhantomData); + } + let (t, high) = redc::(&mut square_wide(&self.0)); + Self(reduce(&t, high, &M::P), PhantomData) + } + /// The `const` forms of [`Self::mul`], [`Self::add`] and [`Self::sub`], + /// for tables built at compile time; their final subtraction is a plain + /// select. + pub(crate) const fn times(&self, b: &Self) -> Self { + let (t, high) = product::(&self.0, &b.0); + Self(reduce_const(&t, high, &M::P), PhantomData) + } + pub(crate) const fn plus(&self, b: &Self) -> Self { + let mut v = [0; N]; + let mut carry = false; + let mut i = 0; + while i < N { + (v[i], carry) = carrying_add(self.0[i], b.0[i], carry); + i += 1; + } + Self(reduce_const(&v, carry as u64, &M::P), PhantomData) + } + pub(crate) const fn minus(&self, b: &Self) -> Self { + let mut v = [0; N]; + let mut borrow = false; + let mut i = 0; + while i < N { + (v[i], borrow) = borrowing_sub(self.0[i], b.0[i], borrow); + i += 1; + } + if borrow { + let mut carry = false; + let mut i = 0; + while i < N { + (v[i], carry) = carrying_add(v[i], M::P[i], carry); + i += 1; + } + } + Self(v, PhantomData) } /// One when self is zero, otherwise zero. pub(crate) fn is_zero(&self) -> u64 { @@ -143,23 +193,261 @@ impl, const N: usize> Residue { } Self(v, PhantomData) } - /// self^(p-2), the inverse of a nonzero residue, by square-and-multiply - /// over the public exponent. + /// self^(2^k), with the masked final subtraction. + fn square_times(&self, k: u32) -> Self { + let mut r = *self; + for _ in 0..k { + r = r.square(); + } + r + } + /// The `const` form of [`Self::square_times`]. + const fn squared_times(&self, k: u32) -> Self { + let mut r = *self; + let mut i = 0; + while i < k { + let (t, high) = redc::(&mut square_wide(&r.0)); + r = Self(reduce_const(&t, high, &M::P), PhantomData); + i += 1; + } + r + } + /// self^(p - 2), the inverse of a nonzero residue. The exponent is + /// public, so the schedule is fixed: P-256's prime by an addition chain of + /// 255 squarings and 12 multiplications, other moduli by 4-bit windows. pub(crate) fn invert(&self) -> Self { - let mut e = M::P; - let mut borrow = false; - (e[0], borrow) = borrowing_sub(e[0], 2, borrow); - for word in &mut e[1..] { - (*word, borrow) = borrowing_sub(*word, 0, borrow); + match M::FORM { + Form::P256 => inversion_chain!(self, square_times, mul), + Form::Generic => { + let e = exponent::(); + let mut powers = [Self::one(); 16]; + for j in 1..16 { + powers[j] = powers[j - 1].mul(self); + } + let mut r = Self::one(); + for nibble in (0..N * 16).rev() { + r = r.square_times(4); + let digit = (e[nibble / 16] >> (4 * (nibble % 16))) & 15; + if digit != 0 { + r = r.mul(&powers[digit as usize]); + } + } + r + } } - let mut r = Self::one(); - for i in (0..N * 64).rev() { - r = r.square(); - if (e[i / 64] >> (i % 64)) & 1 == 1 { - r = r.mul(self); + } + /// The `const` form of [`Self::invert`]. + pub(crate) const fn inverted(&self) -> Self { + match M::FORM { + Form::P256 => inversion_chain!(self, squared_times, times), + Form::Generic => { + let e = exponent::(); + let mut powers = [Self::one(); 16]; + let mut j = 1; + while j < 16 { + powers[j] = powers[j - 1].times(self); + j += 1; + } + let mut r = Self::one(); + let mut nibble = N * 16; + while nibble > 0 { + nibble -= 1; + r = r.squared_times(4); + let digit = (e[nibble / 16] >> (4 * (nibble % 16))) & 15; + if digit != 0 { + r = r.times(&powers[digit as usize]); + } + } + r } } - r + } +} + +/// x^(p - 2) for P-256's prime p: p - 2 = 2^256 - 2^224 + 2^192 + 2^96 - 3 +/// reads, from the top, 32 ones, 31 zeros and a one, 96 zeros, 94 ones, a +/// zero and a one. +macro_rules! inversion_chain { + ($x:expr, $square_times:ident, $mul:ident) => {{ + let x = $x; + let x2 = x.$square_times(1).$mul(x); + let x3 = x2.$square_times(1).$mul(x); + let x6 = x3.$square_times(3).$mul(&x3); + let x12 = x6.$square_times(6).$mul(&x6); + let x15 = x12.$square_times(3).$mul(&x3); + let x30 = x15.$square_times(15).$mul(&x15); + let x32 = x30.$square_times(2).$mul(&x2); + let r = x32.$square_times(32).$mul(x); + let r = r.$square_times(128).$mul(&x32); + let r = r.$square_times(32).$mul(&x32); + let r = r.$square_times(30).$mul(&x30); + r.$square_times(2).$mul(x) + }}; +} +use inversion_chain; + +/// p - 2 for M's p. +const fn exponent, const N: usize>() -> [u64; N] { + let mut e = M::P; + let mut borrow; + (e[0], borrow) = borrowing_sub(e[0], 2, false); + let mut i = 1; + while i < N { + (e[i], borrow) = borrowing_sub(e[i], 0, borrow); + i += 1; + } + e +} + +/// P-256's field arithmetic written out word by word: the shipped profile +/// does not unroll the generic loops, and P-256 is the hot curve. Every step +/// is the generic one's (`product`, `square_wide`, `redc`, `reduce`). +mod p256 { + use super::{Form, Modulus, ct, mac}; + use crate::nist::P256Field; + + const P: [u64; 4] = >::P; + const _: () = assert!(matches!(>::FORM, Form::P256)); + + #[inline(always)] + fn adc(a: u64, b: u64, carry: u64) -> (u64, u64) { + let (s, c1) = a.overflowing_add(b); + let (s, c2) = s.overflowing_add(carry); + (s, (c1 | c2) as u64) + } + + #[inline(always)] + fn sbb(a: u64, b: u64, borrow: u64) -> (u64, u64) { + let (d, b1) = a.overflowing_sub(b); + let (d, b2) = d.overflowing_sub(borrow); + (d, (b1 | b2) as u64) + } + + /// The four words as a residue's array, whose N is four here; indexing + /// rather than a slice copy keeps the words in registers. + #[inline(always)] + fn words(w: [u64; 4]) -> [u64; N] { + let mut out = [0; N]; + out[0] = w[0]; + out[1] = w[1]; + out[2] = w[2]; + out[3] = w[3]; + out + } + + /// v mod p for v + high 2^256 < 2p, p subtracted under a mask. + #[inline(always)] + fn reduce(v: [u64; 4], high: u64) -> [u64; N] { + let (d0, b) = sbb(v[0], P[0], 0); + let (d1, b) = sbb(v[1], P[1], b); + let (d2, b) = sbb(v[2], P[2], b); + let (d3, b) = sbb(v[3], P[3], b); + let m = ct::mask(ct::nonzero(high) | (1 ^ b)); + words([(v[0] & !m) | (d0 & m), (v[1] & !m) | (d1 & m), (v[2] & !m) | (d2 & m), (v[3] & !m) | (d3 & m)]) + } + + #[inline(always)] + pub(super) fn add(a: &[u64; N], b: &[u64; N]) -> [u64; N] { + let (v0, c) = adc(a[0], b[0], 0); + let (v1, c) = adc(a[1], b[1], c); + let (v2, c) = adc(a[2], b[2], c); + let (v3, c) = adc(a[3], b[3], c); + reduce([v0, v1, v2, v3], c) + } + + #[inline(always)] + pub(super) fn sub(a: &[u64; N], b: &[u64; N]) -> [u64; N] { + let (v0, w) = sbb(a[0], b[0], 0); + let (v1, w) = sbb(a[1], b[1], w); + let (v2, w) = sbb(a[2], b[2], w); + let (v3, w) = sbb(a[3], b[3], w); + let m = ct::mask(w); + let (r0, c) = adc(v0, P[0] & m, 0); + let (r1, c) = adc(v1, P[1] & m, c); + let (r2, c) = adc(v2, P[2] & m, c); + let (r3, _) = adc(v3, P[3] & m, c); + words([r0, r1, r2, r3]) + } + + /// CIOS with P-256's reduction step, as `product`. + #[inline(always)] + pub(super) fn mul(a: &[u64; N], b: &[u64; N]) -> [u64; N] { + let (mut t0, mut t1, mut t2, mut t3, mut high) = (0u64, 0u64, 0u64, 0u64, 0u64); + macro_rules! row { + ($ai:expr) => {{ + let ai = $ai; + let (x0, c) = mac(t0, ai, b[0], 0); + let (x1, c) = mac(t1, ai, b[1], c); + let (x2, c) = mac(t2, ai, b[2], c); + let (x3, c) = mac(t3, ai, b[3], c); + let (top, overflow) = adc(high, c, 0); + // q = x0; x0 + q (2^64 - 1) + q (2^32 - 1) 2^64 leaves q 2^32 + // at word one, and q p3 lands at word three. + let (y0, c) = adc(x1, x0 << 32, 0); + let (y1, c) = adc(x2, x0 >> 32, c); + let (y2, c) = mac(x3, x0, P[3], c); + let (y3, carry) = adc(top, c, 0); + (t0, t1, t2, t3, high) = (y0, y1, y2, y3, overflow | carry); + }}; + } + row!(a[0]); + row!(a[1]); + row!(a[2]); + row!(a[3]); + reduce([t0, t1, t2, t3], high) + } + + /// The cross products once, doubled, the squares added, then four + /// reduction steps, as `square_wide` and `redc`. + #[inline(always)] + pub(super) fn square(a: &[u64; N]) -> [u64; N] { + let (t1, c) = mac(0, a[0], a[1], 0); + let (t2, c) = mac(0, a[0], a[2], c); + let (t3, t4) = mac(0, a[0], a[3], c); + let (t3, c) = mac(t3, a[1], a[2], 0); + let (t4, t5) = mac(t4, a[1], a[3], c); + let (t5, t6) = mac(t5, a[2], a[3], 0); + let t7 = t6 >> 63; + let (t6, t5, t4, t3, t2, t1) = + (t6 << 1 | t5 >> 63, t5 << 1 | t4 >> 63, t4 << 1 | t3 >> 63, t3 << 1 | t2 >> 63, t2 << 1 | t1 >> 63, t1 << 1); + let (t0, hi) = mac(0, a[0], a[0], 0); + let (t1, c) = adc(t1, hi, 0); + let (lo, hi) = mac(0, a[1], a[1], 0); + let (t2, c) = adc(t2, lo, c); + let (t3, c) = adc(t3, hi, c); + let (lo, hi) = mac(0, a[2], a[2], 0); + let (t4, c) = adc(t4, lo, c); + let (t5, c) = adc(t5, hi, c); + let (lo, hi) = mac(0, a[3], a[3], 0); + let (t6, c) = adc(t6, lo, c); + let (t7, _) = adc(t7, hi, c); + let mut t = [t0, t1, t2, t3, t4, t5, t6, t7]; + let mut extra = 0u64; + macro_rules! step { + ($i:literal) => {{ + let q = t[$i]; + let (v, c) = adc(t[$i + 1], q << 32, 0); + t[$i + 1] = v; + let (v, c) = adc(t[$i + 2], q >> 32, c); + t[$i + 2] = v; + let (lo, hi) = mac(0, q, P[3], 0); + let (v, c) = adc(t[$i + 3], lo, c); + t[$i + 3] = v; + let (v, mut c) = adc(t[$i + 4], hi, c); + t[$i + 4] = v; + let mut k = $i + 5; + while k < 8 { + (t[k], c) = adc(t[k], c, 0); + k += 1; + } + extra |= c; + }}; + } + step!(0); + step!(1); + step!(2); + step!(3); + reduce([t[4], t[5], t[6], t[7]], extra) } } @@ -233,6 +521,123 @@ const fn product, const N: usize>(a: &[u64; N], b: &[u64; N]) -> ( (t, high) } +/// The 2N-word square of a, in a buffer sized for N <= 6. +#[inline(always)] +const fn square_wide(a: &[u64; N]) -> [u64; 12] { + let mut t = [0u64; 12]; + let mut i = 0; + while i < N { + let mut c = 0; + let mut j = i + 1; + while j < N { + (t[i + j], c) = mac(t[i + j], a[i], a[j], c); + j += 1; + } + t[i + N] = c; + i += 1; + } + // The cross products sum below 2^(128N - 1), so doubling loses no bit. + let mut top = 0; + let mut i = 0; + while i < 2 * N { + let next = t[i] >> 63; + t[i] = t[i] << 1 | top; + top = next; + i += 1; + } + let mut c = 0u64; + let mut i = 0; + while i < N { + let (lo, hi) = mac(0, a[i], a[i], 0); + let v = (t[2 * i] as u128).wrapping_add(lo as u128).wrapping_add(c as u128); + t[2 * i] = v as u64; + let v = (t[2 * i + 1] as u128).wrapping_add(hi as u128).wrapping_add(v >> 64); + t[2 * i + 1] = v as u64; + c = (v >> 64) as u64; + i += 1; + } + t +} + +/// t / R mod 2p for a 2N-word t < p^2: N reduction steps, each adding the +/// multiple of p that clears the lowest remaining word, then the top N words +/// and a top bit. +#[inline(always)] +const fn redc, const N: usize>(t: &mut [u64; 12]) -> ([u64; N], u64) { + // The total stays below 2^(128N + 1), so at most one step carries out of + // the top word. + let mut extra = 0u64; + let mut i = 0; + while i < N { + let mut c: u64; + let mut k; + match M::FORM { + Form::Generic => { + let q = t[i].wrapping_mul(M::N0); + c = 0; + let mut j = 0; + while j < N { + (t[i + j], c) = mac(t[i + j], q, M::P[j], c); + j += 1; + } + k = i + N; + } + Form::P256 => { + // t_i + q (2^64 - 1) = q 2^64 with q = t_i, and with q (2^32 - 1) + // at the next word that is q 2^32 there; q p3 lands two words up. + let q = t[i]; + t[i] = 0; + let v = (t[i + 1] as u128).wrapping_add((q << 32) as u128); + t[i + 1] = v as u64; + let v = (t[i + 2] as u128).wrapping_add((q >> 32) as u128).wrapping_add(v >> 64); + t[i + 2] = v as u64; + let (lo, hi) = mac(0, q, M::P[3], 0); + let v = (t[i + 3] as u128).wrapping_add(lo as u128).wrapping_add(v >> 64); + t[i + 3] = v as u64; + let v = (t[i + 4] as u128).wrapping_add(hi as u128).wrapping_add(v >> 64); + t[i + 4] = v as u64; + c = (v >> 64) as u64; + k = i + 5; + } + } + while k < 2 * N { + let (sum, carry) = t[k].overflowing_add(c); + t[k] = sum; + c = carry as u64; + k += 1; + } + extra |= c; + i += 1; + } + let mut out = [0; N]; + let mut j = 0; + while j < N { + out[j] = t[N + j]; + j += 1; + } + (out, extra) +} + +/// [`reduce`] for constants: the same arithmetic, without the register +/// barrier. +const fn reduce_const(v: &[u64; N], high: u64, p: &[u64; N]) -> [u64; N] { + let mut d = [0; N]; + let mut borrow = false; + let mut i = 0; + while i < N { + (d[i], borrow) = borrowing_sub(v[i], p[i], borrow); + i += 1; + } + let take = 0u64.wrapping_sub((high != 0 || !borrow) as u64); + let mut out = [0; N]; + let mut i = 0; + while i < N { + out[i] = (v[i] & !take) | (d[i] & take); + i += 1; + } + out +} + /// v mod p for v + high * 2^(64N) < 2p: p subtracted under a mask. fn reduce(v: &[u64; N], high: u64, p: &[u64; N]) -> [u64; N] { let mut d = [0; N]; @@ -259,16 +664,16 @@ pub(crate) const fn mac(t: u64, a: u64, b: u64, c: u64) -> (u64, u64) { } #[inline(always)] -fn carrying_add(a: u64, b: u64, carry: bool) -> (u64, bool) { +const fn carrying_add(a: u64, b: u64, carry: bool) -> (u64, bool) { let (s, c1) = a.overflowing_add(b); - let (s, c2) = s.overflowing_add(u64::from(carry)); + let (s, c2) = s.overflowing_add(carry as u64); (s, c1 | c2) } #[inline(always)] -fn borrowing_sub(a: u64, b: u64, borrow: bool) -> (u64, bool) { +const fn borrowing_sub(a: u64, b: u64, borrow: bool) -> (u64, bool) { let (d, b1) = a.overflowing_sub(b); - let (d, b2) = d.overflowing_sub(u64::from(borrow)); + let (d, b2) = d.overflowing_sub(borrow as u64); (d, b1 | b2) } @@ -321,3 +726,4 @@ const fn shift_mod(x: &[u64; N], p: &[u64; N], k: usize) -> [u64 } x } + diff --git a/lib/crypto/src/nist.rs b/lib/crypto/src/nist.rs index 3924148948d..78516578445 100644 --- a/lib/crypto/src/nist.rs +++ b/lib/crypto/src/nist.rs @@ -1,10 +1,12 @@ -//! SEC 1 NIST P-256/P-384 arithmetic. Secret scalar multiplication (ECDH -//! and key generation) runs Booth's signed 5-bit windows on a fixed schedule -//! with masked table reads and complete masked handling of exceptional -//! additions. SEC 1 decoding, DER decoding, ECDSA verification and -//! key-validity refusals depend on public data; verification multiplies in -//! variable time, through `PublicScalar` only. Only uncompressed SEC 1 -//! points are accepted. +//! SEC 1 NIST P-256/P-384 arithmetic. Secret scalar multiplication (ECDH and +//! key generation) runs Booth's signed 5-bit windows on a fixed schedule: +//! P-256's base point through a compile-time comb, other points through a +//! per-call table and doublings. Table reads scan every entry under masks, +//! and the scalar is first moved to an odd representative k + j n, so the +//! running point is never the identity. SEC 1 decoding, DER decoding, ECDSA +//! verification and key-validity refusals depend on public data; +//! verification multiplies in variable time, through `PublicScalar` only. +//! Only uncompressed SEC 1 points are accepted. #![allow( clippy::needless_range_loop, reason = "masked selects index parallel fixed-size arrays, as `montgomery` does" @@ -18,12 +20,23 @@ use crate::{ /// A short Weierstrass curve y^2 = x^3 - 3x + b of prime order and cofactor /// one, with its constants in Montgomery form at compile time. -pub(crate) trait Curve: 'static { +pub(crate) trait Curve: Sized + 'static { type Field: Modulus; type Order: Modulus; const B: [u64; N]; const GX: [u64; N]; const GY: [u64; N]; + /// m of the representative k + j n, j the one of m and m + 1 that makes + /// it odd, that secret multiplication runs on: the smallest m for which + /// every such representative sets a bit of the top Booth window, so the + /// first window's digit is nonzero. + const OFFSET: u64; + /// G, 3G, ..., 63G, for verification's width-7 recoding of u1. + fn odd_multiples_of_g() -> &'static [Affine; 32]; + /// scalar G, on a fixed schedule. + fn base_mul(scalar: &Scalar) -> Point { + Point::base().mul(scalar) + } } type Coord = Residue<>::Field, N>; type ModN = Residue<>::Order, N>; @@ -76,23 +89,46 @@ impl, const N: usize> Scalar { } Ok(scalar) } - /// Booth digit i in [-16, 16]: with u the six bits 5i - 1 ..= 5i + 4 - /// (bit -1 and bits past the top are zero), d = (u + 1) / 2 - 32 u_5, so - /// that the scalar is the sum of d_i 2^(5i). The index is public; the - /// arithmetic does not branch. - fn booth(&self, i: usize) -> i8 { - let bits = |from: usize| -> u64 { - // The six bits from `from`, a public position. - let (word, shift) = (from / 64, from % 64); - let low = if word < N { self.0[word] >> shift } else { 0 }; - let high = if shift > 58 && word + 1 < N { self.0[word + 1] << (64 - shift) } else { 0 }; - (low | high) & 63 - }; - let u = if i == 0 { bits(0) << 1 & 63 } else { bits(5 * i - 1) }; - (((u + 1) >> 1) as i8).wrapping_sub(((u >> 5) << 5) as i8) + /// k + j n for the one j of m = `C::OFFSET` and m + 1 that makes it odd: + /// an integer with the same multiple of every point, below 2^(64N + 2), in + /// N + 1 of the buffer's words. Erased when dropped. + fn representative(&self) -> Workspace<[u64; 7]> { + let multiplier = C::OFFSET.wrapping_add((C::OFFSET ^ self.0[0] ^ 1) & 1); + let mut k = Workspace([0u64; 7]); + let mut carry = 0u64; + for i in 0..N { + let (lo, hi) = montgomery::mac(self.0[i], C::Order::P[i], multiplier, carry); + k.0[i] = lo; + carry = hi; + } + k.0[N] = carry; + k } } +/// Booth digit i of a W-word integer, in [-16, 16]: with u the six bits +/// 5i - 1 ..= 5i + 4 (bit -1 and bits past the top are zero), +/// d = (u + 1) / 2 - 32 u_5, so that the integer is the sum of d_i 2^(5i) +/// while its bit 5i + 4 of the top window is clear. The index is public; the +/// arithmetic does not branch. +fn booth(words: &[u64; W], i: usize) -> i8 { + let bits = |from: usize| -> u64 { + let (word, shift) = (from / 64, from % 64); + let low = if word < W { words[word] >> shift } else { 0 }; + let high = if shift > 58 && word + 1 < W { words[word + 1] << (64 - shift) } else { 0 }; + (low | high) & 63 + }; + let u = if i == 0 { bits(0) << 1 & 63 } else { bits(5 * i - 1) }; + ((u.wrapping_add(1) >> 1) as i8).wrapping_sub(((u >> 5) << 5) as i8) +} + +/// |d| and the sign bit of a Booth digit, by arithmetic. +fn magnitude(d: i8) -> (u64, u64) { + let negative = u64::from((d as u8) >> 7); + let magnitude = ((d ^ 0i8.wrapping_sub(negative as i8)).wrapping_add(negative as i8)) as u64; + (magnitude, negative) +} + /// A public scalar below n, for variable-time verification only; a secret /// [`Scalar`] does not convert to it. struct PublicScalar([u64; N]); @@ -106,12 +142,189 @@ impl PublicScalar { } } -/// The windows of [`Point::mul`]: 5-bit Booth windows over 64N bits, plus -/// the one that absorbs the top window's carry. +/// The 5-bit Booth windows of a representative below 2^(64N + 2), whose top +/// window keeps its bit 5i + 4 clear. const fn windows(n: usize) -> usize { - (64 * n) / 5 + 1 + (64 * n + 2) / 5 + 1 } +/// An affine point (x, y), a table entry. +pub(crate) struct Affine, const N: usize> { + x: Coord, + y: Coord, +} +impl, const N: usize> Clone for Affine { + fn clone(&self) -> Self { + *self + } +} +impl, const N: usize> Copy for Affine {} + +impl, const N: usize> Affine { + /// row[|d| - 1], with y negated by mask when d < 0, reading every entry. + /// For d = 0 it is row[0], which the caller discards by mask: no entry + /// is the identity, so no arithmetic runs on it. + fn lookup(row: &[Self; 16], d: i8) -> Self { + let (magnitude, negative) = magnitude(d); + let mut chosen = row[0]; + // Indexing, not an iterator adaptor, which the shipped profile keeps + // out of line. + for j in 1..16 { + let take = 1 ^ crate::ct::nonzero(magnitude ^ (j as u64 + 1)); + chosen.x = Residue::select(&chosen.x, &row[j].x, take); + chosen.y = Residue::select(&chosen.y, &row[j].y, take); + } + chosen.y = Residue::select(&chosen.y, &Residue::zero().sub(&chosen.y), negative); + chosen + } + fn point(&self) -> Point { + Point { + x: self.x, + y: self.y, + z: Residue::one(), + } + } +} + +/// Compile-time point arithmetic, in Jacobian (X, Y, Z), for building +/// tables. None of its inputs is the identity, and no addition's operands +/// are equal or opposite. +struct Build, const N: usize>(core::marker::PhantomData); + +impl, const N: usize> Build { + const G: [Coord; 3] = [ + Residue::from_montgomery(C::GX), + Residue::from_montgomery(C::GY), + Residue::one(), + ]; + /// dbl-2001-b for a = -3, as [`Point::double`]. + const fn double(p: &[Coord; 3]) -> [Coord; 3] { + let [x, y, z] = p; + let delta = z.times(z); + let gamma = y.times(y); + let beta = x.times(&gamma); + let t = x.minus(&delta).times(&x.plus(&delta)); + let alpha = t.plus(&t).plus(&t); + let beta4 = beta.plus(&beta).plus(&beta.plus(&beta)); + let x3 = alpha.times(&alpha).minus(&beta4.plus(&beta4)); + let z3 = y.plus(z).times(&y.plus(z)).minus(&gamma).minus(&delta); + let gamma2 = gamma.times(&gamma); + let gamma8 = gamma2.plus(&gamma2).plus(&gamma2.plus(&gamma2)); + let y3 = alpha.times(&beta4.minus(&x3)).minus(&gamma8.plus(&gamma8)); + [x3, y3, z3] + } + /// The Jacobian addition of [`Point::add_incomplete`]. + const fn add(p: &[Coord; 3], q: &[Coord; 3]) -> [Coord; 3] { + let (z1, z2) = (p[2].times(&p[2]), q[2].times(&q[2])); + let (u1, u2) = (p[0].times(&z2), q[0].times(&z1)); + let s1 = p[1].times(&q[2].times(&z2)); + let s2 = q[1].times(&p[2].times(&z1)); + let h = u2.minus(&u1); + let r = s2.minus(&s1); + let h2 = h.times(&h); + let h3 = h.times(&h2); + let u = u1.times(&h2); + let x = r.times(&r).minus(&h3).minus(&u.plus(&u)); + let y = r.times(&u.minus(&x)).minus(&s1.times(&h3)); + let z = h.times(&p[2].times(&q[2])); + [x, y, z] + } + /// Jacobian points in affine form, with one inversion for all of them + /// (Montgomery's simultaneous-inversion trick). + const fn normalize(points: &[[Coord; 3]; K]) -> [Affine; K] { + let mut prefix = [Residue::one(); K]; + prefix[0] = points[0][2]; + let mut i = 1; + while i < K { + prefix[i] = prefix[i - 1].times(&points[i][2]); + i += 1; + } + let mut inverse = prefix[K - 1].inverted(); + let mut out = [Affine { + x: Residue::zero(), + y: Residue::zero(), + }; K]; + let mut i = K; + while i > 0 { + i -= 1; + let z_inverse = if i == 0 { inverse } else { inverse.times(&prefix[i - 1]) }; + inverse = inverse.times(&points[i][2]); + let z2 = z_inverse.times(&z_inverse); + out[i] = Affine { + x: points[i][0].times(&z2), + y: points[i][1].times(&z2.times(&z_inverse)), + }; + } + out + } + /// G, 3G, ..., (2K - 1)G. + const fn odd_multiples() -> [Affine; K] { + let twice = Self::double(&Self::G); + let mut points = [Self::G; K]; + let mut i = 1; + while i < K { + points[i] = Self::add(&points[i - 1], &twice); + i += 1; + } + Self::normalize(&points) + } + /// Row i holds j 2^(5i) G for j = 1..=16. + const fn comb() -> [[Affine; 16]; R] { + let mut points = [[Self::G; 16]; R]; + let mut base = Self::G; + let mut i = 0; + while i < R { + points[i][0] = base; + points[i][1] = Self::double(&base); + let mut j = 2; + while j < 16 { + points[i][j] = Self::add(&points[i][j - 1], &base); + j += 1; + } + let mut k = 0; + while k < 5 { + base = Self::double(&base); + k += 1; + } + i += 1; + } + // One simultaneous inversion for every entry, in row order. + let mut prefix = [[Residue::one(); 16]; R]; + let mut running = Residue::one(); + let mut i = 0; + while i < 16 * R { + running = running.times(&points[i / 16][i % 16][2]); + prefix[i / 16][i % 16] = running; + i += 1; + } + let mut inverse = running.inverted(); + let mut out = [[Affine { + x: Residue::zero(), + y: Residue::zero(), + }; 16]; R]; + let mut i = 16 * R; + while i > 0 { + i -= 1; + let p = &points[i / 16][i % 16]; + let z_inverse = if i == 0 { inverse } else { inverse.times(&prefix[(i - 1) / 16][(i - 1) % 16]) }; + inverse = inverse.times(&p[2]); + let z2 = z_inverse.times(&z_inverse); + out[i / 16][i % 16] = Affine { + x: p[0].times(&z2), + y: p[1].times(&z2.times(&z_inverse)), + }; + } + out + } +} + +/// P-256's comb: j 2^(5i) G for j = 1..=16 in 52 rows (53 KiB). +// Building it runs some ten thousand field multiplications at compile time. +#[allow(long_running_const_eval)] +static P256_COMB: [[Affine; 16]; windows(4)] = Build::::comb(); +static P256_G_ODD: [Affine; 32] = Build::::odd_multiples(); +static P384_G_ODD: [Affine; 32] = Build::::odd_multiples(); + /// The odd multiples P, 3P, ..., 15P of a public point, for verification. type OddMultiples = [Point; 8]; @@ -168,7 +381,46 @@ impl, const N: usize> Point { let y = alpha.mul(&beta4.sub(&x)).sub(&gamma2.double().double()); Self { x, y, z } } - fn add(&self, q: &Self) -> Self { + /// scalar * self on a fixed schedule, for a point of prime order n: + /// Booth's signed 5-bit windows of the odd representative k' over the + /// public multiples P, 2P, ..., 16P, each window five doublings and one + /// addition of a masked table read. self is public; neither branches nor + /// addresses depend on the scalar. + /// + /// Exceptional additions. Before the addition of window i the running + /// point is 32 T P, for T the value of the digits above i, and + /// k' = 2^(5i + 5) T + S with |S| <= 2^(5i + 4). As k' > 2^(5 top - 1) + /// (`Curve::OFFSET`), T >= 1 for every i < top; and as k' < 4n, for i >= 1 + /// 32 T <= k' / 2^(5i) + 16 < 4n / 32 + 16. So 0 < 32 T -+ d < n for + /// every |d| <= 16: the running point is not the identity and neither + /// equals nor opposes dP, the addend or the discarded table[0]. Only the + /// last addition, i = 0, can double, and it does so by mask; it cannot + /// cancel, since k' is not 0 mod n. + pub(crate) fn mul(&self, scalar: &Scalar) -> Self { + let mut table = [*self; 16]; + table[1] = self.double(); + for j in 2..16 { + table[j] = table[j - 1].add_public(self); + } + let k = scalar.representative(); + let top = windows(N) - 1; + let mut r = Self::lookup(&table, booth(&k.0, top)); + for i in (0..top).rev() { + for _ in 0..5 { + r = r.double(); + } + let d = booth(&k.0, i); + let (mut sum, equal) = r.add_incomplete(&Self::lookup(&table, d)); + if i == 0 { + sum = Self::select(&sum, &r.double(), equal); + } + r = Self::select(&r, &sum, crate::ct::nonzero(d as u8 as u64)); + } + r + } + /// self + q for points that are not the identity, 12M + 4S, and one when + /// they are equal (h = r = 0), where the formula fails. + fn add_incomplete(&self, q: &Self) -> (Self, u64) { let z1 = self.z.square(); let z2 = q.z.square(); let u1 = self.x.mul(&z2); @@ -183,44 +435,75 @@ impl, const N: usize> Point { let x = r.square().sub(&h3).sub(&u.double()); let y = r.mul(&u.sub(&x)).sub(&s1.mul(&h3)); let z = h.mul(&self.z.mul(&q.z)); - let out = Self { x, y, z }; - let out = Self::select(&out, &self.double(), h.is_zero() & r.is_zero()); - let out = Self::select(&out, self, q.z.is_zero()); - Self::select(&out, q, self.z.is_zero()) - } - /// scalar * self on a fixed schedule: Booth's signed 5-bit windows over - /// the multiples P, 2P, ..., 16P, each window five doublings and one - /// complete addition of a masked table read. self is public; neither - /// branches nor addresses depend on the scalar. - pub(crate) fn mul(&self, scalar: &Scalar) -> Self { - let mut table = [*self; 16]; - table[1] = self.double(); - for j in 2..16 { - table[j] = table[j - 1].add(self); - } - let top = windows(N) - 1; - let mut r = Self::lookup(&table, scalar.booth(top)); - for i in (0..top).rev() { - for _ in 0..5 { - r = r.double(); - } - r = r.add(&Self::lookup(&table, scalar.booth(i))); - } - r + (Self { x, y, z }, h.is_zero() & r.is_zero()) } - /// d P from table[j] = (j + 1) P, for d in [-16, 16], reading every - /// entry: the identity for zero, and Y negated by mask for d < 0. + /// table[|d| - 1] for table[j] = (j + 1) P, with Y negated by mask for + /// d < 0, reading every entry. For d = 0 it is table[0], which the caller + /// discards by mask. fn lookup(table: &[Self; 16], d: i8) -> Self { - let negative = u64::from((d as u8) >> 7); - let magnitude = ((d ^ 0i8.wrapping_sub(negative as i8)).wrapping_add(negative as i8)) as u64; - let mut chosen = Self::identity(); - for (j, entry) in table.iter().enumerate() { - chosen = Self::select(&chosen, entry, 1 ^ crate::ct::nonzero(magnitude ^ (j as u64 + 1))); + let (magnitude, negative) = magnitude(d); + let mut chosen = table[0]; + for j in 1..16 { + chosen = Self::select(&chosen, &table[j], 1 ^ crate::ct::nonzero(magnitude ^ (j as u64 + 1))); } let negated = Residue::zero().sub(&chosen.y); chosen.y = Residue::select(&chosen.y, &negated, negative); chosen } + /// self + q for an affine q, 8M + 3S, and one when the operands are + /// equal (h = r = 0), where the formula fails. self is not the identity. + fn add_affine(&self, q: &Affine) -> (Self, u64) { + let z1 = self.z.square(); + let u2 = q.x.mul(&z1); + let s2 = q.y.mul(&self.z.mul(&z1)); + let h = u2.sub(&self.x); + let r = s2.sub(&self.y); + let h2 = h.square(); + let h3 = h.mul(&h2); + let u = self.x.mul(&h2); + let x = r.square().sub(&h3).sub(&u.double()); + let y = r.mul(&u.sub(&x)).sub(&self.y.mul(&h3)); + let z = h.mul(&self.z); + (Self { x, y, z }, h.is_zero() & r.is_zero()) + } + /// scalar G by a comb, with no doublings: the odd representative's Booth + /// digit i selects from row i, j 2^(5i) G, and one mixed addition adds + /// it. + /// + /// Exceptional additions. Before window i >= 1 the running point is + /// S_i G for the odd S_i = k' mod 2^(5i) - b 2^(5i), b = bit 5i - 1 of + /// the representative k', so 0 < |S_i| <= 2^(5i - 1): never the identity. + /// It equals or opposes the addend d 2^(5i) G only if + /// S_i -+ d 2^(5i) = 0 mod n, while 0 < |S_i -+ d 2^(5i)| <= 16.5 2^(5i), + /// which is below n for i < R - 1. At the top window S + d 2^(5i) = k' is + /// not 0 mod n, but S = d 2^(5i) can hold, for k = 2^256 d mod n; that one + /// addition doubles by mask instead. A zero digit's addition is + /// discarded by mask. + fn comb_mul(table: &[[Affine; 16]; R], scalar: &Scalar) -> Self { + let k = scalar.representative(); + // k' is odd, so digit 0 is nonzero. + let mut r = Affine::lookup(&table[0], booth(&k.0, 0)).point(); + for i in 1..R { + let d = booth(&k.0, i); + let (mut sum, equal) = r.add_affine(&Affine::lookup(&table[i], d)); + if i == R - 1 { + sum = Self::select(&sum, &r.double(), equal); + } + r = Self::select(&r, &sum, crate::ct::nonzero(d as u8 as u64)); + } + r + } + /// self + q for a public affine q, branching on the values. + fn add_affine_public(&self, q: &Affine) -> Self { + if self.z.is_zero() == 1 { + return q.point(); + } + let (sum, equal) = self.add_affine(q); + if sum.z.is_zero() == 1 { + return if equal == 1 { self.double() } else { Self::identity() }; + } + sum + } /// self + q for public points, branching on their values: no masked /// doubling or selects. fn add_public(&self, q: &Self) -> Self { @@ -230,24 +513,11 @@ impl, const N: usize> Point { if q.z.is_zero() == 1 { return *self; } - let z1 = self.z.square(); - let z2 = q.z.square(); - let u1 = self.x.mul(&z2); - let u2 = q.x.mul(&z1); - let s1 = self.y.mul(&q.z.mul(&z2)); - let s2 = q.y.mul(&self.z.mul(&z1)); - let h = u2.sub(&u1); - let r = s2.sub(&s1); - if h.is_zero() == 1 { - return if r.is_zero() == 1 { self.double() } else { Self::identity() }; + let (sum, equal) = self.add_incomplete(q); + if sum.z.is_zero() == 1 { + return if equal == 1 { self.double() } else { Self::identity() }; } - let h2 = h.square(); - let h3 = h.mul(&h2); - let u = u1.mul(&h2); - let x = r.square().sub(&h3).sub(&u.double()); - let y = r.mul(&u.sub(&x)).sub(&s1.mul(&h3)); - let z = h.mul(&self.z.mul(&q.z)); - Self { x, y, z } + sum } fn neg(&self) -> Self { Self { @@ -267,20 +537,29 @@ impl, const N: usize> Point { /// u1 G + u2 self in one doubling chain over both scalars' width-5 NAFs /// (Straus's method), in variable time: every input is public. fn joint(&self, u1: &PublicScalar, u2: &PublicScalar) -> Self { - let (g, q) = (Self::base().odd_multiples(), self.odd_multiples()); - let (d1, d2) = (u1.naf(5), u2.naf(5)); + let (g, q) = (C::odd_multiples_of_g(), self.odd_multiples()); + let (d1, d2) = (u1.naf(7), u2.naf(5)); let Some(top) = (0..64 * N + 1).rev().find(|&i| d1[i] != 0 || d2[i] != 0) else { return Self::identity(); }; let mut r = Self::identity(); for i in (0..=top).rev() { r = r.double(); - for (d, table) in [(d1[i], &g), (d2[i], &q)] { - if d > 0 { - r = r.add_public(&table[d as usize / 2]); - } else if d < 0 { - r = r.add_public(&table[(-d) as usize / 2].neg()); - } + let d = d1[i]; + if d > 0 { + r = r.add_affine_public(&g[d as usize / 2]); + } else if d < 0 { + let entry = g[(-d) as usize / 2]; + r = r.add_affine_public(&Affine { + x: entry.x, + y: Residue::zero().sub(&entry.y), + }); + } + let d = d2[i]; + if d > 0 { + r = r.add_public(&q[d as usize / 2]); + } else if d < 0 { + r = r.add_public(&q[(-d) as usize / 2].neg()); } } r @@ -471,6 +750,14 @@ pub(crate) struct P256; impl Curve<4> for P256 { type Field = P256Field; type Order = P256Order; + // k + n >= n > 2^255 sets bit 254 or 255 of window 51 (bits 254..=259). + const OFFSET: u64 = 1; + fn odd_multiples_of_g() -> &'static [Affine; 32] { + &P256_G_ODD + } + fn base_mul(scalar: &Scalar) -> Point { + Point::comb_mul(&P256_COMB, scalar) + } const B: [u64; 4] = montgomery::encode( &[0x3bce3c3e27d2604b, 0x651d06b0cc53b0f6, 0xb3ebbd55769886bc, 0x5ac635d8aa3a93e7], &P256Field::P, @@ -490,6 +777,12 @@ pub(crate) struct P384; impl Curve<6> for P384 { type Field = P384Field; type Order = P384Order; + // k + 2n >= 2n > 2^384 sets bit 384 or 385 of window 77 (bits 384..=389); + // k + n need not. + const OFFSET: u64 = 2; + fn odd_multiples_of_g() -> &'static [Affine; 32] { + &P384_G_ODD + } const B: [u64; 6] = montgomery::encode( &[ 0x2a85c8edd3ec2aef, @@ -572,22 +865,31 @@ mod tests { ) } - /// Booth digits lie in [-16, 16] and recompose the scalar, for both word - /// counts and at the patterns that give digits of +16 and -16. + /// The representative is odd, congruent to k and sets the top window; + /// its Booth digits lie in [-16, 16] and recompose it, with +16 and -16 + /// both reached. fn check_booth, const N: usize>() { use num_bigint::BigInt; + let n = integer(&C::Order::P); let mut extremes = [false; 2]; for k in scalars(&C::Order::P) { + let k = ModN::::from_wide(&k).integer(); let scalar = Scalar::(k, core::marker::PhantomData); + let words = scalar.representative(); + let representative = integer(&words.0); + assert_eq!(&representative % 2u8, BigInt::from(1)); + assert_eq!(&representative % &n, integer(&k)); + assert!(representative < BigInt::from(1) << (64 * N + 2)); let mut sum = BigInt::from(0); for i in 0..windows(N) { - let d = scalar.booth(i); + let d = booth(&words.0, i); assert!((-16..=16).contains(&d)); extremes[0] |= d == 16; extremes[1] |= d == -16; sum += BigInt::from(d) << (5 * i); } - assert_eq!(sum, integer(&k)); + assert_eq!(sum, representative); + assert_ne!(booth(&words.0, windows(N) - 1), 0); } assert_eq!(extremes, [true, true]); } @@ -651,6 +953,99 @@ mod tests { check_exceptional_sums::(); } + fn affine_of, const N: usize>(p: &Point) -> ([u64; N], [u64; N]) { + let (mut x, mut y) = ([0; N], [0; N]); + p.affine(&mut x, &mut y).unwrap(); + (x, y) + } + + /// k G by the public path, for k below n. + fn reference, const N: usize>(k: &num_bigint::BigInt) -> ([u64; N], [u64; N]) { + let (_, bytes) = k.to_bytes_le(); + let mut words = [0u64; N]; + for (i, chunk) in bytes.chunks(8).enumerate() { + let mut b = [0u8; 8]; + b[..chunk.len()].copy_from_slice(chunk); + words[i] = u64::from_le_bytes(b); + } + affine_of(&Point::::identity().joint(&PublicScalar(words), &PublicScalar([0; N]))) + } + + /// Every comb entry is j 2^(5i) G, and the odd-multiple tables (2j + 1) G. + #[test] + fn static_tables_hold_their_multiples() { + use num_bigint::BigInt; + let n = integer(&P256Order::P); + for (i, row) in P256_COMB.iter().enumerate() { + for (j, entry) in row.iter().enumerate() { + let k = (BigInt::from(j + 1) << (5 * i)) % &n; + assert_eq!(affine_of(&entry.point()), reference::(&k), "row {i}, entry {j}"); + } + } + for (j, entry) in P256_G_ODD.iter().enumerate() { + assert_eq!(affine_of(&entry.point()), reference::(&BigInt::from(2 * j + 1))); + } + for (j, entry) in P384_G_ODD.iter().enumerate() { + assert_eq!(affine_of(&entry.point()), reference::(&BigInt::from(2 * j + 1))); + } + } + + /// The comb equals the public path for every test scalar and for the + /// ones whose top addition doubles, k = 2^256 d mod n with d = 1, 2. + #[test] + fn the_comb_matches_the_public_path() { + use num_bigint::BigInt; + let n = integer(&P256Order::P); + let mut ks: std::vec::Vec = + scalars(&P256Order::P).iter().map(|k| integer(k) % &n).filter(|k| *k != BigInt::from(0)).collect(); + for d in 1..=4u8 { + ks.push((BigInt::from(d) << 256) % &n); + } + for k in ks { + let (_, bytes) = k.to_bytes_le(); + let mut words = [0u64; 4]; + for (i, chunk) in bytes.chunks(8).enumerate() { + let mut b = [0u8; 8]; + b[..chunk.len()].copy_from_slice(chunk); + words[i] = u64::from_le_bytes(b); + } + let comb = P256::base_mul(&Scalar(words, core::marker::PhantomData)); + assert_eq!(affine_of(&comb), reference::(&k), "k = {k:x}"); + } + } + + /// The scalars k = 2d mod n whose representative's first digit is d, + /// for which the last addition of [`Point::mul`] doubles: at least one + /// exists per curve, and each multiplies as the public path does. + fn check_last_addition_doubles, const N: usize>() { + use num_bigint::BigInt; + let n = integer(&C::Order::P); + let mut found = 0; + for d in (-16i32..=16).filter(|&d| d != 0) { + let k = ((BigInt::from(2 * d) % &n) + &n) % &n; + let (_, bytes) = k.to_bytes_le(); + let mut words = [0u64; N]; + for (i, chunk) in bytes.chunks(8).enumerate() { + let mut b = [0u8; 8]; + b[..chunk.len()].copy_from_slice(chunk); + words[i] = u64::from_le_bytes(b); + } + let scalar = Scalar::(words, core::marker::PhantomData); + if i32::from(booth(&scalar.representative().0, 0)) != d { + continue; + } + found += 1; + assert_eq!(affine_of(&Point::::base().mul(&scalar)), reference::(&k), "d = {d}"); + } + assert!(found > 0); + } + + #[test] + fn the_last_addition_doubles_where_it_must() { + check_last_addition_doubles::(); + check_last_addition_doubles::(); + } + #[test] fn booth_multiplication_matches_the_public_path() { check_mul::(); diff --git a/lib/crypto/src/test_arithmetic.rs b/lib/crypto/src/test_arithmetic.rs index 310012f7ec1..b4f5e379ee8 100644 --- a/lib/crypto/src/test_arithmetic.rs +++ b/lib/crypto/src/test_arithmetic.rs @@ -129,6 +129,42 @@ impl, const N: usize> Field for Montgomery { } } +/// The same moduli through the `const` forms that build compile-time tables. +struct MontgomeryConst(PhantomData); + +impl, const N: usize> Field for MontgomeryConst { + type E = Residue; + fn modulus(&self) -> BigUint { + integer(&M::P) + } + fn input_bound(&self) -> BigUint { + BigUint::from(1u8) << (64 * N) + } + fn from(&self, value: &BigUint) -> Residue { + Residue::from_wide(&words(value)) + } + fn value(&self, e: &Residue) -> BigUint { + let value = integer(&e.integer()); + assert!(value < self.modulus(), "a residue at or above the modulus"); + value + } + fn add(&self, a: &Residue, b: &Residue) -> Residue { + a.plus(b) + } + fn sub(&self, a: &Residue, b: &Residue) -> Residue { + a.minus(b) + } + fn mul(&self, a: &Residue, b: &Residue) -> Residue { + a.times(b) + } + fn square(&self, a: &Residue) -> Residue { + a.times(a) + } + fn invert(&self, a: &Residue) -> Residue { + a.inverted() + } +} + struct XorShift(u64); impl XorShift { @@ -305,6 +341,14 @@ fn nist_fields_and_orders_match_integers() { check(&Montgomery::(PhantomData), 4); } +#[test] +fn compile_time_montgomery_forms_match_integers() { + check(&MontgomeryConst::(PhantomData), 5); + check(&MontgomeryConst::(PhantomData), 6); + check(&MontgomeryConst::(PhantomData), 7); + check(&MontgomeryConst::(PhantomData), 8); +} + /// The compile-time Montgomery constants, and the canonical decoding that /// refuses every value from the modulus up. #[test] From 0193e43e3e9cf7935db540d7224d8868cc983980 Mon Sep 17 00:00:00 2001 From: japabu Date: Thu, 22 Oct 2026 15:25:03 +0200 Subject: [PATCH 2/2] crypto: review round 1 of the P-256 comb - The variable-base proof covers i = 0: 32T = k' - d0 is a nonzero multiple of 32 below 4n + 16 and n is odd, so the running point is not the identity. - The comb proof says why the difference is nonzero: S_i is odd and d 2^(5i) even. - The comb is 52 KiB; the comb test's doc names d = 1..=4. - The const forms are documented for compile time and tests only, not for secrets: they lack the register barrier and `minus` branches. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01PXEk8BBMDT8tXu4AnENB9Q --- lib/crypto/src/montgomery.rs | 9 +++++---- lib/crypto/src/nist.rs | 17 ++++++++++------- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/lib/crypto/src/montgomery.rs b/lib/crypto/src/montgomery.rs index ee94471d20a..abc38e0de34 100644 --- a/lib/crypto/src/montgomery.rs +++ b/lib/crypto/src/montgomery.rs @@ -136,8 +136,8 @@ impl, const N: usize> Residue { Self(reduce(&t, high, &M::P), PhantomData) } /// The `const` forms of [`Self::mul`], [`Self::add`] and [`Self::sub`], - /// for tables built at compile time; their final subtraction is a plain - /// select. + /// for tables built at compile time and for tests: not for secrets. They + /// lack the register barrier, and `minus` branches on its borrow. pub(crate) const fn times(&self, b: &Self) -> Self { let (t, high) = product::(&self.0, &b.0); Self(reduce_const(&t, high, &M::P), PhantomData) @@ -201,7 +201,7 @@ impl, const N: usize> Residue { } r } - /// The `const` form of [`Self::square_times`]. + /// The `const` form of [`Self::square_times`], as [`Self::times`]. const fn squared_times(&self, k: u32) -> Self { let mut r = *self; let mut i = 0; @@ -236,7 +236,8 @@ impl, const N: usize> Residue { } } } - /// The `const` form of [`Self::invert`]. + /// The `const` form of [`Self::invert`], for compile time and tests, as + /// [`Self::times`]. pub(crate) const fn inverted(&self) -> Self { match M::FORM { Form::P256 => inversion_chain!(self, squared_times, times), diff --git a/lib/crypto/src/nist.rs b/lib/crypto/src/nist.rs index 78516578445..0d379f0138b 100644 --- a/lib/crypto/src/nist.rs +++ b/lib/crypto/src/nist.rs @@ -318,7 +318,7 @@ impl, const N: usize> Build { } } -/// P-256's comb: j 2^(5i) G for j = 1..=16 in 52 rows (53 KiB). +/// P-256's comb: j 2^(5i) G for j = 1..=16 in 52 rows (52 KiB). // Building it runs some ten thousand field multiplications at compile time. #[allow(long_running_const_eval)] static P256_COMB: [[Affine; 16]; windows(4)] = Build::::comb(); @@ -393,9 +393,11 @@ impl, const N: usize> Point { /// (`Curve::OFFSET`), T >= 1 for every i < top; and as k' < 4n, for i >= 1 /// 32 T <= k' / 2^(5i) + 16 < 4n / 32 + 16. So 0 < 32 T -+ d < n for /// every |d| <= 16: the running point is not the identity and neither - /// equals nor opposes dP, the addend or the discarded table[0]. Only the - /// last addition, i = 0, can double, and it does so by mask; it cannot - /// cancel, since k' is not 0 mod n. + /// equals nor opposes dP, the addend or the discarded table[0]. At i = 0, + /// 32 T = k' - d0 is a nonzero multiple of 32 below 4n + 16, and n is odd, + /// so it is not a multiple of n: the running point is not the identity. + /// That last addition can double, and does so by mask; it cannot cancel, + /// since k' is not 0 mod n. pub(crate) fn mul(&self, scalar: &Scalar) -> Self { let mut table = [*self; 16]; table[1] = self.double(); @@ -474,8 +476,8 @@ impl, const N: usize> Point { /// S_i G for the odd S_i = k' mod 2^(5i) - b 2^(5i), b = bit 5i - 1 of /// the representative k', so 0 < |S_i| <= 2^(5i - 1): never the identity. /// It equals or opposes the addend d 2^(5i) G only if - /// S_i -+ d 2^(5i) = 0 mod n, while 0 < |S_i -+ d 2^(5i)| <= 16.5 2^(5i), - /// which is below n for i < R - 1. At the top window S + d 2^(5i) = k' is + /// S_i -+ d 2^(5i) = 0 mod n, while S_i is odd and d 2^(5i) even, so + /// 0 < |S_i -+ d 2^(5i)| <= 16.5 2^(5i), which is below n for i < R - 1. At the top window S + d 2^(5i) = k' is /// not 0 mod n, but S = d 2^(5i) can hold, for k = 2^256 d mod n; that one /// addition doubles by mask instead. A zero digit's addition is /// discarded by mask. @@ -991,7 +993,8 @@ mod tests { } /// The comb equals the public path for every test scalar and for the - /// ones whose top addition doubles, k = 2^256 d mod n with d = 1, 2. + /// ones k = 2^256 d mod n for d = 1..=4, among them those whose top + /// addition doubles. #[test] fn the_comb_matches_the_public_path() { use num_bigint::BigInt;