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

This commit is contained in:
Mike 2025-09-05 16:47:25 +03:00 committed by GitHub
parent 4dcf79e388
commit 1820a5584b
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
5 changed files with 135 additions and 238 deletions

View file

@ -181,10 +181,7 @@ 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
}
@ -192,7 +189,9 @@ fn divide_digit_array(operand_a []u64, operand_b []u64, mut quotient []u64, mut
// 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.

View file

@ -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)]

View file

@ -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
@[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)
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]
}
}
return true
}
// 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 remainder.len >= (n + 1) {
remainder.trim(n + 1)
} else {
remainder = []u64{len: n + 1}
}
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)
}
a[index] -= (b[index - offset] + carry)
a[index] = a[index] & max_digit
carry = new_carry
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
}
$if debug {
assert carry == 0
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)
}
// 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
}
}
// 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
}
}
// 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)
}

View file

@ -2,92 +2,11 @@ module big
import rand
fn test_left_shift_in_place() {
mut a := [u64(1), 1, 1, 1, 1]
left_shift_in_place(mut a, 1)
assert a == [u64(2), 2, 2, 2, 2]
left_shift_in_place(mut a, 7)
assert a == [u64(256), 256, 256, 256, 256]
mut b := [u64(0x08000000_00000001), 0x0c000000_00000000, 0x08000000_00000000, 0x07ffffff_ffffffff]
left_shift_in_place(mut b, 1)
assert b == [u64(2), 0x08000000_00000001, 1, 0x0fffffff_ffffffff]
mut c := [u64(0x00ffffff_ffffffff), 0x00f0f0f0_f0f0f0f0, 1, 0x03ffffff_ffffffff, 1]
left_shift_in_place(mut c, 2)
assert c == [u64(0x03ffffff_fffffffc), 0x03c3c3c3_c3c3c3c0, 4, 0x0fffffff_fffffffc, 4]
}
fn test_right_shift_in_place() {
mut a := [u64(2), 2, 2, 2, 2]
right_shift_in_place(mut a, 1)
assert a == [u64(1), 1, 1, 1, 1]
a = [u64(256), 256, 256, 256, 256]
right_shift_in_place(mut a, 7)
assert a == [u64(2), 2, 2, 2, 2]
a = [u64(0), 0, 1]
right_shift_in_place(mut a, 1)
assert a == [u64(0), 0x08000000_00000000, 0]
mut b := [u64(3), 0x08000000_00000001, 1, 0x0fffffff_ffffffff]
right_shift_in_place(mut b, 1)
assert b == [u64(0x08000000_00000001), 0x0c000000_00000000, 0x08000000_00000000,
0x07ffffff_ffffffff]
mut c := [u64(0x03ffffff), 0x03c3c3c3_c3c3c3c0, 7, 0xfffffffc, 4]
right_shift_in_place(mut c, 2)
assert c == [u64(0x00ffffff), 0x0cf0f0f0_f0f0f0f0, 1, 0x3fffffff, 1]
}
fn test_subtract_align_last_byte_in_place() {
mut a := [u64(2), 2, 2, 2, 2]
mut b := [u64(1), 1, 2, 1, 1]
subtract_align_last_byte_in_place(mut a, b)
assert a == [u64(1), 1, 0, 1, 1]
a = [u64(0), 0, 0, 0, 1]
b = [u64(0), 0, 1]
subtract_align_last_byte_in_place(mut a, b)
assert a == [u64(0), 0, 0, 0, 0]
a = [u64(0), 0, 0, 0, 1, 13]
b = [u64(1), 0, 1]
mut c := []u64{len: a.len}
mut d := [u64(0), 0, 0]
d << b // to have same length
subtract_digit_array(a, d, mut c)
subtract_align_last_byte_in_place(mut a, b)
assert a == [u64(0), 0, 0, u64(-1) & max_digit, 0, 12]
assert c == a
}
fn test_greater_equal_from_end() {
mut a := [u64(1), 2, 3, 4, 5, 6]
mut b := [u64(3), 4, 5, 6]
assert greater_equal_from_end(a, b) == true
a = [u64(1), 2, 3, 4, 5, 6]
b = [u64(1), 2, 3, 4, 5, 6]
assert greater_equal_from_end(a, b) == true
a = [u64(1), 2, 3, 4, 5, 6]
b = [u64(2), 2, 3, 4, 5, 6]
assert greater_equal_from_end(a, b) == false
a = [u64(0), 0, 0, 4, 5, 6]
b = [u64(4), 5, 6]
assert greater_equal_from_end(a, b) == true
a = [u64(0), 0, 0, 4, 5, 6]
b = [u64(4), 6, 6]
assert greater_equal_from_end(a, b) == false
a = [u64(0), 0, 0, 4, 5, 5]
b = [u64(4), 5, 6]
assert greater_equal_from_end(a, b) == false
}
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)]
@ -97,8 +16,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)]
@ -108,8 +27,8 @@ fn test_divide_digit_array_04() {
fn test_divide_digit_array_05() {
a := [u64(2), 4, 5]
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), 5]
@ -119,8 +38,8 @@ fn test_divide_digit_array_05() {
fn test_divide_digit_array_06() {
a := [u64(2), 4, 5, 3]
b := [u64(0), 0x8000]
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(0xa000_00000000), 0x6000_00000000]

View file

@ -393,8 +393,8 @@ pub fn (multiplicand Integer) * (multiplier Integer) Integer {
//
// DO NOT use this method if the divisor has any chance of being 0.
fn (dividend Integer) div_mod_internal(divisor Integer) (Integer, Integer) {
mut q := []u64{cap: int_max(1, dividend.digits.len - divisor.digits.len + 1)}
mut r := []u64{cap: dividend.digits.len}
mut q := []u64{len: int_max(1, dividend.digits.len - divisor.digits.len + 1)}
mut r := []u64{len: dividend.digits.len}
mut q_signum := 0
mut r_signum := 0