std/math/big: optimize lehmerGCD

Author: mertcandav

std/math/big/int.jule

@@ -545,8 +545,6 @@ impl Int {
 		// temp variables for multiprecision update
 		mut q := Int{}
 		mut r := Int{}
-		mut s := Int{}
-		mut t := Int{}
 
 		// ensure A >= B
 		if cmpW(A.abs, B.abs) < 0 {
@@ -563,24 +561,25 @@ impl Int {
 				// 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)
+				lehmerUpdate(&A, &B, &q, &r, u0, u1, v0, v1, even)
+
 				if extended {
 					// Ua = u0*Ua + v0*Ub
 					// Ub = u1*Ua + v1*Ub
-					lehmerUpdate(&Ua, &Ub, &q, &r, &s, &t, u0, u1, v0, v1, even)
+					lehmerUpdate(&Ua, &Ub, &q, &r, u0, u1, v0, v1, even)
 				}
 			} else {
 				// Single-digit calculations failed to simulate any quotients.
 				// Do a standard Euclidean step.
-				euclidUpdate(&A, &B, &Ua, &Ub, &q, &r, &s, &t, extended)
+				euclidUpdate(&A, &B, &Ua, &Ub, &q, &r, extended)
 			}
 		}
 
 		if len(B.abs) > 0 {
 			// extended Euclidean algorithm base case if B is a single Word
 			if len(A.abs) > 1 {
 				// A is longer than a single Word, so one update is needed.
-				euclidUpdate(&A, &B, &Ua, &Ub, &q, &r, &s, &t, extended)
+				euclidUpdate(&A, &B, &Ua, &Ub, &q, &r, extended)
 			}
 			if len(B.abs) > 0 {
 				// A and B are both a single Word.
@@ -601,15 +600,9 @@ impl Int {
 						even = !even
 					}
 
-					setWV(&t.abs, ua)
-					setWV(&s.abs, va)
-					t.neg = !even
-					s.neg = even
-
-					t.Mul(&Ua, &t)
-					s.Mul(&Ub, &s)
-
-					Ua.Add(&t, &s)
+					mulWI(&Ua, &Ua, !even, ua)
+					mulWI(&Ub, &Ub, even, va)
+					Ua.Add(&Ua, &Ub)
 				} else {
 					for bWord != 0 {
 						aWord, bWord = bWord, aWord%bWord
@@ -1162,41 +1155,33 @@ fn lehmerSimulate(A: Int, B: Int): (u0: Word, u1: Word, v0: Word, v1: Word, 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.
-fn lehmerUpdate(mut &A: *Int, mut &B: *Int, mut &q: *Int, mut &r: *Int, mut &s: *Int, mut &t: *Int,
+fn lehmerUpdate(mut &A: *Int, mut &B: *Int, mut &q: *Int, mut &r: *Int,
 	u0: Word, u1: Word, v0: Word, v1: Word, even: bool) {
-	setWV(&t.abs, u0)
-	setWV(&s.abs, v0)
-	t.neg = !even
-	s.neg = even
-
-	t.Mul(A, t)
-	s.Mul(B, s)
-
-	setWV(&r.abs, u1)
-	setWV(&q.abs, v1)
-	r.neg = even
-	q.neg = !even
-
-	r.Mul(A, r)
-	q.Mul(B, q)
+	mulWI(q, B, even, v0)
+	mulWI(r, A, even, u1)
+	mulWI(A, A, !even, u0)
+	mulWI(B, B, !even, v1)
+	A.Add(A, q)
+	B.Add(B, r)
+}
 
-	A.Add(t, s)
-	B.Add(r, q)
+fn mulWI(mut &z: *Int, mut &x: *Int, neg: bool, w: Word) {
+	mulAddWW(&z.abs, x.abs, w, 0)
+	z.neg = x.neg != neg
 }
 
 // Performs a single step of the Euclidean GCD algorithm
 // if extended is true, it also updates the cosequence Ua, Ub.
 fn euclidUpdate(mut &A: *Int, mut &B: *Int, mut &Ua: *Int, mut &Ub: *Int, mut &q: *Int,
-	mut &r: *Int, mut &s: *Int, mut &t: *Int, extended: bool) {
+	mut &r: *Int, extended: bool) {
 	q.QuoRem(A, B, r)
 
-	*A, *B, *r = *B, *r, *A
-
 	if extended {
-		// Ua, Ub = Ub, Ua - q*Ub
-		t.Set(Ub)
-		s.Mul(Ub, q)
-		Ub.Sub(Ua, s)
-		Ua.Set(t)
+		// Ua, Ub = Ub, Ua-q*Ub
+		q.Mul(q, Ub)
+		*Ua, *Ub = *Ub, *Ua
+		Ub.Sub(Ub, q)
 	}
+
+	*A, *B, *r = *B, *r, *A
 }

std/math/big/worddiv.jule

@@ -278,7 +278,7 @@ a full wide digit for the last step:
 	                       ---------------------
 	                               [    rem    ]
 
-Today, the code in divRecStep works like the first example. Perhaps in
+Today, the code in divRecursiveStep works like the first example. Perhaps in
 the future we will make it work like the alternative, to avoid a special case
 in the final iteration.
 
@@ -621,7 +621,7 @@ fn divLargeW(mut &z: *[]Word, mut &u: *[]Word, uIn: []Word, vIn: []Word) {
 	if n < divRecursiveThreshold {
 		divBasic(z, u, v)
 	} else {
-		divRec(z, u, v)
+		divRecursive(z, u, v)
 	}
 
 	*z = normW(*z)
@@ -631,7 +631,7 @@ fn divLargeW(mut &z: *[]Word, mut &u: *[]Word, uIn: []Word, vIn: []Word) {
 	*u = normW(*u)
 }
 
-// The number of divisor digits at which point divRec is faster than divBasic.
+// The number of divisor digits at which point divRecursive is faster than divBasic.
 const divRecursiveThreshold = 100
 
 // Reports whether the two digit numbers x1 x2 > y1 y2.
@@ -685,16 +685,9 @@ fn divBasic(mut &q: *[]Word, mut &u: *[]Word, v: []Word) {
 				if rhat < prevRhat {
 					break
 				}
-				// No need for a full mulWW: x1, x2 = mulWW(qhat, vn2)
+				// TODO: No need for a full mulWW.
 				// x2 += vn2; if x2 overflows, x1++
-				// TODO: Make sure this is a good exchange with mulWW version.
-				{
-					_x2, carry := bits::Add(uint(x2), uint(vn2), 0)
-					if carry == 0 {
-						x1++
-					}
-					x2 = Word(_x2)
-				}
+				x1, x2 = mulWW(qhat, vn2)
 			}
 		}
 
@@ -732,8 +725,8 @@ fn divBasic(mut &q: *[]Word, mut &u: *[]Word, v: []Word) {
 // It overwrites z with ⌊u/v⌋ and overwrites u with the remainder r.
 // z must be large enough to hold ⌊u/v⌋.
 // This function is just for allocating and freeing temporaries
-// around divRecStep, the real implementation.
-fn divRec(mut &z: *[]Word, mut &u: *[]Word, mut v: []Word) {
+// around divRecursiveStep, the real implementation.
+fn divRecursive(mut &z: *[]Word, mut &u: *[]Word, mut v: []Word) {
 	// Recursion depth is (much) less than 2 log₂(len(v)).
 	// Allocate a slice of temporaries to be reused across recursion,
 	// plus one extra temporary not live across the recursion.
@@ -742,15 +735,15 @@ fn divRec(mut &z: *[]Word, mut &u: *[]Word, mut v: []Word) {
 	mut temps := make([][]Word, recDepth)
 
 	clearW(*z)
-	divRecStep(z, *u, v, 0, tmp, temps)
+	divRecursiveStep(z, *u, v, 0, tmp, temps)
 }
 
 // The actual implementation of recursive division.
 // It adds ⌊u/v⌋ to z and overwrites u with the remainder r.
 // z must be large enough to hold ⌊u/v⌋.
 // It uses temps[depth] (allocating if needed) as a temporary live across
 // the recursive call. It also uses tmp, but not live across the recursion.
-fn divRecStep(mut &z: *[]Word, mut u: []Word, mut v: []Word,
+fn divRecursiveStep(mut &z: *[]Word, mut u: []Word, mut v: []Word,
 	depth: int, mut tmp: []Word, mut temps: [][]Word) {
 	// u is a subsection of the original and may have leading zeros.
 	// TODO: The v = normW(v) is useless and should be removed.
@@ -783,7 +776,7 @@ fn divRecStep(mut &z: *[]Word, mut u: []Word, mut v: []Word,
 
 	// Allocate words for qhat below.
 	if temps[depth] == nil {
-		temps[depth] = getW(n) // TODO: Can be just W+1.
+		temps[depth] = getW(W + 1)
 	} else {
 		temps[depth] = makeW(temps[depth], W+1)
 	}
@@ -816,11 +809,11 @@ fn divRecStep(mut &z: *[]Word, mut u: []Word, mut v: []Word,
 		// Compute the 2-by-1 guess q̂, leaving r̂ in uu[s:W+n].
 		mut qhat := temps[depth]
 		clearW(qhat)
-		divRecStep(&qhat, uu[s:W+n], v[s:], depth+1, tmp, temps)
+		divRecursiveStep(&qhat, uu[s:W+n], v[s:], depth+1, tmp, temps)
 		qhat = normW(qhat)
 
 		// Extend to a 3-by-2 quotient and remainder.
-		// Because divRecStep overwrote the top part of uu with
+		// Because divRecursiveStep overwrote the top part of uu with
 		// the remainder r̂, the full uu already contains the equivalent
 		// of r̂·B + uₙ₋₂ from the “Refining Guesses” discussion.
 		// Subtracting q̂·vₙ₋₂ from it will compute the full-length remainder.
@@ -847,7 +840,7 @@ fn divRecStep(mut &z: *[]Word, mut u: []Word, mut v: []Word,
 			addAtW(uu[s:], v[s:], 0)
 		}
 		if cmpW(qhatv, normW(uu)) > 0 {
-			panic("math/big: divRecStep: impossible")
+			panic("math/big: divRecursiveStep: impossible")
 		}
 		c := subVV(uu[:len(qhatv)], uu[:len(qhatv)], qhatv)
 		if c > 0 {
@@ -864,7 +857,7 @@ fn divRecStep(mut &z: *[]Word, mut u: []Word, mut v: []Word,
 	s := W - 1
 	mut qhat := temps[depth]
 	clearW(qhat)
-	divRecStep(&qhat, normW(u[s:]), v[s:], depth+1, tmp, temps)
+	divRecursiveStep(&qhat, normW(u[s:]), v[s:], depth+1, tmp, temps)
 	qhat = normW(qhat)
 	mut qhatv := makeW(tmp, 3*n)
 	clearW(qhatv)
@@ -882,14 +875,14 @@ fn divRecStep(mut &z: *[]Word, mut u: []Word, mut v: []Word,
 		}
 	}
 	if cmpW(qhatv, normW(u)) > 0 {
-		panic("math/big: divRecStep: impossible")
+		panic("math/big: divRecursiveStep: impossible")
 	}
 	mut c := subVV(u[0:len(qhatv)], u[0:len(qhatv)], qhatv)
 	if c > 0 {
 		c = subVW(u[len(qhatv):], u[len(qhatv):], c)
 	}
 	if c > 0 {
-		panic("math/big: divRecStep: impossible")
+		panic("math/big: divRecursiveStep: impossible")
 	}
 
 	// Done!