From: Christopher Swenson Date: Wed, 13 Jun 2012 16:31:20 +0000 (-0700) Subject: math/big: Implemented binary GCD algorithm X-Git-Tag: go1.1rc2~2938 X-Git-Url: http://www.git.cypherpunks.su/?a=commitdiff_plain;h=38735b957c2e7d1f93021943dafc5b931b9ccab3;p=gostls13.git 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 --- diff --git a/src/pkg/math/big/gcd_test.go b/src/pkg/math/big/gcd_test.go new file mode 100644 index 0000000000..c0b9f58300 --- /dev/null +++ b/src/pkg/math/big/gcd_test.go @@ -0,0 +1,47 @@ +// 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< 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) } diff --git a/src/pkg/math/big/int.go b/src/pkg/math/big/int.go index 276f56708a..16fd9bfa98 100644 --- a/src/pkg/math/big/int.go +++ b/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. diff --git a/src/pkg/math/big/int_test.go b/src/pkg/math/big/int_test.go index 30e55916b5..4ec2ac56f3 100644 --- a/src/pkg/math/big/int_test.go +++ b/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) diff --git a/src/pkg/math/big/rat.go b/src/pkg/math/big/rat.go index 5c2a48654a..eccf34e482 100644 --- a/src/pkg/math/big/rat.go +++ b/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 }