From acbc887df670e4941751b3a6296e86e931008811 Mon Sep 17 00:00:00 2001 From: Bogdan Opanchuk Date: Fri, 22 Dec 2023 16:11:24 -0800 Subject: [PATCH] Add vartime GCD and Jacobi symbol functions --- src/gcd.rs | 132 +++++++++++++++++++++++ src/jacobi.rs | 237 ++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 4 + src/traits.rs | 8 ++ src/uint/boxed/div.rs | 9 +- src/uint/div.rs | 8 +- src/uint/gcd.rs | 3 +- 7 files changed, 397 insertions(+), 4 deletions(-) create mode 100644 src/gcd.rs create mode 100644 src/jacobi.rs diff --git a/src/gcd.rs b/src/gcd.rs new file mode 100644 index 000000000..75399ba77 --- /dev/null +++ b/src/gcd.rs @@ -0,0 +1,132 @@ +use crate::{Integer, Limb, NonZero, Word}; + +/// Compute the greatest common divisor (GCD) of this number and another. +/// Variable time in both arguments. +pub fn gcd_vartime(lhs: &T, rhs: &T) -> T { + // This we can check since it doesn't affect the return type, + // even though `n` will not be 0 either in the application. + if lhs.is_zero().into() { + return rhs.clone(); + } + + if rhs.is_zero().into() { + return lhs.clone(); + } + + // Normalize input: the resulting (a, b) are both small, a >= b, and b != 0. + + let mut a = lhs.clone(); + let mut b = rhs.clone(); + + let mut a_ref = &mut a; + let mut b_ref = &mut b; + + if a_ref > b_ref { + (a_ref, b_ref) = (b_ref, a_ref); + } + + // Euclidean algorithm. + loop { + let r = a_ref.rem_vartime(&NonZero::new(b_ref.clone()).unwrap()); + if r.is_zero().into() { + return b_ref.clone(); + } + + (a_ref, b_ref) = (b_ref, a_ref); + *b_ref = r; + } +} + +/// Compute the greatest common divisor (GCD) of this number and another. +/// Variable time in both arguments. +pub fn gcd_limb_vartime(lhs: &T, rhs: NonZero) -> Limb { + // Allowing `rhs = 0` would require us to have the return type of `T` + // (since `gcd(n, 0) = n`). + + // This we can check since it doesn't affect the return type, + // even though `n` will not be 0 either in the application. + if lhs.is_zero().into() { + return rhs.0; + } + + // Normalize input: the resulting (a, b) are both small, `a >= b`, and `b != 0`. + let (mut a, mut b): (Word, Word) = if lhs.bits_vartime() > Limb::BITS { + let n = lhs.rem_vartime(&NonZero::new(T::from(rhs.0)).unwrap()); + (rhs.0 .0, n.as_ref()[0].0) + } else { + // In this branch `lhs` fits in a `Limb`, + // so we can safely take the first limb and use it. + let n = lhs.as_ref()[0]; + if n > rhs.0 { + (n.0, rhs.0 .0) + } else { + (rhs.0 .0, n.0) + } + }; + + // Euclidean algorithm. + // Binary GCD algorithm could be used here, + // but the performance impact of this code is negligible. + loop { + let r = a % b; + if r == 0 { + return Limb(b); + } + (a, b) = (b, r) + } +} + +#[cfg(test)] +mod tests { + use crate::{Limb, NonZero, Word, U128}; + use num_bigint::BigUint; + use num_integer::Integer; + use proptest::prelude::*; + + use super::gcd_limb_vartime; + + #[test] + fn corner_cases() { + assert_eq!( + gcd_limb_vartime(&U128::from(0u64), NonZero::new(Limb(5)).unwrap()), + Limb(5) + ); + assert_eq!( + gcd_limb_vartime(&U128::from(1u64), NonZero::new(Limb(11 * 13 * 19)).unwrap()), + Limb(1) + ); + assert_eq!( + gcd_limb_vartime(&U128::from(7u64 * 11 * 13), NonZero::new(Limb(1)).unwrap()), + Limb(1) + ); + assert_eq!( + gcd_limb_vartime( + &U128::from(7u64 * 11 * 13), + NonZero::new(Limb(11 * 13 * 19)).unwrap() + ), + Limb(11 * 13) + ); + } + + prop_compose! { + fn uint()(bytes in any::<[u8; 16]>()) -> U128 { + U128::from_le_slice(&bytes) | U128::ONE + } + } + + proptest! { + #[test] + fn fuzzy(m in any::(), n in uint()) { + if m == 0 { + return Ok(()); + } + + let m_bi = BigUint::from(m); + let n_bi = BigUint::from_bytes_be(n.to_be_bytes().as_ref()); + let gcd_ref: Word = n_bi.gcd(&m_bi).try_into().unwrap(); + + let gcd_test = gcd_limb_vartime(&n, NonZero::new(Limb(m)).unwrap()).0; + assert_eq!(gcd_test, gcd_ref); + } + } +} diff --git a/src/jacobi.rs b/src/jacobi.rs new file mode 100644 index 000000000..b9ceb1c7b --- /dev/null +++ b/src/jacobi.rs @@ -0,0 +1,237 @@ +//! Jacobi symbol calculation. + +use crate::{Integer, Limb, NonZero, Odd, Word}; + +/// Possible values of Jacobi symbol. +#[derive(Copy, Clone, Debug, PartialEq, Eq)] +pub enum JacobiSymbol { + /// `0` + Zero, + /// `1` + One, + /// `-1` + MinusOne, +} + +impl core::ops::Neg for JacobiSymbol { + type Output = Self; + fn neg(self) -> Self { + match self { + Self::Zero => Self::Zero, + Self::One => Self::MinusOne, + Self::MinusOne => Self::One, + } + } +} + +// A helper trait to generalize some functions over Word and Uint. +trait SmallMod { + fn mod8(&self) -> Word; + fn mod4(&self) -> Word; +} + +impl SmallMod for Word { + fn mod8(&self) -> Word { + self & 7 + } + fn mod4(&self) -> Word { + self & 3 + } +} + +impl SmallMod for T { + fn mod8(&self) -> Word { + self.as_ref()[0].0 & 7 + } + fn mod4(&self) -> Word { + self.as_ref()[0].0 & 3 + } +} + +/// Transforms `(a/p)` -> `(r/p)` for odd `p`, where the resulting `r` is odd, and `a = r * 2^s`. +/// Takes a Jacobi symbol value, and returns `r` and the new Jacobi symbol, +/// negated if the transformation changes parity. +/// +/// Note that the returned `r` is odd. +fn reduce_numerator(j: JacobiSymbol, a: Word, p: &V) -> (JacobiSymbol, Word) { + let p_mod_8 = p.mod8(); + let s = a.trailing_zeros(); + let j = if (s & 1) == 1 && (p_mod_8 == 3 || p_mod_8 == 5) { + -j + } else { + j + }; + (j, a >> s) +} + +/// Transforms `(a/p)` -> `(p/a)` for odd and coprime `a` and `p`. +/// Takes a Jacobi symbol value, and returns the swapped pair and the new Jacobi symbol, +/// negated if the transformation changes parity. +fn swap(j: JacobiSymbol, a: T, p: V) -> (JacobiSymbol, V, T) { + let j = if a.mod4() == 1 || p.mod4() == 1 { + j + } else { + -j + }; + (j, p, a) +} + +/// Returns the Jacobi symbol `(a/p)` given an odd `p`. Panics on even `p`. +pub fn jacobi_symbol_vartime(a: i32, p_long: &Odd) -> JacobiSymbol { + let p_long = p_long.0.clone(); + + let result = JacobiSymbol::One; // Keep track of all the sign flips here. + + // Deal with a negative `a` first: + // (-a/n) = (-1/n) * (a/n) + // = (-1)^((n-1)/2) * (a/n) + // = (-1 if n = 3 mod 4 else 1) * (a/n) + let (result, a_pos) = { + let result = if a < 0 && p_long.mod4() == 3 { + -result + } else { + result + }; + (result, a.abs_diff(0)) + }; + + // A degenerate case. + if a_pos == 1 || p_long == T::one() { + return result; + } + + let a_limb = Limb::from(a_pos); + + // Normalize input: at the end we want `a < p`, `p` odd, and both fitting into a `Word`. + let (result, a, p): (JacobiSymbol, Word, Word) = if p_long.bits_vartime() <= Limb::BITS { + let a = a_limb.0; + let p = p_long.as_ref()[0].0; + (result, a % p, p) + } else { + let (result, a) = reduce_numerator(result, a_limb.0, &p_long); + if a == 1 { + return result; + } + let (result, a_long, p) = swap(result, a, p_long); + // Can unwrap here, since `p` is swapped with `a`, + // and `a` would be odd after `reduce_numerator()`. + let a = a_long.rem_vartime(&NonZero::new(T::from(p)).unwrap()); + (result, a.as_ref()[0].0, p) + }; + + let mut result = result; + let mut a = a; + let mut p = p; + + loop { + if a == 0 { + return JacobiSymbol::Zero; + } + + // At this point `p` is odd (either coming from outside of the `loop`, + // or from the previous iteration, where a previously reduced `a` + // was swapped into its place), so we can call this. + (result, a) = reduce_numerator(result, a, &p); + + if a == 1 { + return result; + } + + // At this point both `a` and `p` are odd: `p` was odd before, + // and `a` is odd after `reduce_numerator()`. + // Note that technically `swap()` only returns a valid `result` if `a` and `p` are coprime. + // But if they are not, we will return `Zero` eventually, + // which is not affected by any sign changes. + (result, a, p) = swap(result, a, p); + + a %= p; + } +} + +#[cfg(test)] +mod tests { + use crate::{Odd, U128}; + use num_bigint::{BigInt, Sign}; + use num_modular::ModularSymbols; + use proptest::prelude::*; + + use super::{jacobi_symbol_vartime, JacobiSymbol}; + + #[test] + fn jacobi_symbol_neg_zero() { + // This does not happen during normal operation, since we return zero as soon as we get it. + // So just covering it for the completness' sake. + assert_eq!(-JacobiSymbol::Zero, JacobiSymbol::Zero); + } + + // Reference from `num-modular` - supports long `p`, but only positive `a`. + fn jacobi_symbol_ref(a: i32, p: &U128) -> JacobiSymbol { + let a_bi = BigInt::from(a); + let p_bi = BigInt::from_bytes_be(Sign::Plus, p.to_be_bytes().as_ref()); + let j = a_bi.jacobi(&p_bi); + if j == 1 { + JacobiSymbol::One + } else if j == -1 { + JacobiSymbol::MinusOne + } else { + JacobiSymbol::Zero + } + } + + #[test] + fn small_values() { + // Test small values, using a reference implementation. + for a in -31i32..31 { + for p in (1u32..31).step_by(2) { + let p_long = U128::from(p); + let j_ref = jacobi_symbol_ref(a, &p_long); + let j = jacobi_symbol_vartime(a, &Odd::new(p_long).unwrap()); + assert_eq!(j, j_ref); + } + } + } + + #[test] + fn big_values() { + // a = x, p = x * y, where x and y are big primes. Should give 0. + let a = 2147483647i32; // 2^31 - 1, a prime + let p = U128::from_be_hex("000000007ffffffeffffffe28000003b"); // (2^31 - 1) * (2^64 - 59) + assert_eq!( + jacobi_symbol_vartime(a, &Odd::new(p).unwrap()), + JacobiSymbol::Zero + ); + assert_eq!(jacobi_symbol_ref(a, &p), JacobiSymbol::Zero); + + // a = x^2 mod p, should give 1. + let a = 659456i32; // Obtained from x = 2^70 + let p = U128::from_be_hex("ffffffffffffffffffffffffffffff5f"); // 2^128 - 161 - not a prime + assert_eq!( + jacobi_symbol_vartime(a, &Odd::new(p).unwrap()), + JacobiSymbol::One + ); + assert_eq!(jacobi_symbol_ref(a, &p), JacobiSymbol::One); + + let a = i32::MIN; // -2^31, check that no overflow occurs + let p = U128::from_be_hex("000000007ffffffeffffffe28000003b"); // (2^31 - 1) * (2^64 - 59) + assert_eq!( + jacobi_symbol_vartime(a, &Odd::new(p).unwrap()), + JacobiSymbol::One + ); + assert_eq!(jacobi_symbol_ref(a, &p), JacobiSymbol::One); + } + + prop_compose! { + fn odd_uint()(bytes in any::<[u8; 16]>()) -> Odd { + Odd::new(U128::from_le_slice(&bytes) | U128::ONE).unwrap() + } + } + + proptest! { + #[test] + fn fuzzy(a in any::(), p in odd_uint()) { + let j_ref = jacobi_symbol_ref(a, &p); + let j = jacobi_symbol_vartime(a, &p); + assert_eq!(j, j_ref); + } + } +} diff --git a/src/lib.rs b/src/lib.rs index b2682bd25..e9ffd456c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -168,6 +168,8 @@ pub mod modular; mod array; mod checked; mod const_choice; +mod gcd; +mod jacobi; mod limb; mod non_zero; mod odd; @@ -179,6 +181,8 @@ mod wrapping; pub use crate::{ checked::Checked, const_choice::{ConstChoice, ConstCtOption}, + gcd::{gcd_limb_vartime, gcd_vartime}, + jacobi::{jacobi_symbol_vartime, JacobiSymbol}, limb::{Limb, WideWord, Word}, non_zero::NonZero, odd::Odd, diff --git a/src/traits.rs b/src/traits.rs index 97a45c077..48ed6e34e 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -115,6 +115,7 @@ pub trait Integer: + From + From + From + + From + Mul + for<'a> Mul<&'a Self, Output = Self> + MulMod @@ -123,6 +124,7 @@ pub trait Integer: + Ord + Rem, Output = Self> + for<'a> Rem<&'a NonZero, Output = Self> + + RemVartime + Send + Sized + Shl @@ -441,6 +443,12 @@ pub trait SquareAssign { fn square_assign(&mut self); } +/// Support for variable time modulo reduction. +pub trait RemVartime: Sized { + /// Computes `self % other`, possibly in variable time. + fn rem_vartime(&self, other: &NonZero) -> Self; +} + /// Constant-time exponentiation. pub trait Pow { /// Raises to the `exponent` power. diff --git a/src/uint/boxed/div.rs b/src/uint/boxed/div.rs index d230fc1d3..050382986 100644 --- a/src/uint/boxed/div.rs +++ b/src/uint/boxed/div.rs @@ -1,7 +1,8 @@ //! [`BoxedUint`] division operations. use crate::{ - uint::boxed, BoxedUint, CheckedDiv, ConstantTimeSelect, Limb, NonZero, Reciprocal, Wrapping, + uint::boxed, BoxedUint, CheckedDiv, ConstantTimeSelect, Limb, NonZero, Reciprocal, RemVartime, + Wrapping, }; use core::ops::{Div, DivAssign, Rem, RemAssign}; use subtle::{Choice, ConstantTimeEq, ConstantTimeLess, CtOption}; @@ -295,6 +296,12 @@ impl RemAssign> for BoxedUint { } } +impl RemVartime for BoxedUint { + fn rem_vartime(&self, other: &NonZero) -> Self { + Self::rem_vartime(self, other) + } +} + #[cfg(test)] mod tests { use super::{BoxedUint, NonZero}; diff --git a/src/uint/div.rs b/src/uint/div.rs index d650235d8..962f2f62d 100644 --- a/src/uint/div.rs +++ b/src/uint/div.rs @@ -1,7 +1,7 @@ //! [`Uint`] division operations. use super::div_limb::{div_rem_limb_with_reciprocal, Reciprocal}; -use crate::{CheckedDiv, ConstChoice, Limb, NonZero, Uint, Word, Wrapping}; +use crate::{CheckedDiv, ConstChoice, Limb, NonZero, RemVartime, Uint, Word, Wrapping}; use core::ops::{Div, DivAssign, Rem, RemAssign}; use subtle::CtOption; @@ -584,6 +584,12 @@ impl RemAssign<&NonZero>> for Wrapping RemVartime for Uint { + fn rem_vartime(&self, other: &NonZero) -> Self { + Self::rem_vartime(self, other) + } +} + #[cfg(test)] mod tests { use crate::{Limb, NonZero, Uint, Word, U256}; diff --git a/src/uint/gcd.rs b/src/uint/gcd.rs index 308d66342..b164ec7a2 100644 --- a/src/uint/gcd.rs +++ b/src/uint/gcd.rs @@ -1,7 +1,6 @@ //! Support for computing greatest common divisor of two `Uint`s. -use super::Uint; -use crate::{modular::BernsteinYangInverter, ConstCtOption, PrecomputeInverter}; +use crate::{modular::BernsteinYangInverter, ConstCtOption, PrecomputeInverter, Uint}; impl Uint where