math/big: implement Lehmer's extended GCD algorithm

Updates #15833

The extended GCD algorithm can be implemented using
Lehmer's algorithm with additional updates for the
cosequences following Algorithm 10.45 from Cohen et al.
"Handbook of Elliptic and Hyperelliptic Curve Cryptography" pp 192.
This brings the speed of the extended GCD calculation within
~2x of the base GCD calculation.  There is a slight degradation in
the non-extended GCD speed for small inputs (1-2 words) due to the
additional code to handle the extended updates.

name                          old time/op    new time/op    delta
GCD10x10/WithoutXY-4             262ns ± 1%     266ns ± 2%     ~     (p=0.333 n=5+5)
GCD10x10/WithXY-4               1.42µs ± 2%    0.74µs ± 3%  -47.90%  (p=0.008 n=5+5)
GCD10x100/WithoutXY-4            520ns ± 2%     539ns ± 1%   +3.81%  (p=0.008 n=5+5)
GCD10x100/WithXY-4              2.32µs ± 1%    1.67µs ± 0%  -27.80%  (p=0.008 n=5+5)
GCD10x1000/WithoutXY-4          1.40µs ± 1%    1.45µs ± 2%   +3.26%  (p=0.016 n=4+5)
GCD10x1000/WithXY-4             4.78µs ± 1%    3.43µs ± 1%  -28.37%  (p=0.008 n=5+5)
GCD10x10000/WithoutXY-4         10.0µs ± 0%    10.2µs ± 3%   +1.80%  (p=0.008 n=5+5)
GCD10x10000/WithXY-4            20.9µs ± 3%    17.9µs ± 1%  -14.20%  (p=0.008 n=5+5)
GCD10x100000/WithoutXY-4        96.8µs ± 0%    96.3µs ± 1%     ~     (p=0.310 n=5+5)
GCD10x100000/WithXY-4            196µs ± 3%     159µs ± 2%  -18.61%  (p=0.008 n=5+5)
GCD100x100/WithoutXY-4          2.53µs ±15%    2.34µs ± 0%   -7.35%  (p=0.008 n=5+5)
GCD100x100/WithXY-4             19.3µs ± 0%     3.9µs ± 1%  -79.58%  (p=0.008 n=5+5)
GCD100x1000/WithoutXY-4         4.23µs ± 0%    4.17µs ± 3%     ~     (p=0.127 n=5+5)
GCD100x1000/WithXY-4            22.8µs ± 1%     7.5µs ±10%  -67.00%  (p=0.008 n=5+5)
GCD100x10000/WithoutXY-4        19.1µs ± 0%    19.0µs ± 0%     ~     (p=0.095 n=5+5)
GCD100x10000/WithXY-4           75.1µs ± 2%    30.5µs ± 2%  -59.38%  (p=0.008 n=5+5)
GCD100x100000/WithoutXY-4        170µs ± 5%     167µs ± 1%     ~     (p=1.000 n=5+5)
GCD100x100000/WithXY-4           542µs ± 2%     267µs ± 2%  -50.79%  (p=0.008 n=5+5)
GCD1000x1000/WithoutXY-4        28.0µs ± 0%    27.1µs ± 0%   -3.29%  (p=0.008 n=5+5)
GCD1000x1000/WithXY-4            329µs ± 0%      42µs ± 1%  -87.12%  (p=0.008 n=5+5)
GCD1000x10000/WithoutXY-4       47.2µs ± 0%    46.4µs ± 0%   -1.65%  (p=0.016 n=5+4)
GCD1000x10000/WithXY-4           607µs ± 9%     123µs ± 1%  -79.70%  (p=0.008 n=5+5)
GCD1000x100000/WithoutXY-4       260µs ±17%     245µs ± 0%     ~     (p=0.056 n=5+5)
GCD1000x100000/WithXY-4         3.64ms ± 1%    0.93ms ± 1%  -74.41%  (p=0.016 n=4+5)
GCD10000x10000/WithoutXY-4       513µs ± 0%     507µs ± 0%   -1.22%  (p=0.008 n=5+5)
GCD10000x10000/WithXY-4         7.44ms ± 1%    1.00ms ± 0%  -86.58%  (p=0.008 n=5+5)
GCD10000x100000/WithoutXY-4     1.23ms ± 0%    1.23ms ± 1%     ~     (p=0.056 n=5+5)
GCD10000x100000/WithXY-4        37.3ms ± 0%     7.3ms ± 1%  -80.45%  (p=0.008 n=5+5)
GCD100000x100000/WithoutXY-4    24.2ms ± 0%    24.2ms ± 0%     ~     (p=0.841 n=5+5)
GCD100000x100000/WithXY-4        505ms ± 1%      56ms ± 1%  -88.92%  (p=0.008 n=5+5)

