diff --git a/ml-dsa/src/algebra.rs b/ml-dsa/src/algebra.rs index 559b68a..0742c53 100644 --- a/ml-dsa/src/algebra.rs +++ b/ml-dsa/src/algebra.rs @@ -54,8 +54,55 @@ pub(crate) trait Decompose { fn decompose(self) -> (Elem, Elem); } +/// Constant-time division by a compile-time constant divisor. +/// +/// This trait provides a constant-time alternative to the hardware division +/// instruction, which has variable timing based on operand values. +/// Uses Barrett reduction to compute `x / M` where M is a compile-time constant. +pub(crate) trait ConstantTimeDiv: Unsigned { + /// Bit shift for Barrett reduction, chosen to provide sufficient precision + const CT_DIV_SHIFT: usize; + /// Precomputed multiplier: ceil(2^SHIFT / M) + const CT_DIV_MULTIPLIER: u64; + + /// Perform constant-time division of x by `Self::U32` + /// Requires: x < Q (the field modulus, ~2^23) + #[allow(clippy::inline_always)] // Required for constant-time guarantees in crypto code + #[inline(always)] + fn ct_div(x: u32) -> u32 { + // Barrett reduction: q = (x * MULTIPLIER) >> SHIFT + // This gives us floor(x / M) for x < 2^SHIFT / MULTIPLIER * M + let x64 = u64::from(x); + let quotient = (x64 * Self::CT_DIV_MULTIPLIER) >> Self::CT_DIV_SHIFT; + // SAFETY: quotient is guaranteed to fit in u32 because: + // - x < Q (~2^23), so quotient = x / M < x < 2^23 < 2^32 + #[allow(clippy::cast_possible_truncation, clippy::as_conversions)] + let result = quotient as u32; + result + } +} + +impl ConstantTimeDiv for M +where + M: Unsigned, +{ + // Use a shift that provides enough precision for the ML-DSA field (Q ~ 2^23) + // We need SHIFT > log2(Q) + log2(M) to ensure accuracy + // With Q < 2^24 and M < 2^20, SHIFT = 48 is sufficient + const CT_DIV_SHIFT: usize = 48; + + // Precompute the multiplier at compile time + // We add (M-1) before dividing to get ceiling division, ensuring we never underestimate + #[allow(clippy::integer_division_remainder_used)] + const CT_DIV_MULTIPLIER: u64 = (1u64 << Self::CT_DIV_SHIFT).div_ceil(M::U64); +} + impl Decompose for Elem { // Algorithm 36 Decompose + // + // This implementation uses constant-time division to avoid timing side-channels. + // The original algorithm used hardware division which has variable timing based + // on operand values, potentially leaking secret information during signing. fn decompose(self) -> (Elem, Elem) { let r_plus = self.clone(); let r0 = r_plus.mod_plus_minus::(); @@ -63,8 +110,9 @@ impl Decompose for Elem { if r_plus - r0 == Elem::new(BaseField::Q - 1) { (Elem::new(0), r0 - Elem::new(1)) } else { - let mut r1 = r_plus - r0; - r1.0 /= TwoGamma2::U32; + let diff = r_plus - r0; + // Use constant-time division instead of hardware division + let r1 = Elem::new(TwoGamma2::ct_div(diff.0)); (r1, r0) } } diff --git a/ml-dsa/src/ntt.rs b/ml-dsa/src/ntt.rs index a0b0714..da63bb1 100644 --- a/ml-dsa/src/ntt.rs +++ b/ml-dsa/src/ntt.rs @@ -50,28 +50,46 @@ pub(crate) trait Ntt { fn ntt(&self) -> Self::Output; } +/// Constant-time NTT butterfly layer. +/// +/// Uses const generics to ensure loop bounds are compile-time constants, +/// avoiding UDIV instructions from runtime `step_by` calculations. +#[allow(clippy::inline_always)] // Required for constant-time guarantees in crypto code +#[inline(always)] +fn ntt_layer(w: &mut [Elem; 256], m: &mut usize) { + for i in 0..ITERATIONS { + let start = i * 2 * LEN; + *m += 1; + let z = ZETA_POW_BITREV[*m]; + for j in start..(start + LEN) { + let t = z * w[j + LEN]; + w[j + LEN] = w[j] - t; + w[j] = w[j] + t; + } + } +} + impl Ntt for Polynomial { type Output = NttPolynomial; // Algorithm 41 NTT + // + // This implementation uses const-generic helper functions to ensure all loop + // bounds are compile-time constants, avoiding potential UDIV instructions. fn ntt(&self) -> Self::Output { - let mut w = self.0.clone(); - + let mut w: [Elem; 256] = self.0.clone().into(); let mut m = 0; - for len in [128, 64, 32, 16, 8, 4, 2, 1] { - for start in (0..256).step_by(2 * len) { - m += 1; - let z = ZETA_POW_BITREV[m]; - for j in start..(start + len) { - let t = z * w[j + len]; - w[j + len] = w[j] - t; - w[j] = w[j] + t; - } - } - } + ntt_layer::<128, 1>(&mut w, &mut m); + ntt_layer::<64, 2>(&mut w, &mut m); + ntt_layer::<32, 4>(&mut w, &mut m); + ntt_layer::<16, 8>(&mut w, &mut m); + ntt_layer::<8, 16>(&mut w, &mut m); + ntt_layer::<4, 32>(&mut w, &mut m); + ntt_layer::<2, 64>(&mut w, &mut m); + ntt_layer::<1, 128>(&mut w, &mut m); - NttPolynomial::new(w) + NttPolynomial::new(w.into()) } } @@ -89,30 +107,51 @@ pub(crate) trait NttInverse { fn ntt_inverse(&self) -> Self::Output; } +/// Constant-time inverse NTT butterfly layer. +/// +/// Uses const generics to ensure loop bounds are compile-time constants, +/// avoiding UDIV instructions from runtime `step_by` calculations. +#[allow(clippy::inline_always)] // Required for constant-time guarantees in crypto code +#[inline(always)] +fn ntt_inverse_layer( + w: &mut [Elem; 256], + m: &mut usize, +) { + for i in 0..ITERATIONS { + let start = i * 2 * LEN; + *m -= 1; + let z = -ZETA_POW_BITREV[*m]; + for j in start..(start + LEN) { + let t = w[j]; + w[j] = t + w[j + LEN]; + w[j + LEN] = z * (t - w[j + LEN]); + } + } +} + impl NttInverse for NttPolynomial { type Output = Polynomial; // Algorithm 42 NTT^{−1} + // + // This implementation uses const-generic helper functions to ensure all loop + // bounds are compile-time constants, avoiding potential UDIV instructions. fn ntt_inverse(&self) -> Self::Output { const INVERSE_256: Elem = Elem::new(8_347_681); - let mut w = self.0.clone(); - + let mut w: [Elem; 256] = self.0.clone().into(); let mut m = 256; - for len in [1, 2, 4, 8, 16, 32, 64, 128] { - for start in (0..256).step_by(2 * len) { - m -= 1; - let z = -ZETA_POW_BITREV[m]; - for j in start..(start + len) { - let t = w[j]; - w[j] = t + w[j + len]; - w[j + len] = z * (t - w[j + len]); - } - } - } + ntt_inverse_layer::<1, 128>(&mut w, &mut m); + ntt_inverse_layer::<2, 64>(&mut w, &mut m); + ntt_inverse_layer::<4, 32>(&mut w, &mut m); + ntt_inverse_layer::<8, 16>(&mut w, &mut m); + ntt_inverse_layer::<16, 8>(&mut w, &mut m); + ntt_inverse_layer::<32, 4>(&mut w, &mut m); + ntt_inverse_layer::<64, 2>(&mut w, &mut m); + ntt_inverse_layer::<128, 1>(&mut w, &mut m); - INVERSE_256 * &Polynomial::new(w) + INVERSE_256 * &Polynomial::new(w.into()) } }