diff --git a/src/tests/util/numerics/CMakeLists.txt b/src/tests/util/numerics/CMakeLists.txt
index 32498b602..833e3d832 100644
--- a/src/tests/util/numerics/CMakeLists.txt
+++ b/src/tests/util/numerics/CMakeLists.txt
@@ -25,3 +25,6 @@ add_test(primes ${CMAKE_CURRENT_BINARY_DIR}/primes)
 add_executable(numeric_traits numeric_traits.cpp)
 target_link_libraries(numeric_traits ${EXTRA_LIBS})
 add_test(numeric_traits ${CMAKE_CURRENT_BINARY_DIR}/numeric_traits)
+add_executable(gcd gcd.cpp)
+target_link_libraries(gcd ${EXTRA_LIBS})
+add_test(gcd ${CMAKE_CURRENT_BINARY_DIR}/gcd)
diff --git a/src/tests/util/numerics/gcd.cpp b/src/tests/util/numerics/gcd.cpp
new file mode 100644
index 000000000..bee320c83
--- /dev/null
+++ b/src/tests/util/numerics/gcd.cpp
@@ -0,0 +1,47 @@
+/*
+Copyright (c) 2013 Microsoft Corporation. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+
+Author: Leonardo de Moura
+*/
+#include <iostream>
+#include <random>
+#include "util/test.h"
+#include "util/numerics/gcd.h"
+#include "util/numerics/mpz.h"
+using namespace lean;
+
+static void tst1(unsigned num_tests, unsigned max_val) {
+    std::mt19937 rng;
+    std::uniform_int_distribution<unsigned int> uint_dist;
+    for (unsigned i = 0; i < num_tests; i++) {
+        int val1 = uint_dist(rng) % max_val;
+        int val2 = uint_dist(rng) % max_val;
+        int g, a, b;
+        gcdext(g, a, b, val1, val2);
+        mpz _g, _a, _b;
+        gcdext(_g, _a, _b, mpz(val1), mpz(val2));
+        lean_assert_eq(_g.get_int(), g);
+        if (val1 != val2) {
+            lean_assert_eq(_a.get_int(), a);
+            lean_assert_eq(_b.get_int(), b);
+        } else {
+            if (_a.get_int() == a) {
+                lean_assert_eq(_b.get_int(), b);
+            } else {
+                lean_assert_eq(_a.get_int(), b);
+                lean_assert_eq(_b.get_int(), a);
+            }
+        }
+    }
+}
+
+int main() {
+    tst1(1000,    10);
+    tst1(1000,    100);
+    tst1(10000,   1000);
+    tst1(10000,   100);
+    tst1(10000,   1000);
+    tst1(1000000, 10000);
+    return has_violations() ? 1 : 0;
+}
diff --git a/src/util/numerics/gcd.h b/src/util/numerics/gcd.h
new file mode 100644
index 000000000..3383510c7
--- /dev/null
+++ b/src/util/numerics/gcd.h
@@ -0,0 +1,91 @@
+/*
+Copyright (c) 2013 Microsoft Corporation. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+
+Author: Leonardo de Moura
+*/
+#include <utility>
+#include "util/debug.h"
+
+namespace lean {
+/**
+   \brief Template for computing GCD using binary method
+*/
+template<typename T>
+static T gcd(T u, T v) {
+    if (u == 0)
+        return v;
+    if (v == 0)
+        return u;
+
+    int k;
+
+    for (k = 0; ((u | v) & 1) == 0; ++k) {
+        u >>= 1;
+        v >>= 1;
+    }
+
+    while ((u & 1) == 0)
+        u >>= 1;
+
+    do {
+        while ((v & 1) == 0)
+            v >>= 1;
+
+        if (u < v) {
+            v -= u;
+        } else {
+            T diff = u - v;
+            u = v;
+            v = diff;
+        }
+        v >>= 1;
+    } while (v != 0);
+
+    return u << k;
+}
+
+/**
+   \brief Template for the extended GCD procedure
+*/
+template<typename T>
+void gcdext(T & r, T & a, T & b, T const & r1, T const & r2) {
+    using std::swap;
+    T aux, quot;
+    T tmp1 = r1;
+    T tmp2 = r2;
+    a = 1;
+    b = 0;
+    T next_a, next_b;
+    next_a = 0;
+    next_b = 1;
+    if (tmp1 < 0) tmp1 = -tmp1;
+    if (tmp2 < 0) tmp2 = -tmp2;
+
+    if (tmp1 < tmp2) {
+        swap(tmp1, tmp2);
+        swap(next_a, next_b);
+        swap(a, b);
+    }
+
+    // tmp1 >= tmp2 >= 0
+    // quot_rem in one function would be faster.
+    while (tmp2 > 0) {
+        lean_assert(tmp1 >= tmp2);
+        aux    = tmp2;
+        quot   = tmp1 / tmp2;
+        tmp2   = tmp1 % tmp2;
+        tmp1   = aux;
+        aux    = next_a;
+        next_a = a - (quot * next_a);
+        a      = aux;
+        aux    = next_b;
+        next_b = b - (quot * next_b);
+        b      = aux;
+    }
+
+    if (r1 < 0) a = -a;
+    if (r2 < 0) b = -b;
+    r = tmp1;
+}
+}