diff --git a/src/math/greatest_common_divisor.rs b/src/math/greatest_common_divisor.rs
index af1f44d8009..8c88434bade 100644
--- a/src/math/greatest_common_divisor.rs
+++ b/src/math/greatest_common_divisor.rs
@@ -4,6 +4,7 @@
 ///
 /// Wikipedia reference: https://en.wikipedia.org/wiki/Greatest_common_divisor
 /// gcd(a, b) = gcd(a, -b) = gcd(-a, b) = gcd(-a, -b) by definition of divisibility
+use std::cmp::{max, min};
 
 pub fn greatest_common_divisor_recursive(a: i64, b: i64) -> i64 {
     if a == 0 {
@@ -22,6 +23,28 @@ pub fn greatest_common_divisor_iterative(mut a: i64, mut b: i64) -> i64 {
     b.abs()
 }
 
+pub fn greatest_common_divisor_stein(a: u64, b: u64) -> u64 {
+    match ((a, b), (a & 1, b & 1)) {
+        // gcd(x, x) = x
+        ((x, y), _) if x == y => y,
+        // gcd(x, 0) = gcd(0, x) = x
+        ((0, x), _) | ((x, 0), _) => x,
+        // gcd(x, y) = gcd(x / 2, y) if x is even and y is odd
+        // gcd(x, y) = gcd(x, y / 2) if y is even and x is odd
+        ((x, y), (0, 1)) | ((y, x), (1, 0)) => greatest_common_divisor_stein(x >> 1, y),
+        // gcd(x, y) = 2 * gcd(x / 2, y / 2) if x and y are both even
+        ((x, y), (0, 0)) => greatest_common_divisor_stein(x >> 1, y >> 1) << 1,
+        // if x and y are both odd
+        ((x, y), (1, 1)) => {
+            // then gcd(x, y) = gcd((x - y) / 2, y)  if x >= y
+            //      gcd(x, y) = gcd((y - x) / 2, x)  otherwise
+            let (x, y) = (min(x, y), max(x, y));
+            greatest_common_divisor_stein((y - x) >> 1, x)
+        }
+        _ => unreachable!(),
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -44,6 +67,15 @@ mod tests {
         assert_eq!(greatest_common_divisor_iterative(27, 12), 3);
     }
 
+    #[test]
+    fn positive_number_stein() {
+        assert_eq!(greatest_common_divisor_stein(4, 16), 4);
+        assert_eq!(greatest_common_divisor_stein(16, 4), 4);
+        assert_eq!(greatest_common_divisor_stein(3, 5), 1);
+        assert_eq!(greatest_common_divisor_stein(40, 40), 40);
+        assert_eq!(greatest_common_divisor_stein(27, 12), 3);
+    }
+
     #[test]
     fn negative_number_recursive() {
         assert_eq!(greatest_common_divisor_recursive(-32, -8), 8);
diff --git a/src/math/mod.rs b/src/math/mod.rs
index cbd205c89d6..ce5e7e183e3 100644
--- a/src/math/mod.rs
+++ b/src/math/mod.rs
@@ -58,6 +58,7 @@ pub use self::gaussian_elimination::gaussian_elimination;
 pub use self::gcd_of_n_numbers::gcd;
 pub use self::greatest_common_divisor::{
     greatest_common_divisor_iterative, greatest_common_divisor_recursive,
+    greatest_common_divisor_stein,
 };
 pub use self::interest::{compound_interest, simple_interest};
 pub use self::karatsuba_multiplication::multiply;