diff --git a/lib/crypto/src/keccak.rs b/lib/crypto/src/keccak.rs index f88c32385c..95b4c652af 100644 --- a/lib/crypto/src/keccak.rs +++ b/lib/crypto/src/keccak.rs @@ -1,6 +1,6 @@ -//! FIPS 202 sponge used only by ML-KEM. Absorption and permutation have fixed -//! schedules for public lengths. Rate/domain separation are internal constants. -//! State and explicit permutation workspaces are erased, including on failure. +//! FIPS 202 sponge used only by ML-KEM. Its schedule depends only on the rate +//! and the input and output lengths, which are public; the domain and rate are +//! internal constants. The state is erased when the sponge drops. use crate::secret::wipe64; const RC: [u64; 24] = [ 1, @@ -28,7 +28,8 @@ const RC: [u64; 24] = [ 0x80000001, 0x8000000080008008, ]; -const ROT: [u32; 25] = [ +/// ρ's rotation of lane x + 5y (FIPS 202 §3.2.2). +const RHO: [u32; 25] = [ 0, 1, 62, 28, 27, 36, 44, 6, 55, 20, 3, 10, 43, 25, 39, 41, 45, 15, 21, 8, 18, 2, 61, 56, 14, ]; pub(crate) struct Sponge { @@ -49,61 +50,80 @@ impl Sponge { rate, position: 0, }; - for p in parts { - for &b in *p { - s.state[s.position / 8] ^= u64::from(b) << (8 * (s.position % 8)); - s.position += 1; - if s.position == rate { - s.permute(); - s.position = 0 - } - } + for part in parts { + s.absorb(part) } - s.state[s.position / 8] ^= u64::from(domain) << (8 * (s.position % 8)); - let last = rate - 1; - s.state[last / 8] ^= 0x80u64 << (8 * (last % 8)); - s.permute(); + s.xor_byte(s.position, domain); + s.xor_byte(rate - 1, 0x80); + permute(&mut s.state); s.position = 0; s } - fn permute(&mut self) { - let mut c = [0; 5]; - let mut d = [0; 5]; - let mut b = [0; 25]; - for rc in RC { - for x in 0..5 { - c[x] = (0..5).fold(0, |c, y| c ^ self.state[x + 5 * y]) - } - for x in 0..5 { - d[x] = c[(x + 4) % 5] ^ c[(x + 1) % 5].rotate_left(1) - } - for y in 0..5 { - for x in 0..5 { - let i = x + 5 * y; - b[y + 5 * ((2 * x + 3 * y) % 5)] = (self.state[i] ^ d[x]).rotate_left(ROT[i]); + fn xor_byte(&mut self, at: usize, byte: u8) { + self.state[at / 8] ^= u64::from(byte) << (8 * (at % 8)) + } + /// A whole lane at a time where the position is on one; every rate is. + fn absorb(&mut self, mut data: &[u8]) { + while !data.is_empty() { + let at = self.position; + let step = match data.first_chunk::<8>() { + Some(word) if at.is_multiple_of(8) => { + self.state[at / 8] ^= u64::from_le_bytes(*word); + 8 } - } - for y in 0..5 { - for x in 0..5 { - self.state[x + 5 * y] = - b[x + 5 * y] ^ (!b[(x + 1) % 5 + 5 * y] & b[(x + 2) % 5 + 5 * y]); + _ => { + self.xor_byte(at, data[0]); + 1 } + }; + data = &data[step..]; + self.position += step; + if self.position == self.rate { + permute(&mut self.state); + self.position = 0 } - self.state[0] ^= rc; } - wipe64(&mut c); - wipe64(&mut d); - wipe64(&mut b) } - pub(crate) fn read(&mut self, output: &mut [u8]) { - for byte in output { + pub(crate) fn read(&mut self, mut output: &mut [u8]) { + while !output.is_empty() { if self.position == self.rate { - self.permute(); + permute(&mut self.state); self.position = 0 } - *byte = (self.state[self.position / 8] >> (8 * (self.position % 8))) as u8; - self.position += 1 + let at = self.position; + let step = if at.is_multiple_of(8) && output.len() >= 8 { 8 } else { 1 }; + let lane = self.state[at / 8].to_le_bytes(); + let (head, rest) = core::mem::take(&mut output).split_at_mut(step); + head.copy_from_slice(&lane[at % 8..at % 8 + step]); + output = rest; + self.position += step + } + } +} +/// Keccak-f[1600] (FIPS 202 §3.3), lane x + 5y at index x + 5y. Each loop but +/// the rounds' has constant bounds, for the compiler to unroll. +fn permute(a: &mut [u64; 25]) { + for &rc in &RC { + let c: [u64; 5] = core::array::from_fn(|x| a[x] ^ a[x + 5] ^ a[x + 10] ^ a[x + 15] ^ a[x + 20]); + for x in 0..5 { + let d = c[(x + 4) % 5] ^ c[(x + 1) % 5].rotate_left(1); + for y in 0..5 { + a[x + 5 * y] ^= d + } + } + // ρ and π: lane (x, y), rotated, moves to (y, 2x + 3y). + let mut b = [0; 25]; + for y in 0..5 { + for x in 0..5 { + b[y + 5 * ((2 * x + 3 * y) % 5)] = a[x + 5 * y].rotate_left(RHO[x + 5 * y]) + } } + for y in 0..5 { + for x in 0..5 { + a[x + 5 * y] = b[x + 5 * y] ^ (!b[(x + 1) % 5 + 5 * y] & b[(x + 2) % 5 + 5 * y]) + } + } + a[0] ^= rc } } pub(crate) fn h(parts: &[&[u8]]) -> [u8; 32] { @@ -135,4 +155,26 @@ mod tests { ] ); } + + /// The lane-wise and byte-wise paths agree: any split of the input into + /// parts, and of the output into reads, gives the one-piece answer. + #[test] + fn splits_do_not_change_the_output() { + let data: [u8; 400] = core::array::from_fn(|i| (i * 7 + 3) as u8); + for rate in [72, 136, 168] { + let mut want = [0; 400]; + Sponge::new(rate, 31, &[&data]).read(&mut want); + for cut in [1, 3, 8, 9, 15, 64, 71, 135, 167, 169, 399] { + let (a, b) = data.split_at(cut); + let mut got = [0; 400]; + let mut sponge = Sponge::new(rate, 31, &[a, b]); + let (head, tail) = got.split_at_mut(cut); + for piece in head.chunks_mut(5) { + sponge.read(piece) + } + sponge.read(tail); + assert_eq!(got, want, "rate {rate}, cut {cut}"); + } + } + } } diff --git a/lib/crypto/src/mlkem.rs b/lib/crypto/src/mlkem.rs index f6cdb91a65..f82611a9f6 100644 --- a/lib/crypto/src/mlkem.rs +++ b/lib/crypto/src/mlkem.rs @@ -1,19 +1,32 @@ //! FIPS 203 ML-KEM-768, for TLS and SSH hybrid exchange. Secret polynomial -//! operations, CBD sampling, hashes and re-encryption have fixed schedules. -//! Division is by compile-time constants; deployed code must be checked to -//! ensure fixed-latency lowering. Matrix rejection sampling depends only on the +//! operations, CBD sampling, hashes and re-encryption have fixed schedules, +//! and no secret value is divided. Matrix rejection sampling depends only on the //! public seed. Public-key decoding rejects coefficients >=3329. Decapsulation //! implements masked implicit rejection, with no validity bit or error oracle. //! Secret keys can only originate from the caller's entropy source, so external //! expanded decapsulation-key validation is outside this API. +//! +//! Coefficients are signed 16-bit representatives, reduced only where a bound +//! requires it (Seiler, IACR ePrint 2018/039). The reductions, transforms and +//! products state their input bounds. use crate::{ Error, Random, ct, keccak::{self, Sponge}, secret::{Secret, wipe}, }; -const Q: u32 = 3329; -struct Poly([u16; 256]); +const Q: i16 = 3329; +/// q⁻¹ mod 2^16. +const Q_INVERSE: i16 = -3327; +const _: () = assert!(Q.wrapping_mul(Q_INVERSE) == 1); +/// 2^16 mod q, the Montgomery factor. +const R: i32 = (1 << 16) % Q as i32; +/// R² mod q: a Montgomery product by it removes one factor 2^-16. +const R_SQUARED: i16 = (R * R % Q as i32) as i16; +/// R²/128 mod q (3303 is 128⁻¹), Algorithm 10's final scale with that removal. +const INVERSE_SCALE: i16 = (R * R % Q as i32 * 3303 % Q as i32) as i16; +const _: () = assert!(128 * 3303 % Q as i32 == 1); +struct Poly([i16; 256]); impl Drop for Poly { fn drop(&mut self) { for word in &mut self.0 { @@ -23,151 +36,217 @@ impl Drop for Poly { core::sync::atomic::compiler_fence(core::sync::atomic::Ordering::SeqCst) } } -fn add(a: u16, b: u16) -> u16 { - ((u32::from(a) + u32::from(b)) % Q) as u16 +/// ⌊n / q⌋ for n < 2^33 / 623, which holds every compressed value, by a +/// multiply and shift: a divide's latency can depend on its operands +/// (KyberSlash, Bernstein et al., IACR ePrint 2024/1049). With +/// M = ⌈2^33 / q⌉, n·M / 2^33 = n / q + 623·n / (q·2^33), and the excess is +/// below 1/q, too small to carry the quotient past the next integer. +fn divide(n: u32) -> u32 { + const M: u64 = (1u64 << 33).div_ceil(Q as u64); + ((u64::from(n) * M) >> 33) as u32 } -fn sub(a: u16, b: u16) -> u16 { - ((u32::from(a) + Q - u32::from(b)) % Q) as u16 +/// a·2^-16 mod q, in (-q, q), for |a| < q·2^15 (Montgomery reduction, in +/// Seiler's signed form): a − t·q is a multiple of 2^16 for t = a·q⁻¹ mod 2^16. +fn montgomery(a: i32) -> i16 { + let t = (a as i16).wrapping_mul(Q_INVERSE); + ((a - i32::from(t) * i32::from(Q)) >> 16) as i16 } -fn mul(a: u16, b: u16) -> u16 { - ((u32::from(a) * u32::from(b)) % Q) as u16 +/// a·b·2^-16 mod q, in (-q, q), for |a·b| < q·2^15. +fn multiply(a: i16, b: i16) -> i16 { + montgomery(i32::from(a) * i32::from(b)) } -fn power(mut base: u16, mut exponent: u32) -> u16 { - let mut result = 1; - while exponent != 0 { - if exponent & 1 == 1 { - result = mul(result, base) +/// a mod q, in [-(q-1)/2, (q-1)/2], for any a: a less q times the nearest +/// integer to a·⌊2^26/q⌉/2^26 (Barrett reduction). +fn barrett(a: i16) -> i16 { + const V: i32 = ((1 << 26) + Q as i32 / 2) / Q as i32; + let t = (V * i32::from(a) + (1 << 25)) >> 26; + (i32::from(a) - t * i32::from(Q)) as i16 +} +/// a mod q, in [0, q). +fn canonical(a: i16) -> i16 { + let r = barrett(a); + r + ((r >> 15) & Q) +} +/// ζ^BitRev7(i) and, for `odd`, γ_i = ζ^(2·BitRev7(i)+1), with ζ = 17: FIPS 203 +/// Appendix A's two tables, each entry times R, which a Montgomery product +/// by it cancels. +const fn powers(odd: bool) -> [i16; 128] { + let mut out = [0; 128]; + let mut i = 0; + while i < 128 { + let reversed = ((i as u8).reverse_bits() >> 1) as u32; + let exponent = if odd { 2 * reversed + 1 } else { reversed }; + let mut power = R; + let mut k = 0; + while k < exponent { + power = power * 17 % Q as i32; + k += 1 } - base = mul(base, base); - exponent >>= 1 + out[i] = power as i16; + i += 1 } - result -} -fn zeta(i: usize) -> u16 { - power(17, (i as u8).reverse_bits() as u32 >> 1) + out } +const ZETAS: [i16; 128] = powers(false); +const GAMMAS: [i16; 128] = powers(true); impl Poly { fn zero() -> Self { Self([0; 256]) } + /// Algorithm 9, for coefficients below q in magnitude. A layer adds less + /// than q to a magnitude, so seven stay below 8q < 2^15. The result is + /// Barrett-reduced. fn ntt(mut self) -> Self { - let mut k = 1; - let mut len = 128; - while len >= 2 { - for start in (0..256).step_by(2 * len) { - let z = zeta(k); - k += 1; - for i in start..start + len { - let t = mul(z, self.0[i + len]); - self.0[i + len] = sub(self.0[i], t); - self.0[i] = add(self.0[i], t) + let mut zetas = ZETAS[1..].iter(); + for layer in (1..8).rev() { + let len = 1 << layer; + for (block, &z) in self.0.chunks_exact_mut(2 * len).zip(&mut zetas) { + let (low, high) = block.split_at_mut(len); + for (a, b) in low.iter_mut().zip(high) { + let t = multiply(z, *b); + *b = *a - t; + *a += t } } - len /= 2 + } + for c in &mut self.0 { + *c = barrett(*c) } self } + /// Algorithm 10 of [`dot`]'s output, whose factor 2^-16 the final scale + /// removes. Coefficients enter and stay below q in magnitude. fn inverse(mut self) -> Self { - let mut k = 127; - let mut len = 2; - while len <= 128 { - for start in (0..256).step_by(2 * len) { - let z = zeta(k); - k -= 1; - for i in start..start + len { - let t = self.0[i]; - self.0[i] = add(t, self.0[i + len]); - self.0[i + len] = mul(z, sub(self.0[i + len], t)) + let mut zetas = ZETAS[1..].iter().rev(); + for layer in 1..8 { + let len = 1 << layer; + for (block, &z) in self.0.chunks_exact_mut(2 * len).zip(&mut zetas) { + let (low, high) = block.split_at_mut(len); + for (a, b) in low.iter_mut().zip(high) { + let t = *a; + *a = barrett(t + *b); + *b = multiply(z, *b - t) } } - len *= 2 } - for word in &mut self.0 { - *word = mul(*word, 3303) + self.scale(INVERSE_SCALE) + } + /// Each coefficient times factor·2^-16, in (-q, q). + fn scale(mut self, factor: i16) -> Self { + for c in &mut self.0 { + *c = multiply(*c, factor) } self } - fn product(&self, q: &Self) -> Self { - let mut out = Self::zero(); - for i in 0..128 { - let gamma = power(17, 2 * ((i as u8).reverse_bits() as u32 >> 1) + 1); - out.0[2 * i] = add( - mul(self.0[2 * i], q.0[2 * i]), - mul(gamma, mul(self.0[2 * i + 1], q.0[2 * i + 1])), - ); - out.0[2 * i + 1] = add( - mul(self.0[2 * i], q.0[2 * i + 1]), - mul(self.0[2 * i + 1], q.0[2 * i]), - ); + fn plus(mut self, b: &Self) -> Self { + for (c, b) in self.0.iter_mut().zip(&b.0) { + *c += b } - out + self } - fn plus(&self, b: &Self) -> Self { - Self(core::array::from_fn(|i| add(self.0[i], b.0[i]))) + fn minus(mut self, b: &Self) -> Self { + for (c, b) in self.0.iter_mut().zip(&b.0) { + *c -= b + } + self } - fn minus(&self, b: &Self) -> Self { - Self(core::array::from_fn(|i| sub(self.0[i], b.0[i]))) + fn canonical(mut self) -> Self { + for c in &mut self.0 { + *c = canonical(*c) + } + self } - fn encode(&self, d: usize, output: &mut [u8]) { - assert_eq!(output.len(), 32 * d); - output.fill(0); - for i in 0..256 { - for bit in 0..d { - let at = i * d + bit; - output[at / 8] |= (((self.0[i] >> bit) & 1) as u8) << (at % 8) + /// ByteEncode_D (FIPS 203 Algorithm 5) of coefficients in [0, 2^D): each + /// eight fill D bytes. + fn encode(&self, output: &mut [u8]) { + assert_eq!(output.len(), 32 * D); + for (coefficients, bytes) in self.0.as_chunks::<8>().0.iter().zip(output.as_chunks_mut::().0) { + let mut word = 0u128; + for (i, &c) in coefficients.iter().enumerate() { + word |= u128::from(c as u16) << (D * i) } + bytes.copy_from_slice(&word.to_le_bytes()[..D]) } } - fn decode(d: usize, input: &[u8]) -> Self { - assert_eq!(input.len(), 32 * d); + /// ByteDecode_D (Algorithm 6), whose coefficients are in [0, 2^D). + fn decode(input: &[u8]) -> Self { + assert_eq!(input.len(), 32 * D); let mut out = Self::zero(); - for i in 0..256 { - for bit in 0..d { - let at = i * d + bit; - out.0[i] |= u16::from((input[at / 8] >> (at % 8)) & 1) << bit + for (coefficients, bytes) in out.0.as_chunks_mut::<8>().0.iter_mut().zip(input.as_chunks::().0) { + let mut word = [0; 16]; + word[..D].copy_from_slice(bytes); + let word = u128::from_le_bytes(word); + for (i, c) in coefficients.iter_mut().enumerate() { + *c = (word >> (D * i)) as i16 & ((1 << D) - 1) } } out } + /// Compress_d of any representatives, into [0, 2^d). fn compress(mut self, d: usize) -> Self { - for word in &mut self.0 { - *word = ((((u32::from(*word) << d) + Q / 2) / Q) & ((1 << d) - 1)) as u16 + for c in &mut self.0 { + let x = canonical(*c) as u32; + *c = (divide((x << d) + Q as u32 / 2) & ((1 << d) - 1)) as i16 } self } + /// Decompress_d of [0, 2^d), into [0, q). fn decompress(mut self, d: usize) -> Self { - for word in &mut self.0 { - *word = ((Q * u32::from(*word) + (1 << (d - 1))) >> d) as u16 + for c in &mut self.0 { + *c = ((Q as u32 * *c as u32 + (1 << (d - 1))) >> d) as i16 } self } } +/// Σ a_j·b_j in the NTT domain (Algorithms 11 and 12), times 2^-16, in +/// (-q, q). Each a_j's coefficients are below q in magnitude and each b_j's +/// at most (q-1)/2, so a pair's three products sum below 4.5q² < q·2^15 and +/// are reduced once. +fn dot(a: &[Poly; 3], b: &[Poly; 3]) -> Poly { + let mut out = Poly::zero(); + for (i, &gamma) in GAMMAS.iter().enumerate() { + let (mut even, mut odd) = (0, 0); + for (a, b) in a.iter().zip(b) { + let (a0, a1) = (i32::from(a.0[2 * i]), i32::from(a.0[2 * i + 1])); + let (b0, b1) = (i32::from(b.0[2 * i]), i32::from(b.0[2 * i + 1])); + even += a0 * b0 + i32::from(montgomery(a1 * b1)) * i32::from(gamma); + odd += a0 * b1 + a1 * b0 + } + out.0[2 * i] = montgomery(even); + out.0[2 * i + 1] = montgomery(odd) + } + out +} +/// SamplePolyCBD_2 (Algorithm 8), in [-2, 2]. fn noise(seed: &[u8; 32], nonce: u8) -> Poly { let mut bytes = Secret([0; 128]); Sponge::new(136, 31, &[seed, &[nonce]]).read(&mut bytes.0); let mut out = Poly::zero(); for i in 0..256 { let bits = (bytes.0[i / 2] >> (4 * (i % 2))) & 15; - let a = u32::from(bits & 1) + u32::from((bits >> 1) & 1); - let b = u32::from((bits >> 2) & 1) + u32::from((bits >> 3) & 1); - out.0[i] = ((Q + a - b) % Q) as u16 + let a = (bits & 1) + ((bits >> 1) & 1); + let b = ((bits >> 2) & 1) + ((bits >> 3) & 1); + out.0[i] = i16::from(a) - i16::from(b) } out } +/// SampleNTT (Algorithm 7), in [0, q). fn sample(rho: &[u8; 32], row: usize, col: usize) -> Poly { let mut xof = Sponge::new(168, 31, &[rho, &[col as u8, row as u8]]); let mut out = Poly::zero(); let mut pos = 0; + let mut block = [0; 168]; while pos < 256 { - let mut b = [0; 3]; - xof.read(&mut b); - for v in [ - u16::from(b[0]) | ((u16::from(b[1]) & 15) << 8), - (u16::from(b[1]) >> 4) | (u16::from(b[2]) << 4), - ] { - if u32::from(v) < Q && pos < 256 { - out.0[pos] = v; - pos += 1 + xof.read(&mut block); + for b in block.as_chunks::<3>().0 { + for v in [ + i16::from(b[0]) | ((i16::from(b[1]) & 15) << 8), + (i16::from(b[1]) >> 4) | (i16::from(b[2]) << 4), + ] { + if v < Q && pos < 256 { + out.0[pos] = v; + pos += 1 + } } } } @@ -177,6 +256,7 @@ pub struct EncapsulationKey { bytes: [u8; 1184], } pub struct DecapsulationKey { + /// ŝ, Barrett-reduced. secret: [Poly; 3], public: EncapsulationKey, hash: [u8; 32], @@ -194,8 +274,8 @@ impl EncapsulationKey { return Err(Error::Length); } for chunk in input[..1152].as_chunks::<384>().0 { - let p = Poly::decode(12, chunk); - if p.0.iter().any(|&x| u32::from(x) >= Q) { + let p = Poly::decode::<12>(chunk); + if p.0.iter().any(|&x| x >= Q) { return Err(Error::InvalidKey); } } @@ -225,24 +305,20 @@ impl EncapsulationKey { let y: [Poly; 3] = core::array::from_fn(|i| noise(coins, i as u8).ntt()); let mut out = [0; 1088]; for i in 0..3 { - let mut sum = Poly::zero(); - for (j, yp) in y.iter().enumerate() { - sum = sum.plus(&sample(rho, j, i).product(yp)) - } - sum.inverse() + let column: [Poly; 3] = core::array::from_fn(|j| sample(rho, j, i)); + dot(&column, &y) + .inverse() .plus(&noise(coins, (i + 3) as u8)) .compress(10) - .encode(10, &mut out[i * 320..i * 320 + 320]); + .encode::<10>(&mut out[i * 320..i * 320 + 320]); } - let mut sum = Poly::zero(); - for (j, yp) in y.iter().enumerate() { - sum = sum.plus(&Poly::decode(12, &self.bytes[j * 384..j * 384 + 384]).product(yp)) - } - sum.inverse() + let t: [Poly; 3] = core::array::from_fn(|j| Poly::decode::<12>(&self.bytes[j * 384..j * 384 + 384])); + dot(&t, &y) + .inverse() .plus(&noise(coins, 6)) - .plus(&Poly::decode(1, message).decompress(1)) + .plus(&Poly::decode::<1>(message).decompress(1)) .compress(4) - .encode(4, &mut out[960..]); + .encode::<4>(&mut out[960..]); out } } @@ -256,11 +332,12 @@ impl DecapsulationKey { let secret: [Poly; 3] = core::array::from_fn(|i| noise(sigma, i as u8).ntt()); let mut bytes = [0; 1184]; for i in 0..3 { - let mut sum = noise(sigma, (i + 3) as u8).ntt(); - for (j, s) in secret.iter().enumerate() { - sum = sum.plus(&sample(rho, i, j).product(s)) - } - sum.encode(12, &mut bytes[i * 384..i * 384 + 384]) + let row: [Poly; 3] = core::array::from_fn(|j| sample(rho, i, j)); + dot(&row, &secret) + .scale(R_SQUARED) + .plus(&noise(sigma, (i + 3) as u8).ntt()) + .canonical() + .encode::<12>(&mut bytes[i * 384..i * 384 + 384]) } bytes[1152..].copy_from_slice(rho); let hash = keccak::h(&[&bytes]); @@ -280,22 +357,14 @@ impl DecapsulationKey { if ciphertext.len() != 1088 { return Err(Error::Length); } - let mut sum = Poly::zero(); - for (i, s) in self.secret.iter().enumerate() { - sum = sum.plus( - &s.product( - &Poly::decode(10, &ciphertext[i * 320..i * 320 + 320]) - .decompress(10) - .ntt(), - ), - ) - } + let u: [Poly; 3] = + core::array::from_fn(|i| Poly::decode::<10>(&ciphertext[i * 320..i * 320 + 320]).decompress(10).ntt()); let mut message = Secret([0; 32]); - Poly::decode(4, &ciphertext[960..]) + Poly::decode::<4>(&ciphertext[960..]) .decompress(4) - .minus(&sum.inverse()) + .minus(&dot(&self.secret, &u).inverse()) .compress(1) - .encode(1, &mut message.0); + .encode::<1>(&mut message.0); let kr = Secret(keccak::g(&[&message.0, &self.hash])); let rejected = Secret(keccak::j(&[&self.rejection.0, ciphertext])); let mut expected = self @@ -333,16 +402,82 @@ mod tests { .map(|i| u8::from_str_radix(&s[2 * i..2 * i + 2], 16).unwrap()) .collect() } - /// Every pair of reduced coefficients, against 64-bit integers. #[test] - #[cfg_attr(miri, ignore = "3329^2 pairs")] - fn reduction_modulo_q_is_exact_over_its_whole_input_range() { - for a in 0..Q as u16 { - for b in 0..Q as u16 { - let (x, y, q) = (u64::from(a), u64::from(b), u64::from(Q)); - assert_eq!(u64::from(add(a, b)), (x + y) % q, "{a} + {b}"); - assert_eq!(u64::from(sub(a, b)), (x + q - y) % q, "{a} - {b}"); - assert_eq!(u64::from(mul(a, b)), x * y % q, "{a} * {b}"); + #[cfg_attr(miri, ignore = "13.8 M quotients")] + fn division_by_q_is_exact_over_its_whole_input_range() { + for n in 0..=((1u64 << 33) / 623) as u32 { + assert_eq!(divide(n), n / Q as u32, "{n}"); + } + } + + /// Each reduction over its whole input range, against 64-bit integers. + #[test] + #[cfg_attr(miri, ignore = "2^28 reductions")] + fn reductions_are_exact_over_their_whole_input_ranges() { + let q = i64::from(Q); + for a in i16::MIN..=i16::MAX { + let (r, c) = (barrett(a), canonical(a)); + assert!((i64::from(r) - i64::from(a)) % q == 0 && 2 * r.abs() < Q, "barrett({a}) = {r}"); + assert_eq!(i64::from(c), i64::from(a).rem_euclid(q), "canonical({a})"); + } + let bound = i32::from(Q) << 15; + for a in -bound + 1..bound { + let r = montgomery(a); + assert!((i64::from(r) * 65536 - i64::from(a)) % q == 0 && r.abs() < Q, "montgomery({a}) = {r}"); + } + } + + /// x mod (X² − γ_i) for each i, the NTT's definition (FIPS 203 §4.3), with + /// γ_i = 17^(2·BitRev7(i)+1), in plain 64-bit arithmetic. + fn ntt_by_definition(x: &[i16; 256]) -> [i64; 256] { + let q = i64::from(Q); + let mut out = [0; 256]; + for i in 0..128 { + let gamma = (0..2 * ((i as u8).reverse_bits() >> 1) as u32 + 1).fold(1, |p, _| p * 17 % q); + let mut power = 1; + for m in 0..128 { + out[2 * i] = (out[2 * i] + i64::from(x[2 * m]) * power).rem_euclid(q); + out[2 * i + 1] = (out[2 * i + 1] + i64::from(x[2 * m + 1]) * power).rem_euclid(q); + power = power * gamma % q + } + } + out + } + + /// The lazily reduced transforms at the edges of their stated input + /// bounds, where an overflow would be, against the definition. + #[test] + fn the_ntt_domain_holds_its_bounds_and_the_definition() { + let q = i64::from(Q); + let r = (1i64 << 16) % q; + let patterns: [fn(usize, i16) -> i16; 4] = + [|_, m| m, |_, m| -m, |i, m| if i % 2 == 0 { m } else { -m }, |i, m| (i as i16 * 97 % (2 * m + 1)) - m]; + for pattern in patterns { + let wide = Poly(core::array::from_fn(|i| pattern(i, Q - 1))); + let narrow = Poly(core::array::from_fn(|i| pattern(i, (Q - 1) / 2))); + + let want = ntt_by_definition(&wide.0); + let ntt = Poly(wide.0).ntt(); + for (i, &c) in ntt.0.iter().enumerate() { + assert!(2 * c.abs() < Q && i64::from(canonical(c)) == want[i], "ntt[{i}] = {c}"); + } + + let inverse = Poly(wide.0).inverse(); + let back = ntt_by_definition(&inverse.0); + for (i, &c) in inverse.0.iter().enumerate() { + assert!(c.abs() < Q && back[i] == (i64::from(wide.0[i]) * r).rem_euclid(q), "inverse[{i}] = {c}"); + } + + let (a, b) = ([Poly(wide.0), Poly(wide.0), Poly(wide.0)], [Poly(narrow.0), Poly(narrow.0), Poly(narrow.0)]); + let product = dot(&a, &b); + for i in 0..128 { + let gamma = (0..2 * ((i as u8).reverse_bits() >> 1) as u32 + 1).fold(1, |p, _| p * 17 % q); + let (a0, a1) = (i64::from(wide.0[2 * i]), i64::from(wide.0[2 * i + 1])); + let (b0, b1) = (i64::from(narrow.0[2 * i]), i64::from(narrow.0[2 * i + 1])); + for (k, want) in [3 * (a0 * b0 + a1 * b1 % q * gamma), 3 * (a0 * b1 + a1 * b0)].into_iter().enumerate() { + let c = product.0[2 * i + k]; + assert!(c.abs() < Q && (i64::from(c) * 65536 - want) % q == 0, "dot[{}] = {c}", 2 * i + k); + } } } } @@ -354,25 +489,26 @@ mod tests { fn compression_rounds_as_fips_203_for_every_input() { for d in [1, 4, 10] { let mut x = Poly::zero(); - for start in (0..Q as u16).step_by(256) { + for start in (0..Q).step_by(256) { for (i, word) in x.0.iter_mut().enumerate() { - *word = (start + i as u16).min(Q as u16 - 1); + *word = (start + i as i16).min(Q - 1); } let compressed = Poly(x.0).compress(d); for (i, &value) in x.0.iter().enumerate() { - let want = ((u64::from(value) << (d + 1)) + u64::from(Q)) / (2 * u64::from(Q)) % (1 << d); - assert_eq!(u64::from(compressed.0[i]), want, "Compress_{d}({value})"); + let want = ((value as u64) << (d + 1)) + Q as u64; + let want = want / (2 * Q as u64) % (1 << d); + assert_eq!(compressed.0[i] as u64, want, "Compress_{d}({value})"); } } let mut y = Poly::zero(); for start in (0..1u32 << d).step_by(256) { for (i, word) in y.0.iter_mut().enumerate() { - *word = (start + i as u32).min((1 << d) - 1) as u16; + *word = (start + i as u32).min((1 << d) - 1) as i16; } let decompressed = Poly(y.0).decompress(d); for (i, &value) in y.0.iter().enumerate() { - let want = (2 * u64::from(Q) * u64::from(value) + (1 << d)) >> (d + 1); - assert_eq!(u64::from(decompressed.0[i]), want, "Decompress_{d}({value})"); + let want = (2 * Q as u64 * value as u64 + (1 << d)) >> (d + 1); + assert_eq!(decompressed.0[i] as u64, want, "Decompress_{d}({value})"); } } } @@ -398,7 +534,7 @@ mod tests { let expected = hex(test["dk"].as_str().unwrap()); let mut actual = Secret([0; 2400]); for (i, p) in key.secret.iter().enumerate() { - p.encode(12, &mut actual.0[384 * i..384 * i + 384]) + Poly(p.0).canonical().encode::<12>(&mut actual.0[384 * i..384 * i + 384]) } actual.0[1152..2336].copy_from_slice(key.public.as_bytes()); actual.0[2336..2368].copy_from_slice(&key.hash);