math.big: replace division with Knuth, improve performance (#25242)

Author: tankf33der

vlib/math/big/array_ops.v

@@ -181,18 +181,17 @@ fn divide_digit_array(operand_a []u64, operand_b []u64, mut quotient []u64, mut
 	cmp_result := compare_digit_array(operand_a, operand_b)
 	// a == b => q, r = 1, 0
 	if cmp_result == 0 {
-		quotient << 1
-		for quotient.len > 1 {
-			quotient.delete_last()
-		}
+		quotient[0] = 1
 		remainder.clear()
 		return
 	}
 
 	// a < b => q, r = 0, a
 	if cmp_result < 0 {
 		quotient.clear()
-		remainder << operand_a
+		for i in 0 .. operand_a.len {
+			remainder[i] = operand_a[i]
+		}
 		return
 	}
 	if operand_b.len == 1 {
@@ -210,38 +209,36 @@ fn divide_array_by_digit(operand_a []u64, divisor u64, mut quotient []u64, mut r
 		// 1 digit for both dividend and divisor
 		dividend := operand_a[0]
 		q := dividend / divisor
-		if q != 0 {
-			quotient << q
-		}
+		quotient[0] = q
+
 		rem := dividend % divisor
-		if rem != 0 {
-			remainder << rem
-		}
+		remainder[0] = rem
+
+		shrink_tail_zeros(mut quotient)
+		shrink_tail_zeros(mut remainder)
 		return
 	}
 	// Dividend has more digits
 	mut rem := u64(0)
 	mut quo := u64(0)
-	mut qtemp := []u64{len: quotient.cap}
-	divisor64 := u64(divisor)
 
 	// Perform division step by step
 	for index := operand_a.len - 1; index >= 0; index-- {
 		hi := rem >> (64 - digit_bits)
 		lo := rem << digit_bits | operand_a[index]
-		quo, rem = bits.div_64(hi, lo, divisor64)
-		qtemp[index] = quo & max_digit
+		quo, rem = bits.div_64(hi, lo, divisor)
+		quotient[index] = quo & max_digit
 	}
+
 	// Remove leading zeros from quotient
-	shrink_tail_zeros(mut qtemp)
-	quotient << qtemp
-	remainder << rem
+	shrink_tail_zeros(mut quotient)
+	remainder[0] = rem
 	shrink_tail_zeros(mut remainder)
 }
 
 @[inline]
 fn divide_array_by_array(operand_a []u64, operand_b []u64, mut quotient []u64, mut remainder []u64) {
-	binary_divide_array_by_array(operand_a, operand_b, mut quotient, mut remainder)
+	knuth_divide_array_by_array(operand_a, operand_b, mut quotient, mut remainder)
 }
 
 // Shifts the contents of the original array by the given amount of bits to the left.

vlib/math/big/array_ops_test.v

@@ -136,8 +136,8 @@ fn test_compare_digit_array_02() {
 fn test_divide_digit_array_01() {
 	a := [u64(14)]
 	b := [u64(2)]
-	mut q := []u64{cap: 1}
-	mut r := []u64{cap: 1}
+	mut q := []u64{len: 1}
+	mut r := []u64{len: 1}
 
 	divide_digit_array(a, b, mut q, mut r)
 	assert q == [u64(7)]
@@ -147,8 +147,8 @@ fn test_divide_digit_array_01() {
 fn test_divide_digit_array_02() {
 	a := [u64(14)]
 	b := [u64(15)]
-	mut q := []u64{cap: 1}
-	mut r := []u64{cap: 1}
+	mut q := []u64{len: 1}
+	mut r := []u64{len: 1}
 
 	divide_digit_array(a, b, mut q, mut r)
 	assert q == []u64{len: 0}
@@ -158,8 +158,8 @@ fn test_divide_digit_array_02() {
 fn test_divide_digit_array_03() {
 	a := [u64(0), 4]
 	b := [u64(0), 1]
-	mut q := []u64{cap: a.len - b.len + 1}
-	mut r := []u64{cap: a.len}
+	mut q := []u64{len: a.len - b.len + 1}
+	mut r := []u64{len: a.len}
 
 	divide_digit_array(a, b, mut q, mut r)
 	assert q == [u64(4)]
@@ -169,8 +169,8 @@ fn test_divide_digit_array_03() {
 fn test_divide_digit_array_04() {
 	a := [u64(2), 4]
 	b := [u64(0), 1]
-	mut q := []u64{cap: a.len - b.len + 1}
-	mut r := []u64{cap: a.len}
+	mut q := []u64{len: a.len - b.len + 1}
+	mut r := []u64{len: a.len}
 
 	divide_digit_array(a, b, mut q, mut r)
 	assert q == [u64(4)]
@@ -180,8 +180,8 @@ fn test_divide_digit_array_04() {
 fn test_divide_digit_array_05() {
 	a := [u64(3)]
 	b := [u64(2)]
-	mut q := []u64{cap: a.len - b.len + 1}
-	mut r := []u64{cap: a.len}
+	mut q := []u64{len: a.len - b.len + 1}
+	mut r := []u64{len: a.len}
 
 	divide_digit_array(a, b, mut q, mut r)
 	assert q == [u64(1)]

vlib/math/big/division_array_ops.v

@@ -2,57 +2,6 @@ module big
 
 import math.bits
 
-// suppose operand_a bigger than operand_b and both not null.
-// Both quotient and remaider are allocated but of length 0
-@[direct_array_access]
-fn binary_divide_array_by_array(operand_a []u64, operand_b []u64, mut quotient []u64, mut remainder []u64) {
-	remainder << operand_a
-
-	len_diff := operand_a.len - operand_b.len
-	$if debug {
-		assert len_diff >= 0
-	}
-
-	// we must do in place shift and operations.
-	mut divisor := []u64{cap: operand_b.len}
-	for _ in 0 .. len_diff {
-		divisor << u64(0)
-	}
-	divisor << operand_b
-	for _ in 0 .. len_diff + 1 {
-		quotient << u64(0)
-	}
-
-	lead_zer_remainder := u32(bits.leading_zeros_64(remainder.last()) - (64 - digit_bits))
-	lead_zer_divisor := u32(bits.leading_zeros_64(divisor.last()) - (64 - digit_bits))
-	bit_offset := (u32(digit_bits) * u32(len_diff)) + (lead_zer_divisor - lead_zer_remainder)
-
-	// align
-	if lead_zer_remainder < lead_zer_divisor {
-		left_shift_in_place(mut divisor, lead_zer_divisor - lead_zer_remainder)
-	} else if lead_zer_remainder > lead_zer_divisor {
-		left_shift_in_place(mut remainder, lead_zer_remainder - lead_zer_divisor)
-	}
-
-	$if debug {
-		assert left_align_p(divisor[divisor.len - 1], remainder[remainder.len - 1])
-	}
-	for bit_idx := int(bit_offset); bit_idx >= 0; bit_idx-- {
-		if greater_equal_from_end(remainder, divisor) {
-			bit_set(mut quotient, bit_idx)
-			subtract_align_last_byte_in_place(mut remainder, divisor)
-		}
-		right_shift_in_place(mut divisor, 1)
-	}
-
-	// adjust
-	if lead_zer_remainder > lead_zer_divisor {
-		right_shift_in_place(mut remainder, lead_zer_remainder - lead_zer_divisor)
-	}
-	shrink_tail_zeros(mut remainder)
-	shrink_tail_zeros(mut quotient)
-}
-
 // help routines for cleaner code but inline for performance
 // quicker than BitField.set_bit
 @[direct_array_access; inline]
@@ -65,80 +14,112 @@ fn bit_set(mut a []u64, n int) {
 	a[byte_offset] |= mask
 }
 
-// a.len is greater or equal to b.len
-// returns true if a >= b (completed with zeroes)
-@[direct_array_access; inline]
-fn greater_equal_from_end(a []u64, b []u64) bool {
-	$if debug {
-		assert a.len >= b.len
-	}
-	offset := a.len - b.len
-	for index := a.len - 1; index >= offset; index-- {
-		if a[index] > b[index - offset] {
-			return true
-		} else if a[index] < b[index - offset] {
-			return false
-		}
-	}
-	return true
-}
+@[direct_array_access]
+fn knuth_divide_array_by_array(operand_a []u64, operand_b []u64, mut quotient []u64, mut remainder []u64) {
+	m := operand_a.len - operand_b.len
+	n := operand_b.len
+	mut u := []u64{len: operand_a.len + 1}
+	mut v := []u64{len: n}
+	leading_zeros := bits.leading_zeros_64(operand_b.last()) - (64 - digit_bits)
 
-// a := a - b supposed a >= b
-// attention the b operand is align with the a operand before the subtraction
-@[direct_array_access; inline]
-fn subtract_align_last_byte_in_place(mut a []u64, b []u64) {
-	mut carry := u64(0)
-	mut new_carry := u64(0)
-	offset := a.len - b.len
-	for index := a.len - b.len; index < a.len; index++ {
-		if a[index] < (b[index - offset] + carry) || (b[index - offset] == max_digit && carry > 0) {
-			new_carry = 1
-		} else {
-			new_carry = 0
+	if leading_zeros > 0 {
+		mut carry := u64(0)
+		amount := digit_bits - leading_zeros
+		for i in 0 .. operand_a.len {
+			temp := (operand_a[i] << leading_zeros) | carry
+			u[i] = temp & max_digit
+			carry = operand_a[i] >> amount
+		}
+		u[operand_a.len] = carry
+		carry = 0
+		for i in 0 .. operand_b.len {
+			temp := (operand_b[i] << leading_zeros) | carry
+			v[i] = temp & max_digit
+			carry = operand_b[i] >> amount
+		}
+	} else {
+		for i in 0 .. operand_a.len {
+			u[i] = operand_a[i]
+		}
+		for i in 0 .. operand_b.len {
+			v[i] = operand_b[i]
 		}
-		a[index] -= (b[index - offset] + carry)
-		a[index] = a[index] & max_digit
-		carry = new_carry
 	}
-	$if debug {
-		assert carry == 0
+
+	if remainder.len >= (n + 1) {
+		remainder.trim(n + 1)
+	} else {
+		remainder = []u64{len: n + 1}
 	}
-}
 
-// logical left shift
-// there is no overflow. We know that the last bits are zero
-// and that n <= `digit_bits`
-@[direct_array_access; inline]
-fn left_shift_in_place(mut a []u64, n u32) {
-	mut carry := u64(0)
-	mut prec_carry := u64(0)
-	mask := ((u64(1) << n) - 1) << (digit_bits - n)
-	for index in 0 .. a.len {
-		prec_carry = carry >> (digit_bits - n)
-		carry = a[index] & mask
-		a[index] <<= n
-		a[index] = a[index] & max_digit
-		a[index] |= prec_carry
+	v_n_1 := v[n - 1]
+	v_n_2 := v[n - 2]
+	for j := m; j >= 0; j-- {
+		u_j_n := u[j + n]
+		u_j_n_1 := u[j + n - 1]
+		u_j_n_2 := u[j + n - 2]
+
+		mut qhat, mut rhat := bits.div_64(u_j_n >> (64 - digit_bits), (u_j_n << digit_bits) | u_j_n_1,
+			v_n_1)
+		mut x1, mut x2 := bits.mul_64(qhat, v_n_2)
+		x2 = x2 & max_digit
+		x1 = (x1 << (64 - digit_bits)) | (x2 >> digit_bits)
+		for greater_than(x1, x2, rhat, u_j_n_2) {
+			qhat--
+			prev := rhat
+			rhat += v_n_1
+			if rhat < prev {
+				break
+			}
+			x1, x2 = bits.mul_64(qhat, v_n_2)
+			x2 = x2 & max_digit
+			x1 = (x1 << (64 - digit_bits)) | (x2 >> digit_bits)
+		}
+		mut carry := u64(0)
+		for i in 0 .. n {
+			hi, lo := bits.mul_add_64(v[i], qhat, carry)
+			remainder[i] = lo & max_digit
+			carry = (hi << (64 - digit_bits)) | (lo >> digit_bits)
+		}
+		remainder[n] = carry
+
+		mut borrow := u64(0)
+		for i in 0 .. n + 1 {
+			result := u[j + i] - remainder[i] - borrow
+			u[j + i] = result & max_digit
+			borrow = (result >> digit_bits) & 1
+		}
+		if borrow == 1 {
+			qhat--
+			carry = u64(0)
+			for i in 0 .. n {
+				sum := u[j + i] + v[i] + carry
+				u[j + i] = sum & max_digit
+				carry = sum >> digit_bits
+			}
+		}
+		quotient[j] = qhat
 	}
-}
 
-// logical right shift without control because these digits have already been
-// shift left before
-@[direct_array_access; inline]
-fn right_shift_in_place(mut a []u64, n u32) {
-	mut carry := u64(0)
-	mut prec_carry := u64(0)
-	mask := (u64(1) << n) - 1
-	for index := a.len - 1; index >= 0; index-- {
-		carry = a[index] & mask
-		a[index] >>= n
-		a[index] |= prec_carry << (digit_bits - n)
-		prec_carry = carry
+	remainder.delete_last()
+	if leading_zeros > 0 {
+		mut carry := u64(0)
+		max_leading_digit := (u64(1) << leading_zeros) - 1
+		for i := n - 1; i >= 0; i-- {
+			current_limb := u[i]
+			remainder[i] = (current_limb >> leading_zeros) | carry
+			carry = (current_limb & max_leading_digit) << (digit_bits - leading_zeros)
+		}
+	} else {
+		for i in 0 .. n {
+			remainder[i] = u[i]
+		}
 	}
+	shrink_tail_zeros(mut quotient)
+	shrink_tail_zeros(mut remainder)
 }
 
-// for assert
 @[inline]
-fn left_align_p(a u64, b u64) bool {
-	return bits.leading_zeros_64(a) == bits.leading_zeros_64(b)
+fn greater_than(x1 u64, x2 u64, y1 u64, y2 u64) bool {
+	return x1 > y1 || (x1 == y1 && x2 > y2)
 }