Implement Stein's algorithm for gcd

Cargo.toml

@@ -14,6 +14,9 @@ readme = "README.md"
 [[bench]]
 name = "bigint"
 
+[[bench]]
+name = "gcd"
+
 [[bench]]
 harness = false
 name = "shootout-pidigits"

benches/gcd.rs

@@ -0,0 +1,84 @@
+#![feature(test)]
+
+extern crate test;
+extern crate num_bigint;
+extern crate num_integer;
+extern crate num_traits;
+extern crate rand;
+
+use test::Bencher;
+use num_bigint::{BigUint, RandBigInt};
+use num_integer::Integer;
+use num_traits::Zero;
+use rand::{SeedableRng, StdRng};
+
+fn get_rng() -> StdRng {
+    let seed: &[_] = &[1, 2, 3, 4];
+    SeedableRng::from_seed(seed)
+}
+
+fn bench(b: &mut Bencher, bits: usize, gcd: fn(&BigUint, &BigUint) -> BigUint) {
+    let mut rng = get_rng();
+    let x = rng.gen_biguint(bits);
+    let y = rng.gen_biguint(bits);
+
+    assert_eq!(euclid(&x, &y), x.gcd(&y));
+
+    b.iter(|| gcd(&x, &y));
+}
+
+
+fn euclid(x: &BigUint, y: &BigUint) -> BigUint {
+    // Use Euclid's algorithm
+    let mut m = x.clone();
+    let mut n = y.clone();
+    while !m.is_zero() {
+        let temp = m;
+        m = n % &temp;
+        n = temp;
+    }
+    return n;
+}
+
+#[bench]
+fn gcd_euclid_0064(b: &mut Bencher) {
+    bench(b, 64, euclid);
+}
+
+#[bench]
+fn gcd_euclid_0256(b: &mut Bencher) {
+    bench(b, 256, euclid);
+}
+
+#[bench]
+fn gcd_euclid_1024(b: &mut Bencher) {
+    bench(b, 1024, euclid);
+}
+
+#[bench]
+fn gcd_euclid_4096(b: &mut Bencher) {
+    bench(b, 4096, euclid);
+}
+
+
+// Integer for BigUint now uses Stein for gcd
+
+#[bench]
+fn gcd_stein_0064(b: &mut Bencher) {
+    bench(b, 64, BigUint::gcd);
+}
+
+#[bench]
+fn gcd_stein_0256(b: &mut Bencher) {
+    bench(b, 256, BigUint::gcd);
+}
+
+#[bench]
+fn gcd_stein_1024(b: &mut Bencher) {
+    bench(b, 1024, BigUint::gcd);
+}
+
+#[bench]
+fn gcd_stein_4096(b: &mut Bencher) {
+    bench(b, 4096, BigUint::gcd);
+}

src/algorithms.rs

@@ -591,9 +591,12 @@ pub fn biguint_shr(n: Cow<BigUint>, bits: usize) -> BigUint {
     if n_unit >= n.data.len() {
         return Zero::zero();
     }
-    let mut data = match n_unit {
-        0 => n.into_owned().data,
-        _ => n.data[n_unit..].to_vec(),
+    let mut data = match n {
+        Cow::Borrowed(n) => n.data[n_unit..].to_vec(),
+        Cow::Owned(mut n) => {
+            n.data.drain(..n_unit);
+            n.data
+        }
     };
 
     let n_bits = bits % big_digit::BITS;

src/bigint.rs

@@ -4,6 +4,7 @@ use std::str::{self, FromStr};
 use std::fmt;
 use std::cmp::Ordering::{self, Less, Greater, Equal};
 use std::{i64, u64};
+#[allow(unused)]
 use std::ascii::AsciiExt;
 
 #[cfg(feature = "serde")]

src/biguint.rs

@@ -7,9 +7,11 @@ use std::ops::{Add, BitAnd, BitOr, BitXor, Div, Mul, Neg, Rem, Shl, Shr, Sub,
 use std::str::{self, FromStr};
 use std::fmt;
 use std::cmp;
+use std::mem;
 use std::cmp::Ordering::{self, Less, Greater, Equal};
 use std::{f32, f64};
 use std::{u8, u64};
+#[allow(unused)]
 use std::ascii::AsciiExt;
 
 #[cfg(feature = "serde")]
@@ -396,7 +398,8 @@ impl<'a> Shr<usize> for &'a BigUint {
 impl ShrAssign<usize> for BigUint {
     #[inline]
     fn shr_assign(&mut self, rhs: usize) {
-        *self = biguint_shr(Cow::Borrowed(&*self), rhs);
+        let n = mem::replace(self, BigUint::zero());
+        *self = n >> rhs;
     }
 }
 
@@ -940,16 +943,34 @@ impl Integer for BigUint {
     ///
     /// The result is always positive.
     #[inline]
-    fn gcd(&self, other: &BigUint) -> BigUint {
-        // Use Euclid's algorithm
-        let mut m = (*self).clone();
-        let mut n = (*other).clone();
+    fn gcd(&self, other: &Self) -> Self {
+        // Stein's algorithm
+        if self.is_zero() {
+            return other.clone();
+        }
+        if other.is_zero() {
+            return self.clone();
+        }
+        let mut m = self.clone();
+        let mut n = other.clone();
+
+        // find common factors of 2
+        let shift = cmp::min(
+            n.trailing_zeros(),
+            m.trailing_zeros()
+        );
+
+        // divide m and n by 2 until odd
+        // m inside loop
+        n >>= n.trailing_zeros();
+
         while !m.is_zero() {
-            let temp = m;
-            m = n % &temp;
-            n = temp;
+            m >>= m.trailing_zeros();
+            if n > m { mem::swap(&mut n, &mut m) }
+            m -= &n;
         }
-        return n;
+
+        n << shift
     }
 
     /// Calculates the Lowest Common Multiple (LCM) of the number and `other`.
@@ -1607,6 +1628,16 @@ impl BigUint {
         return self.data.len() * big_digit::BITS - zeros as usize;
     }
 
+    // self is assumed to be normalized
+    fn trailing_zeros(&self) -> usize {
+        self.data
+            .iter()
+            .enumerate()
+            .find(|&(_, &digit)| digit != 0)
+            .map(|(i, digit)| i * big_digit::BITS + digit.trailing_zeros() as usize)
+            .unwrap_or(0)
+    }
+
     /// Strips off trailing zero bigdigits - comparisons require the last element in the vector to
     /// be nonzero.
     #[inline]

src/tests/biguint.rs

@@ -1569,6 +1569,7 @@ fn test_from_str_radix() {
 
 #[test]
 fn test_all_str_radix() {
+    #[allow(unused)]
     use std::ascii::AsciiExt;
 
     let n = BigUint::new((0..10).collect());