From 1ff6e0afcbf16e327488beae4eb95b2881b2bb54 Mon Sep 17 00:00:00 2001
From: Mike Nolta <mike@nolta.net>
Date: Mon, 6 Oct 2014 22:13:36 -0400
Subject: [PATCH] faster binary gcd algorithm for 64 & 128 bit ints

---
 base/intfuncs.jl | 20 ++++++++++++++++++++
 test/numbers.jl  |  3 ++-
 2 files changed, 22 insertions(+), 1 deletion(-)

diff --git a/base/intfuncs.jl b/base/intfuncs.jl
index ac1a0a342e1ec..182066d866219 100644
--- a/base/intfuncs.jl
+++ b/base/intfuncs.jl
@@ -9,6 +9,26 @@ function gcd{T<:Integer}(a::T, b::T)
     abs(a)
 end
 
+# binary GCD (aka Stein's) algorithm
+# about 1.7x (2.1x) faster for random Int64s (Int128s)
+function gcd{T<:Union(Int64,Uint64,Int128,Uint128)}(a::T, b::T)
+    a == 0 && return abs(b)
+    b == 0 && return abs(a)
+    za = trailing_zeros(a)
+    zb = trailing_zeros(b)
+    k = min(za, zb)
+    u = abs(a >> za)
+    v = abs(b >> zb)
+    while u != v
+        if u > v
+            u, v = v, u
+        end
+        v -= u
+        v >>= trailing_zeros(v)
+    end
+    u << k
+end
+
 # explicit a==0 test is to handle case of lcm(0,0) correctly
 lcm{T<:Integer}(a::T, b::T) = a == 0 ? a : abs(a * div(b, gcd(b,a)))
 
diff --git a/test/numbers.jl b/test/numbers.jl
index c5ac9bc5e759f..6ad56990b8f3e 100644
--- a/test/numbers.jl
+++ b/test/numbers.jl
@@ -1798,7 +1798,8 @@ end
 @test 3//2 <= typemax(Int)
 
 # check gcd and related functions against GMP
-for i = -20:20, j = -20:20
+for T in (Int32,Int64), ii = -20:20, jj = -20:20
+    i::T, j::T = ii, jj
     local d = gcd(i,j)
     @test d >= 0
     @test lcm(i,j) >= 0