math/big: Implemented binary GCD algorithm

benchmark                    old ns/op    new ns/op    delta
BenchmarkGCD10x10                 4383         2126  -51.49%
BenchmarkGCD10x100                5612         2124  -62.15%
BenchmarkGCD10x1000               8843         2622  -70.35%
BenchmarkGCD10x10000             17025         6576  -61.37%
BenchmarkGCD10x100000           118985        48130  -59.55%
BenchmarkGCD100x100              45328        11683  -74.23%
BenchmarkGCD100x1000             50141        12678  -74.72%
BenchmarkGCD100x10000           110314        26719  -75.78%
BenchmarkGCD100x100000          630000       156041  -75.23%
BenchmarkGCD1000x1000           654809       137973  -78.93%
BenchmarkGCD1000x10000          985683       159951  -83.77%
BenchmarkGCD1000x100000        4920792       366399  -92.55%
BenchmarkGCD10000x10000       16848950      3732062  -77.85%
BenchmarkGCD10000x100000      55401500      4675876  -91.56%
BenchmarkGCD100000x100000   1126775000    258951800  -77.02%

R=gri, rsc, bradfitz, remyoudompheng, mtj
CC=golang-dev
https://golang.org/cl/6305065

src/pkg/math/big/gcd_test.go

+// Copyright 2012 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// This file implements a GCD benchmark.
+// Usage: go test math/big -test.bench GCD
+
+package big
+
+import (
+	"math/rand"
+	"testing"
+)
+
+// randInt returns a pseudo-random Int in the range [1<<(size-1), (1<<size) - 1]
+func randInt(r *rand.Rand, size uint) *Int {
+	n := new(Int).Lsh(intOne, size-1)
+	x := new(Int).Rand(r, n)
+	return x.Add(x, n) // make sure result > 1<<(size-1)
+}
+
+func runGCD(b *testing.B, aSize, bSize uint) {
+	b.StopTimer()
+	var r = rand.New(rand.NewSource(1234))
+	aa := randInt(r, aSize)
+	bb := randInt(r, bSize)
+	b.StartTimer()
+	for i := 0; i < b.N; i++ {
+		new(Int).GCD(nil, nil, aa, bb)
+	}
+}
+
+func BenchmarkGCD10x10(b *testing.B)         { runGCD(b, 10, 10) }
+func BenchmarkGCD10x100(b *testing.B)        { runGCD(b, 10, 100) }
+func BenchmarkGCD10x1000(b *testing.B)       { runGCD(b, 10, 1000) }
+func BenchmarkGCD10x10000(b *testing.B)      { runGCD(b, 10, 10000) }
+func BenchmarkGCD10x100000(b *testing.B)     { runGCD(b, 10, 100000) }
+func BenchmarkGCD100x100(b *testing.B)       { runGCD(b, 100, 100) }
+func BenchmarkGCD100x1000(b *testing.B)      { runGCD(b, 100, 1000) }
+func BenchmarkGCD100x10000(b *testing.B)     { runGCD(b, 100, 10000) }
+func BenchmarkGCD100x100000(b *testing.B)    { runGCD(b, 100, 100000) }
+func BenchmarkGCD1000x1000(b *testing.B)     { runGCD(b, 1000, 1000) }
+func BenchmarkGCD1000x10000(b *testing.B)    { runGCD(b, 1000, 10000) }
+func BenchmarkGCD1000x100000(b *testing.B)   { runGCD(b, 1000, 100000) }
+func BenchmarkGCD10000x10000(b *testing.B)   { runGCD(b, 10000, 10000) }
+func BenchmarkGCD10000x100000(b *testing.B)  { runGCD(b, 10000, 100000) }
+func BenchmarkGCD100000x100000(b *testing.B) { runGCD(b, 100000, 100000) }

src/pkg/math/big/int.go

@@ -596,6 +596,9 @@ func (z *Int) GCD(x, y, a, b *Int) *Int {
 		}
 		return z
 	}
