return z
}
if x == nil && y == nil {
- return z.binaryGCD(a, b)
+ return z.lehmerGCD(a, b)
}
A := new(Int).Set(a)
return z
}
-// binaryGCD sets z to the greatest common divisor of a and b, which both must
-// be > 0, 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):
- // must set v before u since u may be alias for a or b (was issue #11284)
- v.Rem(a, b)
- u.Set(b)
- case len(a.abs) < len(b.abs):
- v.Rem(b, a)
- u.Set(a)
- default:
- v.Set(b)
- u.Set(a)
- }
- // a, b must not be used anymore (may be aliases with u)
+// lehmerGCD sets z to the greatest common divisor of a and b,
+// which both must be > 0, and returns z.
+// See Knuth, The Art of Computer Programming, Vol. 2, Section 4.5.2, Algorithm L.
+// This implementation uses the improved condition by Collins requiring only one
+// quotient and avoiding the possibility of single Word overflow.
+// See Jebelean, "Improving the multiprecision Euclidean algorithm",
+// Design and Implementation of Symbolic Computation Systems, pp 45-58.
+func (z *Int) lehmerGCD(a, b *Int) *Int {
- // v might be 0 now
- if len(v.abs) == 0 {
- return u
+ // ensure a >= b
+ if a.abs.cmp(b.abs) < 0 {
+ a, b = b, a
}
- // u > 0 && v > 0
- // 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)
+ // don't destroy incoming values of a and b
+ B := new(Int).Set(b) // must be set first in case b is an alias of z
+ A := z.Set(a)
- // determine t (we know that u > 0)
+ // temp variables for multiprecision update
t := new(Int)
- if u.abs[0]&1 != 0 {
- // u is odd
- t.Neg(v)
- } else {
- t.Set(u)
- }
+ r := new(Int)
+ s := new(Int)
+ w := new(Int)
+
+ // loop invariant A >= B
+ for len(B.abs) > 1 {
+
+ // initialize the digits
+ var a1, a2, u0, u1, u2, v0, v1, v2 Word
+
+ m := len(B.abs) // m >= 2
+ n := len(A.abs) // n >= m >= 2
+
+ // extract the top Word of bits from A and B
+ h := nlz(A.abs[n-1])
+ a1 = (A.abs[n-1] << h) | (A.abs[n-2] >> (_W - h))
+ // B may have implicit zero words in the high bits if the lengths differ
+ switch {
+ case n == m:
+ a2 = (B.abs[n-1] << h) | (B.abs[n-2] >> (_W - h))
+ case n == m+1:
+ a2 = (B.abs[n-2] >> (_W - h))
+ default:
+ a2 = 0
+ }
+
+ // Since we are calculating with full words to avoid overflow,
+ // we use 'even' to track the sign of the cosequences.
+ // For even iterations: u0, v1 >= 0 && u1, v0 <= 0
+ // For odd iterations: u0, v1 <= 0 && u1, v0 >= 0
+ // The first iteration starts with k=1 (odd).
+ even := false
+ // variables to track the cosequences
+ u0, u1, u2 = 0, 1, 0
+ v0, v1, v2 = 0, 0, 1
+
+ // calculate the quotient and cosequences using Collins' stopping condition
+ for a2 >= v2 && a1-a2 >= v1+v2 {
+ q := a1 / a2
+ a1, a2 = a2, a1-q*a2
+ u0, u1, u2 = u1, u2, u1+q*u2
+ v0, v1, v2 = v1, v2, v1+q*v2
+ even = !even
+ }
+
+ // multiprecision step
+ if v0 != 0 {
+ // simulate the effect of the single precision steps using the cosequences
+ // A = u0*A + v0*B
+ // B = u1*A + v1*B
+
+ t.abs = t.abs.setWord(u0)
+ s.abs = s.abs.setWord(v0)
+ t.neg = !even
+ s.neg = even
+
+ t.Mul(A, t)
+ s.Mul(B, s)
+
+ r.abs = r.abs.setWord(u1)
+ w.abs = w.abs.setWord(v1)
+ r.neg = even
+ w.neg = !even
+
+ r.Mul(A, r)
+ w.Mul(B, w)
+
+ A.Add(t, s)
+ B.Add(r, w)
- for len(t.abs) > 0 {
- // reduce t
- t.Rsh(t, t.abs.trailingZeroBits())
- if t.neg {
- v, t = t, v
- v.neg = len(v.abs) > 0 && !v.neg // 0 has no sign
} else {
- u, t = t, u
+ // single-digit calculations failed to simluate any quotients
+ // do a standard Euclidean step
+ t.Rem(A, B)
+ A, B, t = B, t, A
}
- t.Sub(u, v)
}
- return z.Lsh(u, k)
+ if len(B.abs) > 0 {
+ // standard Euclidean algorithm base case for B a single Word
+ if len(A.abs) > 1 {
+ // A is longer than a single Word
+ t.Rem(A, B)
+ A, B, t = B, t, A
+ }
+ if len(B.abs) > 0 {
+ // A and B are both a single Word
+ a1, a2 := A.abs[0], B.abs[0]
+ for a2 != 0 {
+ a1, a2 = a2, a1%a2
+ }
+ A.abs[0] = a1
+ }
+ }
+ *z = *A
+ return z
}
// Rand sets z to a pseudo-random number in [0, n) and returns z.
return x.Cmp(d) == 0
}
+// euclidGCD is a reference implementation of Euclid's GCD
+// algorithm for testing against optimized algorithms.
+// Requirements: a, b > 0
+func euclidGCD(a, b *Int) *Int {
+
+ A := new(Int).Set(a)
+ B := new(Int).Set(b)
+ t := new(Int)
+
+ for len(B.abs) > 0 {
+ // A, B = B, A % B
+ t.Rem(A, B)
+ A, B, t = B, t, A
+ }
+ return A
+}
+
+func checkLehmerGcd(aBytes, bBytes []byte) bool {
+ a := new(Int).SetBytes(aBytes)
+ b := new(Int).SetBytes(bBytes)
+
+ if a.Sign() <= 0 || b.Sign() <= 0 {
+ return true // can only test positive arguments
+ }
+
+ d := new(Int).lehmerGCD(a, b)
+ d0 := euclidGCD(a, b)
+
+ return d.Cmp(d0) == 0
+}
+
var gcdTests = []struct {
d, x, y, a, b string
}{
D := new(Int).GCD(X, Y, a, b)
if D.Cmp(d) != 0 {
- t.Errorf("GCD(%s, %s): got d = %s, want %s", a, b, D, d)
+ t.Errorf("GCD(%s, %s, %s, %s): got d = %s, want %s", x, y, a, b, D, d)
}
if x != nil && X.Cmp(x) != 0 {
- t.Errorf("GCD(%s, %s): got x = %s, want %s", a, b, X, x)
+ t.Errorf("GCD(%s, %s, %s, %s): got x = %s, want %s", x, y, a, b, X, x)
}
if y != nil && Y.Cmp(y) != 0 {
- t.Errorf("GCD(%s, %s): got y = %s, want %s", a, b, Y, y)
- }
-
- // binaryGCD requires a > 0 && b > 0
- if a.Sign() <= 0 || b.Sign() <= 0 {
- return
- }
-
- D.binaryGCD(a, b)
- if D.Cmp(d) != 0 {
- t.Errorf("binaryGcd(%s, %s): got d = %s, want %s", a, b, D, d)
+ t.Errorf("GCD(%s, %s, %s, %s): got y = %s, want %s", x, y, a, b, Y, y)
}
// check results in presence of aliasing (issue #11284)
a2 := new(Int).Set(a)
b2 := new(Int).Set(b)
- a2.binaryGCD(a2, b2) // result is same as 1st argument
+ a2.GCD(X, Y, a2, b2) // result is same as 1st argument
if a2.Cmp(d) != 0 {
- t.Errorf("binaryGcd(%s, %s): got d = %s, want %s", a, b, a2, d)
+ t.Errorf("aliased z = a GCD(%s, %s, %s, %s): got d = %s, want %s", x, y, a, b, a2, d)
}
a2 = new(Int).Set(a)
b2 = new(Int).Set(b)
- b2.binaryGCD(a2, b2) // result is same as 2nd argument
+ b2.GCD(X, Y, a2, b2) // result is same as 2nd argument
if b2.Cmp(d) != 0 {
- t.Errorf("binaryGcd(%s, %s): got d = %s, want %s", a, b, b2, d)
+ t.Errorf("aliased z = b GCD(%s, %s, %s, %s): got d = %s, want %s", x, y, a, b, b2, d)
}
}
if err := quick.Check(checkGcd, nil); err != nil {
t.Error(err)
}
+
+ if err := quick.Check(checkLehmerGcd, nil); err != nil {
+ t.Error(err)
+ }
}
type intShiftTest struct {
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 {
+ if f := NewInt(0).lehmerGCD(&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 {