faster binary gcd algorithm for 64 & 128 bit ints

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)))
 

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