Add BoxedUint::gcd (#497)

Support for computing the greatest common divisor of two `BoxedUint`s
using Bernstein-Yang

Author: tarcieri

src/modular.rs

@@ -21,7 +21,7 @@ mod monty_form;
 mod reduction;
 
 mod add;
-mod bernstein_yang;
+pub(crate) mod bernstein_yang;
 mod div_by_2;
 mod mul;
 mod pow;

src/modular/bernstein_yang/boxed.rs

@@ -74,6 +74,24 @@ impl Inverter for BoxedBernsteinYangInverter {
     }
 }
 
+/// Returns the greatest common divisor (GCD) of the two given numbers.
+pub(crate) fn gcd(f: &BoxedUint, g: &BoxedUint) -> BoxedUint {
+    let bits_precision = f.bits_precision();
+    let inverse = inv_mod2_62(f.as_words());
+    let f = BoxedInt62L::from(f);
+    let mut g = BoxedInt62L::from(g);
+    let mut d = BoxedInt62L::zero(f.0.len());
+    let e = BoxedInt62L::one(f.0.len());
+
+    let mut f = divsteps(&mut d, &e, &f, &mut g, inverse);
+
+    if f.is_negative() {
+        f = f.neg();
+    }
+
+    f.to_uint(bits_precision)
+}
+
 /// Algorithm `divsteps2` to compute (δₙ, fₙ, gₙ) = divstepⁿ(δ, f, g) as described in Figure 10.1
 /// of <https://eprint.iacr.org/2019/266.pdf>.
 fn divsteps(
@@ -322,6 +340,13 @@ impl BoxedInt62L {
         Self(vec![0; nlimbs].into())
     }
 
+    /// Get the value zero for the given number of limbs.
+    pub fn one(nlimbs: usize) -> Self {
+        let mut ret = Self::zero(nlimbs);
+        ret.0[0] = 0;
+        ret
+    }
+
     /// Widen self to the given number of limbs.
     pub fn widen(self, nlimbs: usize) -> Self {
         let mut limbs = Vec::from(self.0);

src/uint/boxed.rs

@@ -13,6 +13,7 @@ mod div;
 mod div_limb;
 pub(crate) mod encoding;
 mod from;
+mod gcd;
 mod inv_mod;
 mod mul;
 mod mul_mod;

src/uint/boxed/gcd.rs

@@ -0,0 +1,33 @@
+//! Support for computing greatest common divisor of two `BoxedUint`s.
+
+use super::BoxedUint;
+use crate::{modular::bernstein_yang, Odd};
+
+impl Odd<BoxedUint> {
+    /// Compute the greatest common divisor (GCD) of this number and another.
+    pub fn gcd(&self, rhs: &BoxedUint) -> BoxedUint {
+        bernstein_yang::boxed::gcd(self, rhs)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::BoxedUint;
+
+    #[test]
+    fn gcd_relatively_prime() {
+        // Two semiprimes with no common factors
+        let f = BoxedUint::from(59u32 * 67).to_odd().unwrap();
+        let g = BoxedUint::from(61u32 * 71);
+        let gcd = f.gcd(&g);
+        assert_eq!(gcd, BoxedUint::one());
+    }
+
+    #[test]
+    fn gcd_nonprime() {
+        let f = BoxedUint::from(4391633u32).to_odd().unwrap();
+        let g = BoxedUint::from(2022161u32);
+        let gcd = f.gcd(&g);
+        assert_eq!(gcd, BoxedUint::from(1763u32));
+    }
+}

tests/boxed_uint_proptests.rs

@@ -5,6 +5,7 @@
 use core::cmp::Ordering;
 use crypto_bigint::{BoxedUint, CheckedAdd, Integer, Limb, NonZero};
 use num_bigint::{BigUint, ModInverse};
+use num_integer::Integer as _;
 use num_traits::identities::One;
 use proptest::prelude::*;
 
@@ -141,6 +142,22 @@ proptest! {
         prop_assert_eq!(expected_remainder, to_biguint(&actual_remainder));
     }
 
+    #[test]
+    fn gcd((mut f, g) in uint_pair()) {
+        if f.is_even().into() {
+            // Ensure `f` is always odd (required by Bernstein-Yang)
+            f = f.wrapping_add(&BoxedUint::one());
+        }
+
+        let f = f.to_odd().unwrap();
+        let f_bi = to_biguint(&f);
+        let g_bi = to_biguint(&g);
+
+        let expected = to_uint(f_bi.gcd(&g_bi));
+        let actual = f.gcd(&g);
+        assert_eq!(expected, actual);
+    }
+
     #[test]
     fn mod_inv((mut a, mut b) in uint_pair()) {
         if a.is_zero().into() {

tests/uint_proptests.rs

@@ -287,7 +287,6 @@ proptest! {
 
         let expected = to_uint(f_bi.gcd(&g_bi));
         let actual = f.gcd(&g).unwrap();
-
         assert_eq!(expected, actual);
     }