Add stein gcd

benches/stock_functions.rs

@@ -66,16 +66,33 @@ fn bench_gcd(c: &mut Criterion) {
     let small = Gcd::new(14, 15);
 
     for input in [small, med, max].iter() {
-        group.bench_with_input(BenchmarkId::new("Binary", input.left), input, |bench, input| {
+        group.bench_with_input(
+            BenchmarkId::new("Binary", input.left),
+            input,
+            |bench, input| {
                 let a = input.left_fib;
                 let b = input.right_fib;
                 bench.iter(|| u128::gcd(black_box(a), black_box(b)))
-        });
-        group.bench_with_input(BenchmarkId::new("Euclid", input.left), input, |bench, input| {
+            },
+        );
+        group.bench_with_input(
+            BenchmarkId::new("Stein", input.left),
+            input,
+            |bench, input| {
+                let a = input.left_fib;
+                let b = input.right_fib;
+                bench.iter(|| u128::stein_gcd(black_box(a), black_box(b)))
+            },
+        );
+        group.bench_with_input(
+            BenchmarkId::new("Euclid", input.left),
+            input,
+            |bench, input| {
                 let a = input.left_fib;
                 let b = input.right_fib;
                 bench.iter(|| u128::e_gcd(black_box(a), black_box(b)))
-        });
+            },
+        );
     }
     group.finish();
 }

src/num.rs

@@ -2,6 +2,7 @@
 //!
 //! Home to the numeric trait chain of doom, aka `Unsigned`
 #![allow(unused_comparisons)]
+use core::cmp::{max, min};
 use core::convert::TryFrom;
 use core::fmt::Debug;
 use core::ops::{
@@ -113,6 +114,9 @@ pub trait Unsigned:
     /// Euclid gcd algorithm
     fn e_gcd(a: Self, b: Self) -> Self;
 
+    /// Stein gcd algorithm
+    fn stein_gcd(a: Self, b: Self) -> Self;
+
     /// Find the least common multiple of two numbers
     fn lcm(a: Self, b: Self) -> Self;
 }
@@ -221,6 +225,20 @@ macro_rules! impl_unsigned {
                     x
                 }
 
+                fn stein_gcd(a: Self, b: Self) -> Self {
+                    match ((a, b), (a & 1, b & 1)) {
+                        ((x, y), _) if x == y => y,
+                        ((0, x), _) | ((x, 0), _) => x,
+                        ((x, y), (0, 1)) | ((y, x), (1, 0)) => Self::stein_gcd(x >> 1, y),
+                        ((x, y), (0, 0)) => Self::stein_gcd(x >> 1, y >> 1) << 1,
+                        ((x, y), (1, 1)) => {
+                            let (x, y) = (min(x, y), max(x, y));
+                            Self::stein_gcd((y - x) >> 1, x)
+                        }
+                        _ => unreachable!(),
+                    }
+                }
+
                 fn lcm(a: $Type, b: $Type) -> $Type {
                     if (a == 0 && b == 0) {
                         return 0;
@@ -232,7 +250,6 @@ macro_rules! impl_unsigned {
         )*
     };
 }
-
 impl_empty!(
     Num,
     (i8, u8, i16, u16, i32, u32, i64, u64, i128, u128, isize, usize)
@@ -271,6 +288,7 @@ mod tests {
         assert_eq!(u8::gcd(2, 8), 2);
         assert_eq!(u16::gcd(36, 48), 12);
         assert_eq!(u16::e_gcd(36, 48), 12);
+        assert_eq!(u16::stein_gcd(36, 48), 12);
     }
 
     #[test]