Change-Id: I25f42ab8c55033acb83cc32bb03c12c1963925e8
Reviewed-on: https://go-review.googlesource.com/78755Reviewed-by: Robert Griesemer <gri@golang.org>
Run-TryBot: Robert Griesemer <gri@golang.org>
TryBot-Result: Gobot Gobot <gobot@golang.org>

src/math/big/int.go

@@ -485,162 +485,223 @@ func (z *Int) GCD(x, y, a, b *Int) *Int {
 		}
 		return z
 	}
-	if x == nil && y == nil {
-		return z.lehmerGCD(a, b)
-	}
 
-	A := new(Int).Set(a)
-	B := new(Int).Set(b)
+	return z.lehmerGCD(x, y, a, b)
+}
 
-	X := new(Int)
-	lastX := new(Int).SetInt64(1)
+// lehmerSimulate attempts to simulate several Euclidean update steps
+// using the leading digits of A and B.  It returns u0, u1, v0, v1
+// such that A and B can be updated as:
+//		A = u0*A + v0*B
+//		B = u1*A + v1*B
+// Requirements: A >= B and len(B.abs) >= 2
+// 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
+func lehmerSimulate(A, B *Int) (u0, u1, v0, v1 Word, even bool) {
+	// initialize the digits
+	var a1, a2, u2, 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
+	}
 
-	q := new(Int)
-	temp := new(Int)
+	// 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.
+	// Note that overflow of a Word is not possible when computing the remainder
+	// sequence and cosequences since the cosequence size is bounded by the input size.
+	// See section 4.2 of Jebelean for details.
+	for a2 >= v2 && a1-a2 >= v1+v2 {
+		q, r := a1/a2, a1%a2
+		a1, a2 = a2, r
+		u0, u1, u2 = u1, u2, u1+q*u2
+		v0, v1, v2 = v1, v2, v1+q*v2
+		even = !even
+	}
+	return
+}
 
-	r := new(Int)
-	for len(B.abs) > 0 {
-		q, r = q.QuoRem(A, B, r)
+// lehmerUpdate updates the inputs A and B such that:
+//		A = u0*A + v0*B
+//		B = u1*A + v1*B
+// where the signs of u0, u1, v0, v1 are given by even
+// For even == true: u0, v1 >= 0 && u1, v0 <= 0
+// For even == false: u0, v1 <= 0 && u1, v0 >= 0
+// q, r, s, t are temporary variables to avoid allocations in the multiplication
+func lehmerUpdate(A, B, q, r, s, t *Int, u0, u1, v0, v1 Word, even bool) {
+
+	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)
+	q.abs = q.abs.setWord(v1)
+	r.neg = even
+	q.neg = !even
+
+	r.Mul(A, r)
+	q.Mul(B, q)
+
+	A.Add(t, s)
+	B.Add(r, q)
+}
 
-		A, B, r = B, r, A
+// euclidUpdate performs a single step of the Euclidean GCD algorithm
+// if extended is true, it also updates the cosequence Ua, Ub
+func euclidUpdate(A, B, Ua, Ub, q, r, s, t *Int, extended bool) {
+	q, r = q.QuoRem(A, B, r)
 
-		temp.Set(X)
-		X.Mul(X, q)
-		X.Sub(lastX, X)
-		lastX.Set(temp)
-	}
+	*A, *B, *r = *B, *r, *A
 
