Add support for computing Uint::gcd with even modulus (#617)

Uses a similar method to `Uint::inv_mod` for extending `Uint::gcd` with
support for an even modulus, first dividing out `2^k` from both operands
where this operation will ensure at least one value is odd, and then
applying Bernstein-Yang to the odd modulus, multiplying the result by
`2^k`.

src/uint/gcd.rs

@@ -1,18 +1,34 @@
-//! Support for computing greatest common divisor of two `Uint`s.
+//! Support for computing the greatest common divisor of two `Uint`s.
 
-use crate::{modular::BernsteinYangInverter, ConstCtOption, Gcd, Odd, PrecomputeInverter, Uint};
+use crate::{
+    modular::BernsteinYangInverter, ConstChoice, ConstCtOption, Gcd, Odd, PrecomputeInverter, Uint,
+};
 use subtle::CtOption;
 
 impl<const SAT_LIMBS: usize, const UNSAT_LIMBS: usize> Uint<SAT_LIMBS>
 where
     Odd<Self>: PrecomputeInverter<Inverter = BernsteinYangInverter<SAT_LIMBS, UNSAT_LIMBS>>,
 {
     /// Compute the greatest common divisor (GCD) of this number and another.
-    ///
-    /// Returns none in the event that `self` is even (i.e. `self` MUST be odd). However, `rhs` may be even.
     pub const fn gcd(&self, rhs: &Self) -> ConstCtOption<Self> {
-        let ret = <Odd<Self> as PrecomputeInverter>::Inverter::gcd(self, rhs);
-        ConstCtOption::new(ret, self.is_odd())
+        let k1 = self.trailing_zeros();
+        let k2 = rhs.trailing_zeros();
+
+        // Select the smaller of the two `k` values, making 2^k the common even divisor
+        let k = ConstChoice::from_u32_lt(k2, k1).select_u32(k1, k2);
+
+        // Decompose `self` and `rhs` into `s{1, 2} * 2^k` where either `s1` or `s2` is odd
+        let s1 = self.overflowing_shr(k).unwrap_or(Self::ZERO);
+        let s2 = rhs.overflowing_shr(k).unwrap_or(Self::ZERO);
+
+        let f = Self::select(&s1, &s2, s2.is_odd().not());
+        let g = Self::select(&s1, &s2, s2.is_odd());
+
+        let ret = <Odd<Self> as PrecomputeInverter>::Inverter::gcd(&f, &g);
+        ConstCtOption::new(
+            ret.overflowing_shl(k).unwrap_or(Self::ZERO),
+            f.is_nonzero().and(g.is_odd()),
+        )
     }
 }
 
@@ -61,8 +77,9 @@ mod tests {
 
     #[test]
     fn gcd_zero() {
-        let f = U256::ZERO;
-        assert!(f.gcd(&U256::ONE).is_none().is_true_vartime());
+        assert!(U256::ZERO.gcd(&U256::ZERO).is_none().is_true_vartime());
+        assert!(U256::ZERO.gcd(&U256::ONE).is_none().is_true_vartime());
+        assert!(U256::ONE.gcd(&U256::ZERO).is_none().is_true_vartime());
     }
 
     #[test]
@@ -71,4 +88,14 @@ mod tests {
         assert_eq!(U256::ONE, f.gcd(&U256::ONE).unwrap());
         assert_eq!(U256::ONE, f.gcd(&U256::from(2u8)).unwrap());
     }
+
+    #[test]
+    fn gcd_two() {
+        let f = U256::from_u8(2);
+        assert_eq!(f, f.gcd(&f).unwrap());
+
+        let g = U256::from_u8(4);
+        assert_eq!(f, f.gcd(&g).unwrap());
+        assert_eq!(f, g.gcd(&f).unwrap());
+    }
 }

tests/uint.rs

@@ -5,7 +5,7 @@ mod common;
 use common::to_biguint;
 use crypto_bigint::{
     modular::{MontyForm, MontyParams},
-    Encoding, Integer, Limb, NonZero, Odd, Word, U256,
+    Encoding, Limb, NonZero, Odd, Word, U256,
 };
 use num_bigint::BigUint;
 use num_integer::Integer as _;
@@ -275,12 +275,7 @@ proptest! {
     }
 
     #[test]
-    fn gcd(mut f in uint(), g in uint()) {
-        if f.is_even().into() {
-            // Ensure `f` is always odd (required by Bernstein-Yang)
-            f = f.wrapping_add(&U256::ONE);
-        }
-
+    fn gcd(f in uint(), g in uint()) {
         let f_bi = to_biguint(&f);
         let g_bi = to_biguint(&g);