--- /dev/null
+// 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) }
}
return z
}
+ if x == nil && y == nil {
+ return z.binaryGCD(a, b)
+ }
A := new(Int).Set(a)
B := new(Int).Set(b)
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.
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:
// 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
}