diff --git a/src/hazmat/gcd.rs b/src/hazmat/gcd.rs
index 57a0c85..d7ea864 100644
--- a/src/hazmat/gcd.rs
+++ b/src/hazmat/gcd.rs
@@ -3,12 +3,8 @@ 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.
-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>`
-    // (since `gcd(n, 0) = n`).
-    debug_assert!(m != 0);
-
+pub(crate) fn gcd_vartime<T: Integer>(n: &T, m: NonZero<Word>) -> Word {
+    let m = m.get();
     // This we can check since it doesn't affect the return type,
     // even though `n` will not be 0 either in the application.
     if n.is_zero().into() {
@@ -16,7 +12,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,21 +27,54 @@ 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
+fn binary_gcd(mut n: Word, mut m: Word) -> Word {
+    // 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();
     }
 }
 
 #[cfg(test)]
 mod tests {
-    use crypto_bigint::{Word, U128};
+    use crypto_bigint::{NonZero, Word, U128};
     use num_bigint::BigUint;
     use num_integer::Integer;
     use proptest::prelude::*;
@@ -54,11 +83,20 @@ mod tests {
 
     #[test]
     fn corner_cases() {
-        assert_eq!(gcd_vartime(&U128::from(0u64), 5), 5);
-        assert_eq!(gcd_vartime(&U128::from(1u64), 11 * 13 * 19), 1);
-        assert_eq!(gcd_vartime(&U128::from(7u64 * 11 * 13), 1), 1);
+        assert_eq!(gcd_vartime(&U128::from(0u64), NonZero::new(5).unwrap()), 5);
+        assert_eq!(
+            gcd_vartime(&U128::from(1u64), NonZero::new(11 * 13 * 19).unwrap()),
+            1
+        );
+        assert_eq!(
+            gcd_vartime(&U128::from(7u64 * 11 * 13), NonZero::new(1).unwrap()),
+            1
+        );
         assert_eq!(
-            gcd_vartime(&U128::from(7u64 * 11 * 13), 11 * 13 * 19),
+            gcd_vartime(
+                &U128::from(7u64 * 11 * 13),
+                NonZero::new(11 * 13 * 19).unwrap()
+            ),
             11 * 13
         );
     }
@@ -80,7 +118,7 @@ mod tests {
             let n_bi = BigUint::from_bytes_be(n.to_be_bytes().as_ref());
             let gcd_ref: Word = n_bi.gcd(&m_bi).try_into().unwrap();
 
-            let gcd_test = gcd_vartime(&n, m);
+            let gcd_test = gcd_vartime(&n, NonZero::new(m).unwrap());
             assert_eq!(gcd_test, gcd_ref);
         }
     }
diff --git a/src/hazmat/lucas.rs b/src/hazmat/lucas.rs
index f4ad47b..d8a54a5 100644
--- a/src/hazmat/lucas.rs
+++ b/src/hazmat/lucas.rs
@@ -1,5 +1,5 @@
 //! Lucas primality test.
-use crypto_bigint::{Integer, Limb, Monty, Odd, Square, Word};
+use crypto_bigint::{Integer, Limb, Monty, NonZero, Odd, Square, Word};
 
 use super::{
     gcd::gcd_vartime,
@@ -334,7 +334,10 @@ pub fn lucas_test<T: Integer>(
     // we check that gcd(n, Q) = 1 anyway - again, since `Q` is small,
     // it does not noticeably affect the performance.
     if abs_q != 1
-        && gcd_vartime(candidate.as_ref(), abs_q) != 1
+        && gcd_vartime(
+            candidate.as_ref(),
+            NonZero::new(abs_q).expect("q is not zero by construction"),
+        ) != 1
         && candidate.as_ref() > &to_integer(abs_q)
     {
         return Primality::Composite;