-	if x != nil {
-		*x = *lastX
+	if extended {
+		// Ua, Ub = Ub, Ua - q*Ub
+		t.Set(Ub)
+		s.Mul(Ub, q)
+		Ub.Sub(Ua, s)
+		Ua.Set(t)
 	}
-
-	if y != nil {
-		// y = (z - a*x)/b
-		y.Mul(a, lastX)
-		y.Sub(A, y)
-		y.Div(y, b)
-	}
-
-	*z = *A
-	return z
 }
 
 // lehmerGCD sets z to the greatest common divisor of a and b,
 // which both must be > 0, and returns z.
+// If x or y are not nil, their values are set such that z = a*x + b*y.
 // 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 {
-	// ensure a >= b
-	if a.abs.cmp(b.abs) < 0 {
-		a, b = b, a
-	}
+// The cosequences are updated according to Algorithm 10.45 from
+// Cohen et al. "Handbook of Elliptic and Hyperelliptic Curve Cryptography" pp 192.
+func (z *Int) lehmerGCD(x, y, a, b *Int) *Int {
+	var A, B, Ua, Ub *Int
+
+	A = new(Int).Set(a)
+	B = new(Int).Set(b)
+
+	extended := x != nil || y != nil
 
-	// 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)
+	if extended {
+		// Ua (Ub) tracks how many times input a has been accumulated into A (B).
+		Ua = new(Int).SetInt64(1)
+		Ub = new(Int)
+	}
 
 	// temp variables for multiprecision update
-	t := new(Int)
+	q := new(Int)
 	r := new(Int)
 	s := new(Int)
-	w := new(Int)
+	t := new(Int)
+
+	// ensure A >= B
+	if A.abs.cmp(B.abs) < 0 {
+		A, B = B, A
+		Ub, Ua = Ua, Ub
+	}
 
 	// 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
-		}
+		// Attempt to calculate in single-precision using leading words of A and B.
+		u0, u1, v0, v1, even := lehmerSimulate(A, B)
 
-		// 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.
-		// Note that overflow of a Word is not possible when computing the remainder
-		// sequence and cosequences since the cosequence size is bounded by the input size.
-		// See section 4.2 of Jebelean for details.
-		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
+		// multiprecision Step
 		if v0 != 0 {
-			// simulate the effect of the single precision steps using the cosequences
+			// Simulate the effect of the single-precision steps using the cosequences.
 			// A = u0*A + v0*B
 			// B = u1*A + v1*B
+			lehmerUpdate(A, B, q, r, s, t, u0, u1, v0, v1, even)
 
-			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)
+			if extended {
+				// Ua = u0*Ua + v0*Ub
+				// Ub = u1*Ua + v1*Ub
+				lehmerUpdate(Ua, Ub, q, r, s, t, u0, u1, v0, v1, even)
+			}
 
 		} else {
-			// single-digit calculations failed to simulate any quotients
-			// do a standard Euclidean step
-			t.Rem(A, B)
-			A, B, t = B, t, A
+			// Single-digit calculations failed to simulate any quotients.
+			// Do a standard Euclidean step.
+			euclidUpdate(A, B, Ua, Ub, q, r, s, t, extended)
 		}
 	}
 
 	if len(B.abs) > 0 {
-		// standard Euclidean algorithm base case for B a single Word
+		// extended Euclidean algorithm base case if B is 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
+			// A is longer than a single Word, so one update is needed.
+			euclidUpdate(A, B, Ua, Ub, q, r, s, t, extended)
 		}
 		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 and B are both a single Word.
+			aWord, bWord := A.abs[0], B.abs[0]
+			if extended {
+				var ua, ub, va, vb Word
+				ua, ub = 1, 0
+				va, vb = 0, 1
+				even := true
+				for bWord != 0 {
+					q, r := aWord/bWord, aWord%bWord
+					aWord, bWord = bWord, r
+					ua, ub = ub, ua+q*ub
+					va, vb = vb, va+q*vb
+					even = !even
+				}
+
+				t.abs = t.abs.setWord(ua)
+				s.abs = s.abs.setWord(va)
+				t.neg = !even
+				s.neg = even
+
+				t.Mul(Ua, t)
+				s.Mul(Ub, s)
+
+				Ua.Add(t, s)
+			} else {
+				for bWord != 0 {
+					aWord, bWord = bWord, aWord%bWord
+				}
 			}