+	if x == nil && y == nil {
+		return z.binaryGCD(a, b)
+	}
 
 	A := new(Int).Set(a)
 	B := new(Int).Set(b)
@@ -640,6 +643,56 @@ func (z *Int) GCD(x, y, a, b *Int) *Int {
 	return z
 }
 
+// binaryGCD sets z to the greatest common divisor of a and b, which must be
+// positive, and returns z.
+// See Knuth, The Art of Computer Programming, Vol. 2, Section 4.5.2, Algorithm B.
+func (z *Int) binaryGCD(a, b *Int) *Int {
+	u := z
+	v := new(Int)
+	// use one Euclidean iteration to ensure that u and v are approx. the same size
+	switch {
+	case len(a.abs) > len(b.abs):
+		u.Set(b)
+		v.Rem(a, b)
+	case len(a.abs) < len(b.abs):
+		u.Set(a)
+		v.Rem(b, a)
+	default:
+		u.Set(a)
+		v.Set(b)
+	}
+
+	// determine largest k such that u = u' << k, v = v' << k
+	k := u.abs.trailingZeroBits()
+	if vk := v.abs.trailingZeroBits(); vk < k {
+		k = vk
+	}
+	u.Rsh(u, k)
+	v.Rsh(v, k)
+
+	// determine t (we know that u > 0)
+	t := new(Int)
+	if u.abs[0]&1 != 0 {
+		// u is odd
+		t.Neg(v)
+	} else {
+		t.Set(u)
+	}
+
+	for len(t.abs) > 0 {
+		// reduce t
+		t.Rsh(t, t.abs.trailingZeroBits())
+		if t.neg {
+			v.Neg(t)
+		} else {
+			u.Set(t)
+		}
+		t.Sub(u, v)
+	}
+
+	return u.Lsh(u, k)
+}
+
 // ProbablyPrime performs n Miller-Rabin tests to check whether x is prime.
 // If it returns true, x is prime with probability 1 - 1/4^n.
 // If it returns false, x is not prime.

src/pkg/math/big/int_test.go

@@ -860,6 +860,13 @@ func TestGcd(t *testing.T) {
 			expectedD.Cmp(d) != 0 {
 			t.Errorf("#%d got (%s %s %s) want (%s %s %s)", i, x, y, d, expectedX, expectedY, expectedD)
 		}
+
+		d.binaryGCD(a, b)
+
+		if expectedD.Cmp(d) != 0 {
+			t.Errorf("#%d got (%s) want (%s)", i, d, expectedD)
+		}
+
 	}
 
 	quick.Check(checkGcd, nil)

src/pkg/math/big/rat.go

@@ -145,20 +145,6 @@ func (x *Rat) Denom() *Int {
 	return &x.b
 }
 
-func gcd(x, y nat) nat {
-	// Euclidean algorithm.
-	var a, b nat
-	a = a.set(x)
-	b = b.set(y)
-	for len(b) != 0 {
-		var q, r nat
-		_, r = q.div(r, a, b)
-		a = b
-		b = r
-	}
-	return a
-}
-
 func (z *Rat) norm() *Rat {
 	switch {
 	case len(z.a.abs) == 0:
@@ -171,14 +157,18 @@ func (z *Rat) norm() *Rat {
 		// z is int - normalize denominator
 		z.b.abs = z.b.abs.make(0)
 	default:
-		if f := gcd(z.a.abs, z.b.abs); f.cmp(natOne) != 0 {
-			z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f)
-			z.b.abs, _ = z.b.abs.div(nil, z.b.abs, f)
+		neg := z.a.neg
+		z.a.neg = false
+		z.b.neg = false
+		if f := NewInt(0).binaryGCD(&z.a, &z.b); f.Cmp(intOne) != 0 {
+			z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f.abs)
+			z.b.abs, _ = z.b.abs.div(nil, z.b.abs, f.abs)
 			if z.b.abs.cmp(natOne) == 0 {
 				// z is int - normalize denominator
 				z.b.abs = z.b.abs.make(0)
 			}
 		}
+		z.a.neg = neg
 	}
 	return z
 }