Add binary GCD

src/hazmat/gcd.rs

@@ -3,6 +3,7 @@ use crypto_bigint::{Integer, Limb, NonZero, Word};
 /// Calculates the greatest common divisor of `n` and `m`.
 /// By definition, `gcd(0, m) == m`.
 /// `n` must be non-zero.
+#[inline]
 pub(crate) fn gcd_vartime<T: Integer>(n: &T, m: Word) -> Word {
     // This is an internal function, and it will never be called with `m = 0`.
     // Allowing `m = 0` would require us to have the return type of `Uint<L>`
@@ -16,7 +17,7 @@ pub(crate) fn gcd_vartime<T: Integer>(n: &T, m: Word) -> Word {
     }
 
     // Normalize input: the resulting (a, b) are both small, a >= b, and b != 0.
-    let (mut a, mut b): (Word, Word) = if n.bits() > Word::BITS {
+    let (a, b): (Word, Word) = if n.bits() > Word::BITS {
         // `m` is non-zero, so we can unwrap.
         let r = n.rem_limb(NonZero::new(Limb::from(m)).expect("divisor should be non-zero here"));
         (m, r.0)
@@ -31,15 +32,49 @@ pub(crate) fn gcd_vartime<T: Integer>(n: &T, m: Word) -> Word {
         }
     };
 
-    // Euclidean algorithm.
-    // Binary GCD algorithm could be used here,
-    // but the performance impact of this code is negligible.
+    // Now do standard binary GCD on the two u64s
+    binary_gcd(a, b)
+}
+
+// Binary GCD lifted verbatim from [1], minus the base checks.
+// The identities mentioned in the comments are the following:
+// 1. `gcd(n, 0) = n`: everything divides 0 and n is the largest number that divides n.
+// 2. `gcd(2n, 2m) = 2*gcd(n, m)`: 2 is a common divisor
+// 3. `gcd(n, 2m) = gcd(n, m)` if m is odd: 2 is then not a common divisor.
+// 4. `gcd(n, m) = gcd(n, m - n)` if n, m odd and n <= m
+//
+// As GCD is commutative `gcd(n, m) = gcd(m, n)` those identities still apply if the operands are swapped.
+//
+// [1]: https://en.wikipedia.org/wiki/Binary_GCD_algorithm
+#[inline(always)]
+fn binary_gcd(mut n: u64, mut m: u64) -> u64 {
+    // Using identities 2 and 3:
+    // gcd(2ⁱn, 2ʲm) = 2ᵏ gcd(n, m) with n, m odd and k = min(i, j)
+    // 2ᵏ is the greatest power of two that divides both 2ⁱn and 2ʲm
+    let i = n.trailing_zeros();
+    n >>= i;
+    let j = m.trailing_zeros();
+    m >>= j;
+    let k = core::cmp::min(i, j);
+
     loop {
-        let r = a % b;
-        if r == 0 {
-            return b;
+        // Swap if necessary so n ≤ m
+        if n > m {
+            core::mem::swap(&mut n, &mut m);
         }
-        (a, b) = (b, r)
+
+        // Identity 4: gcd(n, m) = gcd(n, m-n) as n ≤ m and n, m are both odd
+        m -= n;
+        // m is now even
+
+        if m == 0 {
+            // Identity 1: gcd(n, 0) = n
+            // The shift by k is necessary to add back the 2ᵏ factor that was removed before the loop
+            return n << k;
+        }
+
+        // Identity 3: gcd(n, 2ʲ m) = gcd(n, m) as n is odd
+        m >>= m.trailing_zeros();
     }
 }