-			A.abs[0] = a1
+			A.abs[0] = aWord
 		}
 	}
+
+	if x != nil {
+		*x = *Ua
+	}
+
+	if y != nil {
+		// y = (z - a*x)/b
+		y.Mul(a, Ua)
+		y.Sub(A, y)
+		y.Div(y, b)
+	}
+
 	*z = *A
+
 	return z
 }
 

src/math/big/int_test.go

@@ -675,21 +675,43 @@ func checkGcd(aBytes, bBytes []byte) bool {
 	return x.Cmp(d) == 0
 }
 
-// euclidGCD is a reference implementation of Euclid's GCD
-// algorithm for testing against optimized algorithms.
+// euclidExtGCD is a reference implementation of Euclid's
+// extended GCD algorithm for testing against optimized algorithms.
 // Requirements: a, b > 0
-func euclidGCD(a, b *Int) *Int {
-
+func euclidExtGCD(a, b *Int) (g, x, y *Int) {
 	A := new(Int).Set(a)
 	B := new(Int).Set(b)
-	t := new(Int)
 
+	// A = Ua*a + Va*b
+	// B = Ub*a + Vb*b
+	Ua := new(Int).SetInt64(1)
+	Va := new(Int)
+
+	Ub := new(Int)
+	Vb := new(Int).SetInt64(1)
+
+	q := new(Int)
+	temp := new(Int)
+
+	r := new(Int)
 	for len(B.abs) > 0 {
-		// A, B = B, A % B
-		t.Rem(A, B)
-		A, B, t = B, t, A
+		q, r = q.QuoRem(A, B, r)
+
+		A, B, r = B, r, A
+
+		// Ua, Ub = Ub, Ua-q*Ub
+		temp.Set(Ub)
+		Ub.Mul(Ub, q)
+		Ub.Sub(Ua, Ub)
+		Ua.Set(temp)
+
+		// Va, Vb = Vb, Va-q*Vb
+		temp.Set(Vb)
+		Vb.Mul(Vb, q)
+		Vb.Sub(Va, Vb)
+		Va.Set(temp)
 	}
-	return A
+	return A, Ua, Va
 }
 
 func checkLehmerGcd(aBytes, bBytes []byte) bool {
@@ -700,12 +722,28 @@ func checkLehmerGcd(aBytes, bBytes []byte) bool {
 		return true // can only test positive arguments
 	}
 
-	d := new(Int).lehmerGCD(a, b)
-	d0 := euclidGCD(a, b)
+	d := new(Int).lehmerGCD(nil, nil, a, b)
+	d0, _, _ := euclidExtGCD(a, b)
 
 	return d.Cmp(d0) == 0
 }
 
+func checkLehmerExtGcd(aBytes, bBytes []byte) bool {
+	a := new(Int).SetBytes(aBytes)
+	b := new(Int).SetBytes(bBytes)
+	x := new(Int)
+	y := new(Int)
+
+	if a.Sign() <= 0 || b.Sign() <= 0 {
+		return true // can only test positive arguments
+	}
+
+	d := new(Int).lehmerGCD(x, y, a, b)
+	d0, x0, y0 := euclidExtGCD(a, b)
+
+	return d.Cmp(d0) == 0 && x.Cmp(x0) == 0 && y.Cmp(y0) == 0
+}
+
 var gcdTests = []struct {
 	d, x, y, a, b string
 }{
@@ -785,6 +823,10 @@ func TestGcd(t *testing.T) {
 	if err := quick.Check(checkLehmerGcd, nil); err != nil {
 		t.Error(err)
 	}
+
+	if err := quick.Check(checkLehmerExtGcd, nil); err != nil {
+		t.Error(err)
+	}
 }
 
 type intShiftTest struct {

src/math/big/rat.go

@@ -422,7 +422,7 @@ func (z *Rat) norm() *Rat {
 		neg := z.a.neg
 		z.a.neg = false
 		z.b.neg = false
-		if f := NewInt(0).lehmerGCD(&z.a, &z.b); f.Cmp(intOne) != 0 {
+		if f := NewInt(0).lehmerGCD(nil, nil, &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 {