From a1b131c99aa9f3c0a96827d471c35bb3cead6eed Mon Sep 17 00:00:00 2001 From: kbkpbot Date: Thu, 14 Aug 2025 14:59:35 +0800 Subject: [PATCH] math.big: move from u32 to u60 digits (#25018) --- vlib/math/big/array_ops.v | 121 +++--- vlib/math/big/array_ops_test.v | 198 ++++----- .../bench_big_dividing_big_numbers_big_pow.v | 31 ++ vlib/math/big/big.v | 8 +- vlib/math/big/big_test.v | 2 +- vlib/math/big/division_array_ops.v | 59 +-- vlib/math/big/division_array_ops_test.v | 129 +++--- vlib/math/big/exponentiation.v | 8 +- vlib/math/big/integer.v | 403 ++++++------------ vlib/math/big/special_array_ops.v | 131 ++++-- vlib/math/big/special_array_ops_test.v | 140 +++--- 11 files changed, 593 insertions(+), 637 deletions(-) create mode 100644 vlib/math/big/bench/bench_big_dividing_big_numbers_big_pow.v diff --git a/vlib/math/big/array_ops.v b/vlib/math/big/array_ops.v index 6b1da19c6..90397cf8c 100644 --- a/vlib/math/big/array_ops.v +++ b/vlib/math/big/array_ops.v @@ -1,10 +1,12 @@ module big +import math.bits + // Compares the magnitude of the two unsigned integers represented the given // digit arrays. Returns -1 if a < b, 0 if a == b and +1 if a > b. Here // a is operand_a and b is operand_b (for brevity). @[direct_array_access] -fn compare_digit_array(operand_a []u32, operand_b []u32) int { +fn compare_digit_array(operand_a []u64, operand_b []u64) int { a_len := operand_a.len b_len := operand_b.len if a_len != b_len { @@ -26,7 +28,7 @@ fn compare_digit_array(operand_a []u32, operand_b []u32) int { // This function does not perform any allocation and assumes that the storage is // large enough. It may affect the last element, based on the presence of a carry @[direct_array_access] -fn add_digit_array(operand_a []u32, operand_b []u32, mut sum []u32) { +fn add_digit_array(operand_a []u64, operand_b []u64, mut sum []u64) { // Zero length cases if operand_a.len == 0 { for index in 0 .. operand_b.len { @@ -50,20 +52,20 @@ fn add_digit_array(operand_a []u32, operand_b []u32, mut sum []u32) { mut carry := u64(0) for index in 0 .. smaller_limit { partial := carry + a[index] + b[index] - sum[index] = u32(partial) - carry = u32(partial >> 32) + sum[index] = u64(partial) & max_digit + carry = u64(partial >> digit_bits) } for index in smaller_limit .. larger_limit { partial := carry + a[index] - sum[index] = u32(partial) - carry = u32(partial >> 32) + sum[index] = u64(partial) & max_digit + carry = u64(partial >> digit_bits) } if carry == 0 { sum.delete_last() } else { - sum[larger_limit] = u32(carry) + sum[larger_limit] = carry } } @@ -71,7 +73,7 @@ fn add_digit_array(operand_a []u32, operand_b []u32, mut sum []u32) { // It assumes operand_a contains the larger "integer" and that storage is // the same size as operand_a and is 0 @[direct_array_access] -fn subtract_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) { +fn subtract_digit_array(operand_a []u64, operand_b []u64, mut storage []u64) { // Zero length cases if operand_a.len == 0 { // nothing to subtract from @@ -86,23 +88,23 @@ fn subtract_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) { mut carry := false for index in 0 .. operand_b.len { - mut a_digit := u64(operand_a[index]) + mut a_digit := operand_a[index] b_digit := operand_b[index] + if carry { u64(1) } else { u64(0) } carry = a_digit < b_digit if carry { - a_digit += 0x100000000 + a_digit = a_digit | (u64(1) << digit_bits) } - storage[index] = u32(a_digit - b_digit) + storage[index] = a_digit - b_digit } for index in operand_b.len .. operand_a.len { - mut a_digit := u64(operand_a[index]) + mut a_digit := operand_a[index] b_digit := if carry { u64(1) } else { u64(0) } carry = a_digit < b_digit if carry { - a_digit += 0x100000000 + a_digit = a_digit | (u64(1) << digit_bits) } - storage[index] = u32(a_digit - b_digit) + storage[index] = a_digit - b_digit } shrink_tail_zeros(mut storage) @@ -113,7 +115,7 @@ const karatsuba_multiplication_limit = 70 const toom3_multiplication_limit = 360 @[inline] -fn multiply_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) { +fn multiply_digit_array(operand_a []u64, operand_b []u64, mut storage []u64) { max_len := if operand_a.len >= operand_b.len { operand_a.len } else { @@ -132,17 +134,18 @@ fn multiply_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) { // stored in storage. It assumes that storage has length equal to the sum of lengths // of a and b. Length refers to length of array, that is, digit count. @[direct_array_access] -fn simple_multiply_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) { +fn simple_multiply_digit_array(operand_a []u64, operand_b []u64, mut storage []u64) { for b_index in 0 .. operand_b.len { - mut carry := u64(0) + mut hi := u64(0) + mut lo := u64(0) for a_index in 0 .. operand_a.len { - partial_product := u64(storage[a_index + b_index]) + carry + - u64(operand_a[a_index]) * u64(operand_b[b_index]) - storage[a_index + b_index] = u32(partial_product) - carry = partial_product >> 32 + hi, lo = bits.mul_add_64(operand_a[a_index], operand_b[b_index], storage[a_index + + b_index] + hi) + storage[a_index + b_index] = lo & max_digit + hi = (hi << (64 - digit_bits)) | (lo >> digit_bits) } - if carry != 0 { - storage[b_index + operand_a.len] = u32(carry) + if hi != 0 { + storage[b_index + operand_a.len] = hi } } shrink_tail_zeros(mut storage) @@ -151,7 +154,7 @@ fn simple_multiply_digit_array(operand_a []u32, operand_b []u32, mut storage []u // Stores the product of the unsigned (non-negative) integer represented in a and the digit in value // in the storage array. It assumes storage is pre-initialised and populated with 0's @[direct_array_access] -fn multiply_array_by_digit(operand_a []u32, value u32, mut storage []u32) { +fn multiply_array_by_digit(operand_a []u64, value u64, mut storage []u64) { if value == 0 { storage.clear() return @@ -163,18 +166,16 @@ fn multiply_array_by_digit(operand_a []u32, value u32, mut storage []u32) { shrink_tail_zeros(mut storage) return } - mut carry := u32(0) + mut hi := u64(0) + mut lo := u64(0) for index in 0 .. operand_a.len { - product := u64(operand_a[index]) * value + carry - storage[index] = u32(product) - carry = u32(product >> 32) - } - if carry > 0 { - if storage.last() == 0 { - storage[operand_a.len] = carry - } else { - storage << carry - } + hi, lo = bits.mul_add_64(operand_a[index], value, hi) + storage[index] = lo & max_digit + hi = hi << (64 - digit_bits) + (lo >> digit_bits) + } + + if hi > 0 { + storage[operand_a.len] = hi } shrink_tail_zeros(mut storage) } @@ -184,7 +185,7 @@ fn multiply_array_by_digit(operand_a []u32, value u32, mut storage []u32) { // because it assumes that quotient and remainder are empty zero length arrays. They can be // made to have appropriate capacity though @[direct_array_access] -fn divide_digit_array(operand_a []u32, operand_b []u32, mut quotient []u32, mut remainder []u32) { +fn divide_digit_array(operand_a []u64, operand_b []u64, mut quotient []u64, mut remainder []u64) { cmp_result := compare_digit_array(operand_a, operand_b) // a == b => q, r = 1, 0 if cmp_result == 0 { @@ -212,7 +213,7 @@ fn divide_digit_array(operand_a []u32, operand_b []u32, mut quotient []u32, mut // Performs division on the non-negative dividend in a by the single digit divisor b. It assumes // quotient and remainder are empty zero length arrays without previous allocation @[direct_array_access] -fn divide_array_by_digit(operand_a []u32, divisor u32, mut quotient []u32, mut remainder []u32) { +fn divide_array_by_digit(operand_a []u64, divisor u64, mut quotient []u64, mut remainder []u64) { if operand_a.len == 1 { // 1 digit for both dividend and divisor dividend := operand_a[0] @@ -228,26 +229,28 @@ fn divide_array_by_digit(operand_a []u32, divisor u32, mut quotient []u32, mut r } // Dividend has more digits mut rem := u64(0) - mut qtemp := []u32{len: quotient.cap} + 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-- { - dividend := (rem << 32) + operand_a[index] - qtemp[index] = u32(dividend / divisor64) - rem = dividend % divisor64 + 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 } // Remove leading zeros from quotient shrink_tail_zeros(mut qtemp) quotient << qtemp - remainder << u32(rem) + remainder << rem shrink_tail_zeros(mut remainder) } const newton_division_limit = 1_000_000 @[inline] -fn divide_array_by_array(operand_a []u32, operand_b []u32, mut quotient []u32, mut remainder []u32) { +fn divide_array_by_array(operand_a []u64, operand_b []u64, mut quotient []u64, mut remainder []u64) { if operand_a.len >= newton_division_limit { newton_divide_array_by_array(operand_a, operand_b, mut quotient, mut remainder) } else { @@ -256,15 +259,15 @@ fn divide_array_by_array(operand_a []u32, operand_b []u32, mut quotient []u32, m } // Shifts the contents of the original array by the given amount of bits to the left. -// This function assumes that the amount is less than 32. The storage is expected to +// This function assumes that the amount is less than `digit_bits`. The storage is expected to // allocated with zeroes. @[direct_array_access] -fn shift_digits_left(original []u32, amount u32, mut storage []u32) { - mut leftover := u32(0) - offset := 32 - amount +fn shift_digits_left(original []u64, amount u32, mut storage []u64) { + mut leftover := u64(0) + offset := digit_bits - amount for index in 0 .. original.len { - value := leftover | (original[index] << amount) - leftover = (original[index] & (u32(-1) << offset)) >> offset + value := (leftover | (original[index] << amount)) & max_digit + leftover = (original[index] & (u64(-1) << offset)) >> offset storage[index] = value } if leftover != 0 { @@ -273,13 +276,13 @@ fn shift_digits_left(original []u32, amount u32, mut storage []u32) { } // Shifts the contents of the original array by the given amount of bits to the right. -// This function assumes that the amount is less than 32. The storage is expected to +// This function assumes that the amount is less than `digit_bits`. The storage is expected to // be allocated with zeroes. @[direct_array_access] -fn shift_digits_right(original []u32, amount u32, mut storage []u32) { - mut moveover := u32(0) - mask := (u32(1) << amount) - 1 - offset := 32 - amount +fn shift_digits_right(original []u64, amount u32, mut storage []u64) { + mut moveover := u64(0) + mask := (u64(1) << amount) - 1 + offset := digit_bits - amount for index := original.len - 1; index >= 0; index-- { value := (moveover << offset) | (original[index] >> amount) moveover = original[index] & mask @@ -289,7 +292,7 @@ fn shift_digits_right(original []u32, amount u32, mut storage []u32) { } @[direct_array_access] -fn bitwise_or_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) { +fn bitwise_or_digit_array(operand_a []u64, operand_b []u64, mut storage []u64) { lower, upper, bigger := if operand_a.len < operand_b.len { operand_a.len, operand_b.len, operand_b } else { @@ -305,7 +308,7 @@ fn bitwise_or_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) { } @[direct_array_access] -fn bitwise_and_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) { +fn bitwise_and_digit_array(operand_a []u64, operand_b []u64, mut storage []u64) { lower := imin(operand_a.len, operand_b.len) for index in 0 .. lower { storage[index] = operand_a[index] & operand_b[index] @@ -314,7 +317,7 @@ fn bitwise_and_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) } @[direct_array_access] -fn bitwise_xor_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) { +fn bitwise_xor_digit_array(operand_a []u64, operand_b []u64, mut storage []u64) { lower, upper, bigger := if operand_a.len < operand_b.len { operand_a.len, operand_b.len, operand_b } else { @@ -330,9 +333,9 @@ fn bitwise_xor_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) } @[direct_array_access] -fn bitwise_not_digit_array(original []u32, mut storage []u32) { +fn bitwise_not_digit_array(original []u64, mut storage []u64) { for index in 0 .. original.len { - storage[index] = ~original[index] + storage[index] = (~original[index]) & max_digit } shrink_tail_zeros(mut storage) } diff --git a/vlib/math/big/array_ops_test.v b/vlib/math/big/array_ops_test.v index 5c7448b74..88441a610 100644 --- a/vlib/math/big/array_ops_test.v +++ b/vlib/math/big/array_ops_test.v @@ -2,120 +2,120 @@ module big fn test_add_digit_array_01() { - a := [u32(1), 1, 1] - b := [u32(1), 1, 1] - mut c := []u32{len: 4} + a := [u64(1), 1, 1] + b := [u64(1), 1, 1] + mut c := []u64{len: 4} add_digit_array(a, b, mut c) - assert c == [u32(2), 2, 2] + assert c == [u64(2), 2, 2] } fn test_add_digit_array_02() { - a := [u32(1), u32(1) << 31, 1] - b := [u32(1), u32(1) << 31, 1] - mut c := []u32{len: 4} + a := [u64(1), u64(1) << (digit_bits - 1), 1] + b := [u64(1), u64(1) << (digit_bits - 1), 1] + mut c := []u64{len: 4} add_digit_array(a, b, mut c) - assert c == [u32(2), 0, 3] + assert c == [u64(2), 0, 3] } fn test_add_digit_array_03() { - a := [u32(1), (u32(1) << 31) + u32(34), 1] - b := [u32(242), u32(1) << 31, 1] - mut c := []u32{len: 4} + a := [u64(1), (u64(1) << (digit_bits - 1)) + u64(34), 1] + b := [u64(242), u64(1) << (digit_bits - 1), 1] + mut c := []u64{len: 4} add_digit_array(a, b, mut c) - assert c == [u32(243), 34, 3] + assert c == [u64(243), 34, 3] } fn test_add_digit_array_04() { - a := [u32(0)] - b := [u32(1), 3, 4] - mut c := []u32{len: 4} + a := [u64(0)] + b := [u64(1), 3, 4] + mut c := []u64{len: 4} add_digit_array(a, b, mut c) - assert c == [u32(1), 3, 4] + assert c == [u64(1), 3, 4] } fn test_add_digit_array_05() { - a := [u32(1), 3, 4] - b := [u32(0)] - mut c := []u32{len: 4} + a := [u64(1), 3, 4] + b := [u64(0)] + mut c := []u64{len: 4} add_digit_array(a, b, mut c) - assert c == [u32(1), 3, 4] + assert c == [u64(1), 3, 4] } fn test_add_digit_array_06() { - a := [u32(46), 13, 462, 13] - b := [u32(1), 3, 4] - mut c := []u32{len: 5} + a := [u64(46), 13, 462, 13] + b := [u64(1), 3, 4] + mut c := []u64{len: 5} add_digit_array(a, b, mut c) - assert c == [u32(47), 16, 466, 13] + assert c == [u64(47), 16, 466, 13] } fn test_subtract_digit_array_01() { - a := [u32(2), 2, 2, 2, 2] - b := [u32(1), 1, 2, 1, 1] - mut c := []u32{len: a.len} + a := [u64(2), 2, 2, 2, 2] + b := [u64(1), 1, 2, 1, 1] + mut c := []u64{len: a.len} subtract_digit_array(a, b, mut c) - assert c == [u32(1), 1, 0, 1, 1] + assert c == [u64(1), 1, 0, 1, 1] } fn test_subtract_digit_array_02() { - a := [u32(0), 0, 0, 0, 1] - b := [u32(0), 0, 1] - mut c := []u32{len: a.len} + a := [u64(0), 0, 0, 0, 1] + b := [u64(0), 0, 1] + mut c := []u64{len: a.len} subtract_digit_array(a, b, mut c) - assert c == [u32(0), 0, u32(-1), u32(-1)] + assert c == [u64(0), 0, u64(-1) & max_digit, u64(-1) & max_digit] } fn test_subtract_digit_array_03() { - a := [u32(0), 0, 0, 0, 1, 13] - b := [u32(0), 0, 1] - mut c := []u32{len: a.len} + a := [u64(0), 0, 0, 0, 1, 13] + b := [u64(0), 0, 1] + mut c := []u64{len: a.len} subtract_digit_array(a, b, mut c) - assert c == [u32(0), 0, u32(-1), u32(-1), 0, 13] + assert c == [u64(0), 0, u64(-1) & max_digit, u64(-1) & max_digit, 0, 13] } fn test_subtract_digit_array_04() { - a := [u32(0x2), 0x4, 0x5, 0x3] - b := [u32(0x0), 0x0, 0x5, 0x3] - mut c := []u32{len: a.len} + a := [u64(0x2), 0x4, 0x5, 0x3] + b := [u64(0x0), 0x0, 0x5, 0x3] + mut c := []u64{len: a.len} subtract_digit_array(a, b, mut c) - assert c == [u32(0x2), 0x4] + assert c == [u64(0x2), 0x4] } fn test_multiply_digit_array_01() { - a := [u32(0), 0, 0, 1] - b := [u32(0), 0, 1] - mut c := []u32{len: a.len + b.len} + a := [u64(0), 0, 0, 1] + b := [u64(0), 0, 1] + mut c := []u64{len: a.len + b.len} multiply_digit_array(a, b, mut c) - assert c == [u32(0), 0, 0, 0, 0, 1] + assert c == [u64(0), 0, 0, 0, 0, 1] } fn test_multiply_digit_array_02() { - a := []u32{len: 0} - b := [u32(0), 0, 1] - mut c := []u32{len: a.len + b.len} + a := []u64{len: 0} + b := [u64(0), 0, 1] + mut c := []u64{len: a.len + b.len} multiply_digit_array(a, b, mut c) assert c == [] - c = []u32{len: a.len + b.len} + c = []u64{len: a.len + b.len} multiply_digit_array(b, a, mut c) assert c == [] } fn test_compare_digit_array_01() { - a := [u32(0), 0, 2] - b := [u32(0), 0, 4] + a := [u64(0), 0, 2] + b := [u64(0), 0, 4] assert compare_digit_array(a, b) < 0 assert compare_digit_array(b, a) > 0 @@ -124,8 +124,8 @@ fn test_compare_digit_array_01() { } fn test_compare_digit_array_02() { - a := [u32(0), 0, 2324, 0, 124] - b := [u32(0), 0, 4, 0, 0, 1] + a := [u64(0), 0, 2324, 0, 124] + b := [u64(0), 0, 4, 0, 0, 1] assert compare_digit_array(a, b) < 0 assert compare_digit_array(b, a) > 0 @@ -134,64 +134,64 @@ fn test_compare_digit_array_02() { } fn test_divide_digit_array_01() { - a := [u32(14)] - b := [u32(2)] - mut q := []u32{cap: 1} - mut r := []u32{cap: 1} + a := [u64(14)] + b := [u64(2)] + mut q := []u64{cap: 1} + mut r := []u64{cap: 1} divide_digit_array(a, b, mut q, mut r) - assert q == [u32(7)] - assert r == []u32{len: 0} + assert q == [u64(7)] + assert r == []u64{len: 0} } fn test_divide_digit_array_02() { - a := [u32(14)] - b := [u32(15)] - mut q := []u32{cap: 1} - mut r := []u32{cap: 1} + a := [u64(14)] + b := [u64(15)] + mut q := []u64{cap: 1} + mut r := []u64{cap: 1} divide_digit_array(a, b, mut q, mut r) - assert q == []u32{len: 0} + assert q == []u64{len: 0} assert r == a } fn test_divide_digit_array_03() { - a := [u32(0), 4] - b := [u32(0), 1] - mut q := []u32{cap: a.len - b.len + 1} - mut r := []u32{cap: a.len} + a := [u64(0), 4] + b := [u64(0), 1] + mut q := []u64{cap: a.len - b.len + 1} + mut r := []u64{cap: a.len} divide_digit_array(a, b, mut q, mut r) - assert q == [u32(4)] - assert r == []u32{len: 0} + assert q == [u64(4)] + assert r == []u64{len: 0} } fn test_divide_digit_array_04() { - a := [u32(2), 4] - b := [u32(0), 1] - mut q := []u32{cap: a.len - b.len + 1} - mut r := []u32{cap: a.len} + a := [u64(2), 4] + b := [u64(0), 1] + mut q := []u64{cap: a.len - b.len + 1} + mut r := []u64{cap: a.len} divide_digit_array(a, b, mut q, mut r) - assert q == [u32(4)] - assert r == [u32(2)] + assert q == [u64(4)] + assert r == [u64(2)] } fn test_divide_digit_array_05() { - a := [u32(3)] - b := [u32(2)] - mut q := []u32{cap: a.len - b.len + 1} - mut r := []u32{cap: a.len} + a := [u64(3)] + b := [u64(2)] + mut q := []u64{cap: a.len - b.len + 1} + mut r := []u64{cap: a.len} divide_digit_array(a, b, mut q, mut r) - assert q == [u32(1)] - assert r == [u32(1)] + assert q == [u64(1)] + assert r == [u64(1)] } fn test_left_and_right_shift() { - a := [u32(1), 1, 1] - mut r := [u32(2), 2, 2] - mut b := []u32{len: 3, init: 0} + a := [u64(1), 1, 1] + mut r := [u64(2), 2, 2] + mut b := []u64{len: 3, init: 0} shift_digits_left(a, 1, mut b) assert r == b shift_digits_right(r, 1, mut r) @@ -199,32 +199,34 @@ fn test_left_and_right_shift() { shift_digits_left(r, 1, mut r) assert r == b - mut c := [u32(0xffff_ffff)] + mut c := [u64(0x0fff_ffff_ffff_ffff)] shift_digits_left(c, 16, mut c) - assert c == [u32(0xffff_0000), 0xffff] + assert c == [u64(0x0fff_ffff_ffff_0000), u64(0xffff)] shift_digits_right(c, 8, mut c) - assert c == [u32(0xffff_ff00), 0xff] + assert c == [u64(0x0fff_ffff_ffff_ff00), 0xff] shift_digits_right(c, 16, mut c) - assert c == [u32(0x00ff_ffff)] + assert c == [u64(0x000f_ffff_ffff_ffff)] shift_digits_right(c, 16, mut c) - assert c == [u32(0xff)] + assert c == [u64(0x000f_ffff_ffff)] shift_digits_right(c, 16, mut c) - assert c == []u32{len: 0} + assert c == [u64(0x000f_ffff)] + shift_digits_right(c, 40, mut c) + assert c == []u64{len: 0} } fn test_or_digit_array() { - a := [u32(10), 10, 10] - b := [u32(5), 5, 5] - mut c := []u32{len: 3, init: 0} + a := [u64(10), 10, 10] + b := [u64(5), 5, 5] + mut c := []u64{len: 3, init: 0} bitwise_or_digit_array(a, b, mut c) - assert c == [u32(15), 15, 15] + assert c == [u64(15), 15, 15] bitwise_or_digit_array(a, a, mut c) assert c == a - x := [u32(10), 10, 10, 42, 42] - y := [u32(2), 2, 5, 2] - mut d := []u32{len: 5, init: 0} + x := [u64(10), 10, 10, 42, 42] + y := [u64(2), 2, 5, 2] + mut d := []u64{len: 5, init: 0} bitwise_or_digit_array(y, x, mut d) - assert d == [u32(10), 10, 15, 42, 42] + assert d == [u64(10), 10, 15, 42, 42] } diff --git a/vlib/math/big/bench/bench_big_dividing_big_numbers_big_pow.v b/vlib/math/big/bench/bench_big_dividing_big_numbers_big_pow.v new file mode 100644 index 000000000..aeb1ba2ac --- /dev/null +++ b/vlib/math/big/bench/bench_big_dividing_big_numbers_big_pow.v @@ -0,0 +1,31 @@ +module main + +import math.big +import benchmark + +const n = 2500_000 + +fn dividing_big_numbers_big_pow() { + d_12 := big.integer_from_int(12) + d_4 := big.integer_from_int(4) + d_10 := big.integer_from_int(10) + + mut clock := benchmark.start() + base := d_10.pow(u32(n - 1)) + clock.measure('math.big> n: ${n} | base.bit_len: ${base.bit_len()}') + a := d_12 * base + b := d_4 * base + clock.measure('math.big> a.bit_len: ${a.bit_len()} | b.bit_len: ${b.bit_len()}') + c := a / b // c should be 3 + clock.measure('math.big> c: ${c} | c.bit_len(): ${c.bit_len()}') + if c.str() != '3' { + println('math.big:error! c.str() == ${c.str()}') + } + if c.bit_len() != 2 { + println('math.big:error! c.bit_len() == ${c.bit_len()}') + } +} + +fn main() { + dividing_big_numbers_big_pow() +} diff --git a/vlib/math/big/big.v b/vlib/math/big/big.v index ddcd4cc42..2c0cb7ab4 100644 --- a/vlib/math/big/big.v +++ b/vlib/math/big/big.v @@ -1,22 +1,22 @@ module big pub const zero_int = Integer{ - digits: []u32{len: 0} + digits: []u64{len: 0} signum: 0 is_const: true } pub const one_int = Integer{ - digits: [u32(1)] + digits: [u64(1)] signum: 1 is_const: true } pub const two_int = Integer{ - digits: [u32(2)] + digits: [u64(2)] signum: 1 is_const: true } pub const three_int = Integer{ - digits: [u32(3)] + digits: [u64(3)] signum: 1 is_const: true } diff --git a/vlib/math/big/big_test.v b/vlib/math/big/big_test.v index e0abf7f17..2097b6a4e 100644 --- a/vlib/math/big/big_test.v +++ b/vlib/math/big/big_test.v @@ -913,7 +913,7 @@ fn test_isqrt() { } fn test_bitwise_ops() { - a := big.integer_from_int(1).left_shift(512) + a := big.integer_from_int(1).left_shift(big.digit_bits * 10) b := a - big.one_int assert a.bitwise_and(b) == big.zero_int assert b.bitwise_xor(b) == big.zero_int diff --git a/vlib/math/big/division_array_ops.v b/vlib/math/big/division_array_ops.v index e7f277498..81802f9ed 100644 --- a/vlib/math/big/division_array_ops.v +++ b/vlib/math/big/division_array_ops.v @@ -5,7 +5,7 @@ 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 []u32, operand_b []u32, mut quotient []u32, mut remainder []u32) { +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 @@ -14,18 +14,18 @@ fn binary_divide_array_by_array(operand_a []u32, operand_b []u32, mut quotient [ } // we must do in place shift and operations. - mut divisor := []u32{cap: operand_b.len} + mut divisor := []u64{cap: operand_b.len} for _ in 0 .. len_diff { - divisor << u32(0) + divisor << u64(0) } divisor << operand_b for _ in 0 .. len_diff + 1 { - quotient << u32(0) + quotient << u64(0) } - lead_zer_remainder := u32(bits.leading_zeros_32(remainder.last())) - lead_zer_divisor := u32(bits.leading_zeros_32(divisor.last())) - bit_offset := (u32(32) * u32(len_diff)) + (lead_zer_divisor - lead_zer_remainder) + 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 { @@ -47,7 +47,6 @@ fn binary_divide_array_by_array(operand_a []u32, operand_b []u32, mut quotient [ // adjust if lead_zer_remainder > lead_zer_divisor { - // right_shift_in_place(mut quotient, lead_zer_remainder - lead_zer_divisor) right_shift_in_place(mut remainder, lead_zer_remainder - lead_zer_divisor) } shrink_tail_zeros(mut remainder) @@ -57,9 +56,9 @@ fn binary_divide_array_by_array(operand_a []u32, operand_b []u32, mut quotient [ // help routines for cleaner code but inline for performance // quicker than BitField.set_bit @[direct_array_access; inline] -fn bit_set(mut a []u32, n int) { - byte_offset := n >> 5 - mask := u32(1) << u32(n % 32) +fn bit_set(mut a []u64, n int) { + byte_offset := n / digit_bits + mask := u64(1) << u32(n % digit_bits) $if debug { assert a.len >= byte_offset } @@ -69,7 +68,7 @@ fn bit_set(mut a []u32, n int) { // 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 []u32, b []u32) bool { +fn greater_equal_from_end(a []u64, b []u64) bool { $if debug { assert a.len >= b.len } @@ -87,17 +86,18 @@ fn greater_equal_from_end(a []u32, b []u32) bool { // 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 []u32, b []u32) { - mut carry := u32(0) - mut new_carry := u32(0) +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_u32 && carry > 0) { + if a[index] < (b[index - offset] + carry) || (b[index - offset] == max_digit && carry > 0) { new_carry = 1 } else { new_carry = 0 } a[index] -= (b[index - offset] + carry) + a[index] = a[index] & max_digit carry = new_carry } $if debug { @@ -107,16 +107,17 @@ fn subtract_align_last_byte_in_place(mut a []u32, b []u32) { // logical left shift // there is no overflow. We know that the last bits are zero -// and that n <= 32 +// and that n <= `digit_bits` @[direct_array_access; inline] -fn left_shift_in_place(mut a []u32, n u32) { - mut carry := u32(0) - mut prec_carry := u32(0) - mask := ((u32(1) << n) - 1) << (32 - n) +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 >> (32 - n) + prec_carry = carry >> (digit_bits - n) carry = a[index] & mask a[index] <<= n + a[index] = a[index] & max_digit a[index] |= prec_carry } } @@ -124,20 +125,20 @@ fn left_shift_in_place(mut a []u32, n u32) { // 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 []u32, n u32) { - mut carry := u32(0) - mut prec_carry := u32(0) - mask := u32((1 << n) - 1) +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 << (32 - n) + a[index] |= prec_carry << (digit_bits - n) prec_carry = carry } } // for assert @[inline] -fn left_align_p(a u32, b u32) bool { - return bits.leading_zeros_32(a) == bits.leading_zeros_32(b) +fn left_align_p(a u64, b u64) bool { + return bits.leading_zeros_64(a) == bits.leading_zeros_64(b) } diff --git a/vlib/math/big/division_array_ops_test.v b/vlib/math/big/division_array_ops_test.v index de8e0bb10..6f81c9427 100644 --- a/vlib/math/big/division_array_ops_test.v +++ b/vlib/math/big/division_array_ops_test.v @@ -3,127 +3,128 @@ module big import rand fn test_left_shift_in_place() { - mut a := [u32(1), 1, 1, 1, 1] + mut a := [u64(1), 1, 1, 1, 1] left_shift_in_place(mut a, 1) - assert a == [u32(2), 2, 2, 2, 2] + assert a == [u64(2), 2, 2, 2, 2] left_shift_in_place(mut a, 7) - assert a == [u32(256), 256, 256, 256, 256] - mut b := [u32(0x80000001), 0xc0000000, 0x80000000, 0x7fffffff] + 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 == [u32(2), 0x80000001, 1, 0xffffffff] - mut c := [u32(0x00ffffff), 0xf0f0f0f0, 1, 0x3fffffff, 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 == [u32(0x3fffffc), 0xc3c3c3c0, 7, 0xfffffffc, 4] + assert c == [u64(0x03ffffff_fffffffc), 0x03c3c3c3_c3c3c3c0, 4, 0x0fffffff_fffffffc, 4] } fn test_right_shift_in_place() { - mut a := [u32(2), 2, 2, 2, 2] + mut a := [u64(2), 2, 2, 2, 2] right_shift_in_place(mut a, 1) - assert a == [u32(1), 1, 1, 1, 1] - a = [u32(256), 256, 256, 256, 256] + assert a == [u64(1), 1, 1, 1, 1] + a = [u64(256), 256, 256, 256, 256] right_shift_in_place(mut a, 7) - assert a == [u32(2), 2, 2, 2, 2] - a = [u32(0), 0, 1] + assert a == [u64(2), 2, 2, 2, 2] + a = [u64(0), 0, 1] right_shift_in_place(mut a, 1) - assert a == [u32(0), 0x80000000, 0] - mut b := [u32(3), 0x80000001, 1, 0xffffffff] + 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 == [u32(0x80000001), 0xc0000000, 0x80000000, 0x7fffffff] - mut c := [u32(0x03ffffff), 0xc3c3c3c0, 7, 0xfffffffc, 4] + 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 == [u32(0x00ffffff), 0xf0f0f0f0, 1, 0x3fffffff, 1] + assert c == [u64(0x00ffffff), 0x0cf0f0f0_f0f0f0f0, 1, 0x3fffffff, 1] } fn test_subtract_align_last_byte_in_place() { - mut a := [u32(2), 2, 2, 2, 2] - mut b := [u32(1), 1, 2, 1, 1] + 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 == [u32(1), 1, 0, 1, 1] + assert a == [u64(1), 1, 0, 1, 1] - a = [u32(0), 0, 0, 0, 1] - b = [u32(0), 0, 1] + a = [u64(0), 0, 0, 0, 1] + b = [u64(0), 0, 1] subtract_align_last_byte_in_place(mut a, b) - assert a == [u32(0), 0, 0, 0, 0] + assert a == [u64(0), 0, 0, 0, 0] - a = [u32(0), 0, 0, 0, 1, 13] - b = [u32(1), 0, 1] - mut c := []u32{len: a.len} - mut d := [u32(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 == [u32(0), 0, 0, u32(-1), 0, 12] + assert a == [u64(0), 0, 0, u64(-1) & max_digit, 0, 12] assert c == a } fn test_greater_equal_from_end() { - mut a := [u32(1), 2, 3, 4, 5, 6] - mut b := [u32(3), 4, 5, 6] + 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 = [u32(1), 2, 3, 4, 5, 6] - b = [u32(1), 2, 3, 4, 5, 6] + 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 = [u32(1), 2, 3, 4, 5, 6] - b = [u32(2), 2, 3, 4, 5, 6] + 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 = [u32(0), 0, 0, 4, 5, 6] - b = [u32(4), 5, 6] + a = [u64(0), 0, 0, 4, 5, 6] + b = [u64(4), 5, 6] assert greater_equal_from_end(a, b) == true - a = [u32(0), 0, 0, 4, 5, 6] - b = [u32(4), 6, 6] + a = [u64(0), 0, 0, 4, 5, 6] + b = [u64(4), 6, 6] assert greater_equal_from_end(a, b) == false - a = [u32(0), 0, 0, 4, 5, 5] - b = [u32(4), 5, 6] + 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 := [u32(0), 4] - b := [u32(0), 1] - mut q := []u32{cap: a.len - b.len + 1} - mut r := []u32{cap: a.len} + a := [u64(0), 4] + b := [u64(0), 1] + mut q := []u64{cap: a.len - b.len + 1} + mut r := []u64{cap: a.len} divide_digit_array(a, b, mut q, mut r) - assert q == [u32(4)] - assert r == []u32{len: 0} + assert q == [u64(4)] + assert r == []u64{len: 0} } fn test_divide_digit_array_04() { - a := [u32(2), 4] - b := [u32(0), 1] - mut q := []u32{cap: a.len - b.len + 1} - mut r := []u32{cap: a.len} + a := [u64(2), 4] + b := [u64(0), 1] + mut q := []u64{cap: a.len - b.len + 1} + mut r := []u64{cap: a.len} divide_digit_array(a, b, mut q, mut r) - assert q == [u32(4)] - assert r == [u32(2)] + assert q == [u64(4)] + assert r == [u64(2)] } fn test_divide_digit_array_05() { - a := [u32(2), 4, 5] - b := [u32(0), 1] - mut q := []u32{cap: a.len - b.len + 1} - mut r := []u32{cap: a.len} + a := [u64(2), 4, 5] + b := [u64(0), 1] + mut q := []u64{cap: a.len - b.len + 1} + mut r := []u64{cap: a.len} divide_digit_array(a, b, mut q, mut r) - assert q == [u32(4), 5] - assert r == [u32(2)] + assert q == [u64(4), 5] + assert r == [u64(2)] } fn test_divide_digit_array_06() { - a := [u32(2), 4, 5, 3] - b := [u32(0), 0x8000] - mut q := []u32{cap: a.len - b.len + 1} - mut r := []u32{cap: a.len} + 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} divide_digit_array(a, b, mut q, mut r) - assert q == [u32(0xa0000), 0x60000] - assert r == [u32(2), 4] + assert q == [u64(0xa000_00000000), 0x6000_00000000] + assert r == [u64(2), 4] } fn test_many_divisions() { diff --git a/vlib/math/big/exponentiation.v b/vlib/math/big/exponentiation.v index 6142ceb2f..33f79f744 100644 --- a/vlib/math/big/exponentiation.v +++ b/vlib/math/big/exponentiation.v @@ -64,12 +64,12 @@ fn (a Integer) mont_odd(x Integer, m Integer) Integer { table[i] = table[i - 1].mont_mul(d, ctx) } } - mut r := if m.digits.last() & (1 << (32 - 1)) != 0 { - mut rdigits := []u32{len: m.digits.len} + mut r := if m.digits.last() & (u64(1) << (digit_bits - 1)) != 0 { + mut rdigits := []u64{len: m.digits.len} - rdigits[0] = -m.digits[0] + rdigits[0] = (-m.digits[0]) & max_digit for i := 1; i < m.digits.len; i++ { - rdigits[i] = ~m.digits[i] + rdigits[i] = (~m.digits[i]) & max_digit } Integer{ diff --git a/vlib/math/big/integer.v b/vlib/math/big/integer.v index d2bcb0287..15736297b 100644 --- a/vlib/math/big/integer.v +++ b/vlib/math/big/integer.v @@ -1,32 +1,23 @@ module big -import math import math.bits -import strings -import strconv const digit_array = '0123456789abcdefghijklmnopqrstuvwxyz'.bytes() -// vfmt off -const radix_options = { - 2: 31, 3: 20, 4: 15, 5: 13, 6: 12, 7: 11, 8: 10, 9: 10, 10: 9, - 11: 9, 12: 8, 13: 8, 14: 8, 15: 8, 16: 7, 17: 7, 18: 7, 19: 7, - 20: 7, 21: 7, 22: 7, 23: 7, 24: 6, 25: 6, 26: 6, 27: 6, 28: 6, - 29: 6, 30: 6, 31: 6, 32: 6, 33: 6, 34: 6, 35: 6, 36: 6 -} -// vfmt on - +pub const digit_bits = 60 // 60bits +const max_digit = (u64(1) << digit_bits) - u64(1) // big.Integer // ----------- // It has the following properties: -// 1. Every "digit" is an integer in the range [0, 2^32). +// 1. Every "digit" is an integer in the range [0, 2^digit_bits-1). // 2. The signum can be one of three values: -1, 0, +1 for // negative, zero, and positive values, respectively. // 3. There should be no leading zeros in the digit array. // 4. The digits are stored in little endian format, that is, // the digits with a lower positional value (towards the right // when represented as a string) have a lower index, and vice versa. +// 5. zero's signum is zero, digits.len = 0 pub struct Integer { - digits []u32 + digits []u64 // in one u64, use only `digit_bits` store a digit pub: signum int is_const bool @@ -61,7 +52,7 @@ pub fn integer_from_int(value int) Integer { return zero_int } return Integer{ - digits: [u32(iabs(value))] + digits: [u64(iabs(value))] signum: int_signum(value) } } @@ -72,7 +63,7 @@ pub fn integer_from_u32(value u32) Integer { return zero_int } return Integer{ - digits: [value] + digits: [u64(value)] signum: 1 } } @@ -84,10 +75,14 @@ pub fn integer_from_i64(value i64) Integer { } signum_value := if value < 0 { -1 } else { 1 } - abs_value := u64(value * signum_value) + abs_value := if value == i64(-9223372036854775808) { + u64(0x8000000000000000) + } else { + u64(value * signum_value) + } - lower := u32(abs_value) - upper := u32(abs_value >> 32) + lower := u64(abs_value & max_digit) + upper := u64(abs_value >> digit_bits) if upper == 0 { return Integer{ @@ -108,8 +103,8 @@ pub fn integer_from_u64(value u64) Integer { return zero_int } - lower := u32(value & 0x00000000ffffffff) - upper := u32((value & 0xffffffff00000000) >> 32) + lower := u64(value & max_digit) + upper := u64(value >> digit_bits) if upper == 0 { return Integer{ @@ -138,7 +133,7 @@ pub: pub fn integer_from_bytes(oinput []u8, config IntegerConfig) Integer { // Thank you to Miccah (@mcastorina) for this implementation and relevant unit tests. if oinput.len == 0 { - return integer_from_int(0) + return zero_int } // Ignore leading 0 bytes: mut first_non_zero_index := -1 @@ -149,23 +144,40 @@ pub fn integer_from_bytes(oinput []u8, config IntegerConfig) Integer { } } if first_non_zero_index == -1 { - return integer_from_int(0) + return zero_int } input := oinput[first_non_zero_index..] - // pad input - mut padded_input := []u8{len: int_max(0, ((input.len + 3) & ~0x3) - input.len), cap: ( - input.len + 3) & ~0x3} - padded_input << input - mut digits := []u32{len: padded_input.len / 4} - // combine every 4 bytes into a u32 and insert into n.digits - for i := 0; i < padded_input.len; i += 4 { - x3 := u32(padded_input[i]) - x2 := u32(padded_input[i + 1]) - x1 := u32(padded_input[i + 2]) - x0 := u32(padded_input[i + 3]) - val := (x3 << 24) | (x2 << 16) | (x1 << 8) | x0 - digits[(padded_input.len - i) / 4 - 1] = val + + mut carry_bits := 0 + mut carry_value := u64(0) + mut digits := []u64{} + + for i := input.len - 1; i >= 0; i-- { + byte_value := input[i] + for shift in 0 .. 8 { + bit := (byte_value >> u8(shift)) & 1 + carry_value = (carry_value >> 1) | (u64(bit) << (digit_bits - 1)) + carry_bits++ + if carry_bits == digit_bits { + digits << carry_value + carry_value = 0 + carry_bits = 0 + } + } + } + + if carry_bits > 0 { + remaining_shift := digit_bits - carry_bits + digits << (carry_value >> remaining_shift) } + + // Remove trailing zeros + shrink_tail_zeros(mut digits) + + if digits.len == 0 { + return zero_int + } + return Integer{ digits: digits signum: config.signum @@ -184,18 +196,8 @@ pub fn integer_from_radix(all_characters string, radix u32) !Integer { return error('math.big: Radix must be between 2 and 36 (inclusive)') } characters := all_characters.to_lower() - validate_string(characters, radix) or { return err } - return match radix { - 2 { - integer_from_special_string(characters, 1) - } - 16 { - integer_from_special_string(characters, 4) - } - else { - integer_from_regular_string(characters, radix) - } - } + validate_string(characters, radix)! + return integer_from_regular_string(characters, radix) } @[direct_array_access] @@ -217,48 +219,6 @@ fn validate_string(characters string, radix u32) ! { } } -@[direct_array_access] -fn integer_from_special_string(characters string, chunk_size int) Integer { - sign_present := characters[0] == `+` || characters[0] == `-` - - signum := if sign_present { - if characters[0] == `-` { -1 } else { 1 } - } else { - 1 - } - - start_index := if sign_present { 1 } else { 0 } - - mut big_digits := []u32{cap: ((characters.len * chunk_size) >> 5) + 1} - mut current := u32(0) - mut offset := 0 - for index := characters.len - 1; index >= start_index; index-- { - digit := characters[index] - value := u32(digit_array.index(digit)) - - current |= value << offset - offset += chunk_size - - if offset == 32 { - big_digits << current - current = u32(0) - offset = 0 - } - } - - // Store the accumulated value into the digit array - if current != 0 { - big_digits << current - } - - shrink_tail_zeros(mut big_digits) - - return Integer{ - digits: big_digits - signum: if big_digits.len == 0 { 0 } else { signum } - } -} - @[direct_array_access] fn integer_from_regular_string(characters string, radix u32) Integer { sign_present := characters.len > 0 && (characters[0] == `+` || characters[0] == `-`) @@ -273,17 +233,13 @@ fn integer_from_regular_string(characters string, radix u32) Integer { mut result := zero_int radix_int := integer_from_u32(radix) - pow := radix_options[int(radix)] - radix_pow := radix_int.pow(u32(pow)) - for i := start_index; i < characters.len; i += pow { - end := math.min(i + pow, characters.len) - num_str := characters[i..end] - if num_str.len == pow { - result *= radix_pow - } else { - result *= radix_int.pow(u32(num_str.len)) - } - result += integer_from_u32(regular_string_to_radix(num_str, radix)) + + for index := start_index; index < characters.len; index++ { + digit := characters[index] + value := digit_array.index(digit) + + result *= radix_int + result += integer_from_int(value) } return Integer{ @@ -292,15 +248,6 @@ fn integer_from_regular_string(characters string, radix u32) Integer { } } -fn regular_string_to_radix(characters string, radix u32) u32 { - mut result := u32(0) - - for c in characters { - result = result * radix + u32(digit_array.index(c)) - } - return result -} - // abs returns the absolute value of the integer `a`. pub fn (a Integer) abs() Integer { return if a.signum == 0 { @@ -366,7 +313,7 @@ pub fn (minuend Integer) - (subtrahend Integer) Integer { fn (integer Integer) add(addend Integer) Integer { a := integer.digits b := addend.digits - mut storage := []u32{len: imax(a.len, b.len) + 1} + mut storage := []u64{len: imax(a.len, b.len) + 1} add_digit_array(a, b, mut storage) return Integer{ signum: integer.signum @@ -380,7 +327,7 @@ fn (integer Integer) subtract(subtrahend Integer) Integer { return zero_int } a, b := if cmp > 0 { integer, subtrahend } else { subtrahend, integer } - mut storage := []u32{len: a.digits.len} + mut storage := []u64{len: a.digits.len} subtract_digit_array(a.digits, b.digits, mut storage) return Integer{ signum: cmp * a.signum @@ -401,7 +348,7 @@ pub fn (multiplicand Integer) * (multiplier Integer) Integer { return multiplicand.clone() } // The final sign is the product of the signs - mut storage := []u32{len: multiplicand.digits.len + multiplier.digits.len} + mut storage := []u64{len: multiplicand.digits.len + multiplier.digits.len} multiply_digit_array(multiplicand.digits, multiplier.digits, mut storage) return Integer{ signum: multiplicand.signum * multiplier.signum @@ -415,8 +362,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 := []u32{cap: int_max(1, dividend.digits.len - divisor.digits.len + 1)} - mut r := []u32{cap: dividend.digits.len} + mut q := []u64{cap: int_max(1, dividend.digits.len - divisor.digits.len + 1)} + mut r := []u64{cap: dividend.digits.len} mut q_signum := 0 mut r_signum := 0 @@ -545,8 +492,8 @@ fn (a Integer) mask_bits(n u32) Integer { return zero_int } - w := n / 32 - b := n % 32 + w := n / digit_bits + b := n % digit_bits if w >= a.digits.len { return a @@ -554,17 +501,17 @@ fn (a Integer) mask_bits(n u32) Integer { return Integer{ digits: if b == 0 { - mut storage := []u32{len: int(w)} + mut storage := []u64{len: int(w)} for i := 0; i < storage.len; i++ { storage[i] = a.digits[i] } storage } else { - mut storage := []u32{len: int(w) + 1} + mut storage := []u64{len: int(w) + 1} for i := 0; i < storage.len; i++ { storage[i] = a.digits[i] } - storage[w] &= ~(u32(-1) << b) + storage[w] &= ~(u64(-1) << b) storage } signum: 1 @@ -593,7 +540,7 @@ pub fn (base Integer) pow(exponent u32) Integer { } // mod_pow returns the integer `base` raised to the power of the u32 `exponent` modulo the integer `modulus`. -pub fn (base Integer) mod_pow(exponent u32, modulus Integer) Integer { +pub fn (base Integer) mod_pow(exponent u64, modulus Integer) Integer { if exponent == 0 { return one_int } @@ -707,8 +654,8 @@ pub fn (a Integer) < (b Integer) bool { // get_bit checks whether the bit at the given index is set. @[direct_array_access] pub fn (a Integer) get_bit(i u32) bool { - target_index := i / 32 - offset := i % 32 + target_index := i / digit_bits + offset := i % digit_bits if target_index >= a.digits.len { return false } @@ -717,8 +664,8 @@ pub fn (a Integer) get_bit(i u32) bool { // set_bit sets the bit at the given index to the given value. pub fn (mut a Integer) set_bit(i u32, value bool) { - target_index := i / 32 - offset := i % 32 + target_index := i / digit_bits + offset := i % digit_bits if target_index >= a.digits.len { if value { @@ -730,9 +677,9 @@ pub fn (mut a Integer) set_bit(i u32, value bool) { mut copy := a.digits.clone() if value { - copy[target_index] |= 1 << offset + copy[target_index] |= u64(1) << offset } else { - copy[target_index] &= ~(1 << offset) + copy[target_index] &= ~(u64(1) << offset) } a = Integer{ @@ -745,7 +692,7 @@ pub fn (mut a Integer) set_bit(i u32, value bool) { // // Note: both operands are treated as absolute values. pub fn (a Integer) bitwise_or(b Integer) Integer { - mut result := []u32{len: imax(a.digits.len, b.digits.len)} + mut result := []u64{len: imax(a.digits.len, b.digits.len)} bitwise_or_digit_array(a.digits, b.digits, mut result) return Integer{ digits: result @@ -757,7 +704,7 @@ pub fn (a Integer) bitwise_or(b Integer) Integer { // // Note: both operands are treated as absolute values. pub fn (a Integer) bitwise_and(b Integer) Integer { - mut result := []u32{len: imax(a.digits.len, b.digits.len)} + mut result := []u64{len: imax(a.digits.len, b.digits.len)} bitwise_and_digit_array(a.digits, b.digits, mut result) return Integer{ digits: result @@ -769,7 +716,7 @@ pub fn (a Integer) bitwise_and(b Integer) Integer { // // Note: the integer is treated as an absolute value. pub fn (a Integer) bitwise_not() Integer { - mut result := []u32{len: a.digits.len} + mut result := []u64{len: a.digits.len} bitwise_not_digit_array(a.digits, mut result) return Integer{ digits: result @@ -792,7 +739,7 @@ pub fn (a Integer) bitwise_com() Integer { // // Note: both operands are treated as absolute values. pub fn (a Integer) bitwise_xor(b Integer) Integer { - mut result := []u32{len: imax(a.digits.len, b.digits.len)} + mut result := []u64{len: imax(a.digits.len, b.digits.len)} bitwise_xor_digit_array(a.digits, b.digits, mut result) return Integer{ digits: result @@ -809,9 +756,9 @@ pub fn (a Integer) left_shift(amount u32) Integer { if amount == 0 { return a } - normalised_amount := amount & 31 - digit_offset := int(amount >> 5) - mut new_array := []u32{len: a.digits.len + digit_offset} + normalised_amount := amount % digit_bits + digit_offset := int(amount / digit_bits) + mut new_array := []u64{len: a.digits.len + digit_offset} for index in 0 .. a.digits.len { new_array[index + digit_offset] = a.digits[index] } @@ -833,12 +780,12 @@ pub fn (a Integer) right_shift(amount u32) Integer { if amount == 0 { return a } - normalised_amount := amount & 31 - digit_offset := int(amount >> 5) + normalised_amount := amount % digit_bits + digit_offset := int(amount / digit_bits) if digit_offset >= a.digits.len { return zero_int } - mut new_array := []u32{len: a.digits.len - digit_offset} + mut new_array := []u64{len: a.digits.len - digit_offset} for index in 0 .. new_array.len { new_array[index] = a.digits[index + digit_offset] } @@ -854,46 +801,13 @@ pub fn (a Integer) right_shift(amount u32) Integer { // bin_str returns the binary string representation of the integer `a`. @[direct_array_access] pub fn (integer Integer) bin_str() string { - // We have the zero integer - if integer.signum == 0 { - return '0' - } - // Add the sign if present - sign_needed := integer.signum == -1 - mut result_builder := strings.new_builder(integer.bit_len() + if sign_needed { 1 } else { 0 }) - if sign_needed { - result_builder.write_string('-') - } - - result_builder.write_string(u32_to_binary_without_lz(integer.digits[integer.digits.len - 1])) - - for index := integer.digits.len - 2; index >= 0; index-- { - result_builder.write_string(u32_to_binary_with_lz(integer.digits[index])) - } - return result_builder.str() + return integer.radix_str(2) } // hex returns the hexadecimal string representation of the integer `a`. @[direct_array_access] pub fn (integer Integer) hex() string { - // We have the zero integer - if integer.signum == 0 { - return '0' - } - // Add the sign if present - sign_needed := integer.signum == -1 - mut result_builder := strings.new_builder(integer.digits.len * 8 + - if sign_needed { 1 } else { 0 }) - if sign_needed { - result_builder.write_string('-') - } - - result_builder.write_string(u32_to_hex_without_lz(integer.digits[integer.digits.len - 1])) - - for index := integer.digits.len - 2; index >= 0; index-- { - result_builder.write_string(u32_to_hex_with_lz(integer.digits[index])) - } - return result_builder.str() + return integer.radix_str(16) } // radix_str returns the string representation of the integer `a` in the specified radix. @@ -901,67 +815,32 @@ pub fn (integer Integer) radix_str(radix u32) string { if integer.signum == 0 || radix == 0 { return '0' } - return match radix { - 2 { - integer.bin_str() - } - 16 { - integer.hex() - } - else { - integer.general_radix_str(int(radix)) - } - } + return integer.general_radix_str(radix) } -fn (integer Integer) general_radix_str(radix int) string { +fn (integer Integer) general_radix_str(radix u32) string { $if debug { assert radix != 0 } - divisor := integer_from_int(radix).pow(u32(radix_options[radix])) + divisor := integer_from_u32(radix) mut current := integer.abs() mut new_current := zero_int mut digit := zero_int - mut sb := strings.new_builder(integer.digits.len * radix_options[radix]) - mut st := []string{cap: integer.digits.len * radix_options[radix]} + mut rune_array := []rune{cap: current.digits.len * 4} for current.signum > 0 { new_current, digit = current.div_mod_internal(divisor) - st << general_str(new_current, digit, radix) + rune_array << digit_array[digit.int()] + unsafe { digit.free() } + unsafe { current.free() } current = new_current } if integer.signum == -1 { - sb.write_string('-') - } - for st.len > 0 { - sb.write_string(st.pop()) - } - return sb.str() -} - -fn general_str(quotient Integer, remainder Integer, radix int) string { - if quotient.signum == 0 && remainder.signum == 0 { - return '0' + rune_array << `-` } - divisor := integer_from_int(radix) - mut current := remainder.abs() - mut new_current := zero_int - mut digit := zero_int - mut sb := strings.new_builder(radix_options[radix]) - mut st := []u8{cap: radix_options[radix]} - for current.signum > 0 { - new_current, digit = current.div_mod_internal(divisor) - st << digit_array[digit.int()] - current = new_current - } - if quotient.signum > 0 { - sb.write_string(strings.repeat(48, radix_options[radix] - st.len)) - } - for st.len > 0 { - sb.write_u8(st.pop()) - } - return sb.str() + rune_array.reverse_in_place() + return rune_array.string() } // str returns the decimal string representation of the integer `a`. @@ -969,34 +848,6 @@ pub fn (integer Integer) str() string { return integer.radix_str(10) } -fn u32_to_binary_without_lz(value u32) string { - return strconv.format_uint(value, 2) -} - -fn u32_to_binary_with_lz(value u32) string { - mut result_builder := strings.new_builder(32) - binary_result := strconv.format_uint(value, 2) - - result_builder.write_string(strings.repeat(`0`, 32 - binary_result.len)) - result_builder.write_string(binary_result) - - return result_builder.str() -} - -fn u32_to_hex_without_lz(value u32) string { - return strconv.format_uint(value, 16) -} - -fn u32_to_hex_with_lz(value u32) string { - mut result_builder := strings.new_builder(8) - hex_result := strconv.format_uint(value, 16) - - result_builder.write_string(strings.repeat(`0`, 8 - hex_result.len)) - result_builder.write_string(hex_result) - - return result_builder.str() -} - // int returns the integer value of the integer `a`. // NOTE: This may cause loss of precision. @[direct_array_access] @@ -1005,7 +856,7 @@ pub fn (a Integer) int() int { return 0 } // Check for minimum value int - if a.digits[0] == 2147483648 && a.signum == -1 { + if a.digits[0] >= 2147483648 && a.signum == -1 { return -2147483648 } // Rest of the values should be fine @@ -1020,25 +871,45 @@ pub fn (a Integer) bytes() ([]u8, int) { if a.signum == 0 { return []u8{len: 0}, 0 } - mut result := []u8{cap: a.digits.len * 4} - mut mask := u32(0xff000000) - mut offset := 24 - mut non_zero_found := false - for index := a.digits.len - 1; index >= 0; { - value := u8((a.digits[index] & mask) >> offset) - non_zero_found = non_zero_found || value != 0 - if non_zero_found { - result << value + bit_len := a.bit_len() + mut bytes := []u8{cap: bit_len / 8 + 1} + mut current_byte := u8(0) + mut bits_in_byte := 0 + mut digit := a.digits.last() + mut bit := u8(0) + + // pad first byte + bits_in_byte = 8 - bit_len % 8 + if bits_in_byte == 8 { + bits_in_byte = 0 + } + // MSB digit + for i := bit_len % digit_bits - 1; i >= 0; i-- { + bit = u8((digit >> i) & 1) + current_byte = (current_byte << 1) | u8(bit) + bits_in_byte++ + if bits_in_byte == 8 { + bytes << current_byte + current_byte = 0 + bits_in_byte = 0 } - mask >>= 8 - offset -= 8 - if offset < 0 { - mask = u32(0xff000000) - offset = 24 - index-- + } + + for i := a.digits.len - 2; i >= 0; i-- { + digit = a.digits[i] + for shift := digit_bits - 1; shift >= 0; shift-- { + bit = u8((digit >> shift) & 1) + current_byte = (current_byte << 1) | bit + bits_in_byte++ + if bits_in_byte == 8 { + bytes << current_byte + current_byte = 0 + bits_in_byte = 0 + } } } - return result, a.signum + + return bytes, a.signum } // factorial returns the factorial of the integer `a`. @@ -1221,7 +1092,7 @@ fn (a Integer) mod_inv(m Integer) Integer { tmp := if q == one_int { x } else if q.digits.len == 1 && q.digits[0] & (q.digits[0] - 1) == 0 { - x.left_shift(u32(bits.trailing_zeros_32(q.digits[0]))) + x.left_shift(u32(bits.trailing_zeros_64(q.digits[0]))) } else { q * x } + y @@ -1259,7 +1130,7 @@ fn (x Integer) rsh_to_set_bit() (Integer, u32) { for x.digits[n] == 0 { n++ } - n = (n << 5) + u32(bits.trailing_zeros_32(x.digits[n])) + n = (n * digit_bits) + u32(bits.trailing_zeros_64(x.digits[n])) return x.right_shift(n), n } @@ -1284,7 +1155,7 @@ pub fn (x Integer) is_power_of_2() bool { } } n := x.digits.last() - return n & (n - u32(1)) == 0 + return n & (n - u64(1)) == 0 } // bit_len returns the number of bits required to represent the integer `a`. @@ -1296,5 +1167,5 @@ pub fn (x Integer) bit_len() int { if x.digits.len == 0 { return 0 } - return x.digits.len * 32 - bits.leading_zeros_32(x.digits.last()) + return x.digits.len * digit_bits - (bits.leading_zeros_64(x.digits.last()) - (64 - digit_bits)) } diff --git a/vlib/math/big/special_array_ops.v b/vlib/math/big/special_array_ops.v index c374e7d05..f86be3133 100644 --- a/vlib/math/big/special_array_ops.v +++ b/vlib/math/big/special_array_ops.v @@ -3,7 +3,7 @@ module big import strings @[direct_array_access; inline] -fn shrink_tail_zeros(mut a []u32) { +fn shrink_tail_zeros(mut a []u64) { mut alen := a.len for alen > 0 && a[alen - 1] == 0 { alen-- @@ -28,7 +28,7 @@ fn (i &Integer) shrink_tail_zeros() { // Both quotient and remaider are already allocated but of length 0 // TODO: the manualfree tag here is a workaround for compilation with -autofree. Remove it, when the -autofree bug is fixed. @[manualfree] -fn newton_divide_array_by_array(operand_a []u32, operand_b []u32, mut quotient []u32, mut remainder []u32) { +fn newton_divide_array_by_array(operand_a []u64, operand_b []u64, mut quotient []u64, mut remainder []u64) { // transform back to Integers (on the stack without allocation) a := Integer{ signum: 1 @@ -42,7 +42,10 @@ fn newton_divide_array_by_array(operand_a []u32, operand_b []u32, mut quotient [ k := a.bit_len() + b.bit_len() // a*b < 2**k mut x := integer_from_int(2) // 0 < x < 2**(k+1)/b // initial guess for convergence // https://en.wikipedia.org/wiki/Division_algorithm#Newton%E2%80%93Raphson_division - // use 48/17 - 32/17.D (divisor) + // use 48/17 - (32/17)*D (divisor) + // where D - Divisor B adjusted to fit 0.5-1 range by replacing the exponent field with 8'd126 + // 0x0f0f0f0f0f0f0f0f * 0x11 == 0xffffffff_ffffffff = -1 + // 0x0f0f0f0f0f0f0f0f = -1/17 initial_guess := (((integer_from_int(48) - (integer_from_int(32) * b)) * integer_from_i64(0x0f0f0f0f0f0f0f0f)).right_shift(64)).neg() // / 17 == 0x11 if initial_guess > zero_int { x = initial_guess @@ -72,8 +75,9 @@ fn newton_divide_array_by_array(operand_a []u32, operand_b []u32, mut quotient [ shrink_tail_zeros(mut remainder) } -@[direct_array_access; inline] -fn debug_u32_str(a []u32) string { +// debug_u64_str output a `[]u64` +@[direct_array_access] +fn debug_u64_str(a []u64) string { mut sb := strings.new_builder(30) sb.write_string('[') mut first := true @@ -81,7 +85,50 @@ fn debug_u32_str(a []u32) string { if !first { sb.write_string(', ') } - sb.write_string('0x${a[i].hex()}') + sb.write_string('0x${a[i].hex():016}') + first = false + } + sb.write_string(']') + return sb.str() +} + +// debug_u32_str for 32bit bignum test only, convert a `[]u64` to `[]u32`. +@[direct_array_access] +fn debug_u32_str(a []u64) string { + mut b := []u32{cap: a.len * 2} + mut curr_u32 := u32(0) + mut bits_collected := 0 + for w in a { + for i in 0 .. digit_bits { + bit := (w >> i) & 1 + curr_u32 |= u32(bit) << bits_collected + bits_collected++ + if bits_collected == 32 { + b << curr_u32 + curr_u32 = 0 + bits_collected = 0 + } + } + } + if bits_collected > 0 { + b << curr_u32 + } + + mut blen := b.len + for blen > 0 && b[blen - 1] == 0 { + blen-- + } + unsafe { + b.len = blen + } + mut sb := strings.new_builder(30) + sb.write_string('[') + mut first := true + for i in 0 .. b.len { + if !first { + sb.write_string(', ') + } + sb.write_string('0x${b[i].hex():08}') first = false } sb.write_string(']') @@ -89,7 +136,7 @@ fn debug_u32_str(a []u32) string { } @[direct_array_access; inline] -fn found_multiplication_base_case(operand_a []u32, operand_b []u32, mut storage []u32) bool { +fn found_multiplication_base_case(operand_a []u64, operand_b []u64, mut storage []u64) bool { // base case necessary to end recursion if operand_a.len == 0 || operand_b.len == 0 { storage.clear() @@ -112,7 +159,7 @@ fn found_multiplication_base_case(operand_a []u32, operand_b []u32, mut storage // possible optimisations: // - transform one or all the recurrences in loops @[direct_array_access] -fn karatsuba_multiply_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) { +fn karatsuba_multiply_digit_array(operand_a []u64, operand_b []u64, mut storage []u64) { if found_multiplication_base_case(operand_a, operand_b, mut storage) { return } @@ -121,14 +168,14 @@ fn karatsuba_multiply_digit_array(operand_a []u32, operand_b []u32, mut storage half := imax(operand_a.len, operand_b.len) / 2 mut a_l := unsafe { operand_a[0..half] } mut a_h := unsafe { operand_a[half..] } - mut b_l := []u32{} - mut b_h := []u32{} + mut b_l := []u64{} + mut b_h := []u64{} if half <= operand_b.len { b_l = unsafe { operand_b[0..half] } b_h = unsafe { operand_b[half..] } } else { b_l = unsafe { operand_b } - // b_h = []u32{} + // b_h = []u64{} } shrink_tail_zeros(mut a_l) shrink_tail_zeros(mut a_h) @@ -138,15 +185,15 @@ fn karatsuba_multiply_digit_array(operand_a []u32, operand_b []u32, mut storage // use storage for p_1 to avoid allocation and copy later multiply_digit_array(a_h, b_h, mut storage) - mut p_3 := []u32{len: a_l.len + b_l.len + 1} + mut p_3 := []u64{len: a_l.len + b_l.len + 1} multiply_digit_array(a_l, b_l, mut p_3) - mut tmp_1 := []u32{len: imax(a_h.len, a_l.len) + 1} - mut tmp_2 := []u32{len: imax(b_h.len, b_l.len) + 1} + mut tmp_1 := []u64{len: imax(a_h.len, a_l.len) + 1} + mut tmp_2 := []u64{len: imax(b_h.len, b_l.len) + 1} add_digit_array(a_h, a_l, mut tmp_1) add_digit_array(b_h, b_l, mut tmp_2) - mut p_2 := []u32{len: operand_a.len + operand_b.len + 1} + mut p_2 := []u64{len: operand_a.len + operand_b.len + 1} multiply_digit_array(tmp_1, tmp_2, mut p_2) subtract_in_place(mut p_2, storage) // p_1 subtract_in_place(mut p_2, p_3) @@ -162,14 +209,14 @@ fn karatsuba_multiply_digit_array(operand_a []u32, operand_b []u32, mut storage // TODO: the manualfree tag here is a workaround for compilation with -autofree. Remove it, when the -autofree bug is fixed. @[direct_array_access; manualfree] -fn toom3_multiply_digit_array(operand_a []u32, operand_b []u32, mut storage []u32) { +fn toom3_multiply_digit_array(operand_a []u64, operand_b []u64, mut storage []u64) { if found_multiplication_base_case(operand_a, operand_b, mut storage) { return } // After the base case, we have operand_a as the larger integer in terms of digit length - // k is the length (in u32 digits) of the lower order slices + // k is the length (in u64 digits) of the lower order slices k := (operand_a.len + 2) / 3 k2 := 2 * k @@ -214,32 +261,32 @@ fn toom3_multiply_digit_array(operand_a []u32, operand_b []u32, mut storage []u3 } else if operand_b.len < k2 { if !operand_b[..k].all(it == 0) { b0 = Integer{ - digits: operand_b[..k] + digits: operand_b[..k].clone() signum: 1 } } b0.shrink_tail_zeros() b1 = Integer{ - digits: operand_b[k..] + digits: operand_b[k..].clone() signum: 1 } } else { if !operand_b[..k].all(it == 0) { b0 = Integer{ - digits: operand_b[..k] + digits: operand_b[..k].clone() signum: 1 } } b0.shrink_tail_zeros() if !operand_b[k..k2].all(it == 0) { b1 = Integer{ - digits: operand_b[k..k2] + digits: operand_b[k..k2].clone() signum: 1 } } b1.shrink_tail_zeros() b2 = Integer{ - digits: operand_b[k2..] + digits: operand_b[k2..].clone() signum: 1 } } @@ -266,7 +313,7 @@ fn toom3_multiply_digit_array(operand_a []u32, operand_b []u32, mut storage []u3 tm1 = tm1 - t2 // shift amount - s := u32(k) << 5 + s := u32(k) * digit_bits result := (((pinf.left_shift(s) + t2).left_shift(s) + t1).left_shift(s) + tm1).left_shift(s) + p0 @@ -276,7 +323,7 @@ fn toom3_multiply_digit_array(operand_a []u32, operand_b []u32, mut storage []u3 @[inline] fn pow2(k int) Integer { - mut ret := []u32{len: (k >> 5) + 1} + mut ret := []u64{len: (k / digit_bits) + 1} bit_set(mut ret, k) return Integer{ signum: 1 @@ -285,9 +332,9 @@ fn pow2(k int) Integer { } // optimized left shift in place. amount must be positive -fn left_shift_digits_in_place(mut a []u32, amount int) { +fn left_shift_digits_in_place(mut a []u64, amount int) { // this is actual in builtin/array.v, prepend_many (private fn) - // x := []u32{ len : amount } + // x := []u64{ len : amount } // a.prepend_many(&x[0], amount) old_len := a.len elem_size := a.element_size @@ -301,14 +348,14 @@ fn left_shift_digits_in_place(mut a []u32, amount int) { } // optimized right shift in place. amount must be positive -fn right_shift_digits_in_place(mut a []u32, amount int) { +fn right_shift_digits_in_place(mut a []u64, amount int) { a.drop(amount) } // operand b can be greater than operand a // the capacity of both array is supposed to be sufficient @[direct_array_access; inline] -fn add_in_place(mut a []u32, b []u32) { +fn add_in_place(mut a []u64, b []u64) { len_a := a.len len_b := b.len max := imax(len_a, len_b) @@ -316,30 +363,30 @@ fn add_in_place(mut a []u32, b []u32) { mut carry := u64(0) for index in 0 .. min { partial := carry + a[index] + b[index] - a[index] = u32(partial) - carry = u32(partial >> 32) + a[index] = u64(partial) & max_digit + carry = u64(partial >> digit_bits) } if len_a >= len_b { for index in min .. max { partial := carry + a[index] - a[index] = u32(partial) - carry = u32(partial >> 32) + a[index] = u64(partial) & max_digit + carry = u64(partial >> digit_bits) } } else { for index in min .. max { partial := carry + b[index] - a << u32(partial) - carry = u32(partial >> 32) + a << u64(partial) & max_digit + carry = u64(partial >> digit_bits) } } if carry > 0 { - a << u32(carry) + a << carry } } // a := a - b supposed a >= b @[direct_array_access; inline] -fn subtract_in_place(mut a []u32, b []u32) { +fn subtract_in_place(mut a []u64, b []u64) { len_a := a.len len_b := b.len max := imax(len_a, len_b) @@ -347,24 +394,24 @@ fn subtract_in_place(mut a []u32, b []u32) { mut borrow := false for index in 0 .. min { - mut a_digit := u64(a[index]) + mut a_digit := a[index] b_digit := b[index] + if borrow { u64(1) } else { u64(0) } borrow = a_digit < b_digit if borrow { - a_digit += 0x100000000 + a_digit = a_digit | (u64(1) << digit_bits) } - a[index] = u32(a_digit - b_digit) + a[index] = a_digit - b_digit } if len_a >= len_b { for index in min .. max { - mut a_digit := u64(a[index]) + mut a_digit := a[index] b_digit := if borrow { u64(1) } else { u64(0) } borrow = a_digit < b_digit if borrow { - a_digit += 0x100000000 + a_digit = a_digit | (u64(1) << digit_bits) } - a[index] = u32(a_digit - b_digit) + a[index] = a_digit - b_digit } } else { // if len.b > len.a return zero a.clear() diff --git a/vlib/math/big/special_array_ops_test.v b/vlib/math/big/special_array_ops_test.v index 21a2c284a..fe3f51186 100644 --- a/vlib/math/big/special_array_ops_test.v +++ b/vlib/math/big/special_array_ops_test.v @@ -1,75 +1,75 @@ module big fn test_add_in_place() { - mut a := [u32(1), 2, 3] - mut b := [u32(5), 6, 7] + mut a := [u64(1), 2, 3] + mut b := [u64(5), 6, 7] add_in_place(mut a, b) - assert a == [u32(6), 8, 10] + assert a == [u64(6), 8, 10] - a = [u32(11), 10, 11, 12] - b = [u32(1), 2] + a = [u64(11), 10, 11, 12] + b = [u64(1), 2] add_in_place(mut a, b) - assert a == [u32(12), 12, 11, 12] + assert a == [u64(12), 12, 11, 12] - a = []u32{cap: 4} - a << u32(1) - a << u32(2) - b = [u32(3), 4, 5, 6] + a = []u64{cap: 4} + a << u64(1) + a << u64(2) + b = [u64(3), 4, 5, 6] add_in_place(mut a, b) - assert a == [u32(4), 6, 5, 6] + assert a == [u64(4), 6, 5, 6] - a = [u32(0x3ce9124b), 0x1438] - b = [u32(0xdb166062)] + a = [u64(0x0423efab_3ce9124b), 0x1438] + b = [u64(0x0d1dd122_db166062)] add_in_place(mut a, b) - assert a == [u32(0x17ff72ad), 0x1439] + assert a == [u64(0x0141c0ce_17ff72ad), 0x1439] } fn test_left_shift_digits_in_place() { - mut a := [u32(5), 6, 7, 8] + mut a := [u64(5), 6, 7, 8] left_shift_digits_in_place(mut a, 2) - assert a == [u32(0), 0, 5, 6, 7, 8] + assert a == [u64(0), 0, 5, 6, 7, 8] } fn test_multiply_karatsuba_01() { - mut a := [u32(3)] - mut b := []u32{} - mut c := []u32{len: a.len + b.len, init: 0} + mut a := [u64(3)] + mut b := []u64{} + mut c := []u64{len: a.len + b.len, init: 0} karatsuba_multiply_digit_array(a, b, mut c) - assert c == []u32{} + assert c == []u64{} - a = []u32{} - b = [u32(4)] - c = []u32{len: a.len + b.len, init: 0} + a = []u64{} + b = [u64(4)] + c = []u64{len: a.len + b.len, init: 0} karatsuba_multiply_digit_array(a, b, mut c) - assert c == []u32{} + assert c == []u64{} - a = [u32(3)] - b = [u32(1)] - c = []u32{len: a.len + b.len, init: 0} + a = [u64(3)] + b = [u64(1)] + c = []u64{len: a.len + b.len, init: 0} karatsuba_multiply_digit_array(a, b, mut c) assert c == a - a = [u32(1)] - b = [u32(5)] - c = []u32{len: a.len + b.len, init: 0} + a = [u64(1)] + b = [u64(5)] + c = []u64{len: a.len + b.len, init: 0} karatsuba_multiply_digit_array(a, b, mut c) assert c == b - a = [u32(1234)] - b = [u32(567)] - c = []u32{len: a.len + b.len + 1, init: 0} + a = [u64(1234)] + b = [u64(567)] + c = []u64{len: a.len + b.len + 1, init: 0} karatsuba_multiply_digit_array(a, b, mut c) - assert c == [u32(699_678)] + assert c == [u64(699_678)] - a = [u32(0x17ff72ad), 0x1439] - b = [u32(0x30df2ea6)] - c = []u32{len: a.len + b.len + 1, init: 0} + a = [u64(0x17ff72ad), 0x1439] + b = [u64(0x30df2ea6)] + c = []u64{len: a.len + b.len + 1, init: 0} karatsuba_multiply_digit_array(a, b, mut c) - assert c == [u32(0xcaf2722e), 0x55eb2c5a, 0x3dc] + assert c == [u64(0x0494d164_caf2722e), 0x3dc_51565af6] a_operand := integer_from_string('95484736384949495947362') or { panic(err) } b_operand := integer_from_string('39474638493') or { panic(err) } - c = []u32{len: a_operand.digits.len + b_operand.digits.len, init: 0} + c = []u64{len: a_operand.digits.len + b_operand.digits.len, init: 0} karatsuba_multiply_digit_array(a_operand.digits, b_operand.digits, mut c) expected := integer_from_string('3769225450395285038584683507005466') or { panic(err) } assert c == expected.digits @@ -82,7 +82,7 @@ fn test_multiply_karatsuba_02() { b := integer_from_string('977091076523237491790970633699383779582771973038531457285598238843271083830214915826312193418602834034688531898668229388286706296786321423078510899614439367') or { panic(err) } - mut c := []u32{len: a.digits.len + b.digits.len + 1, init: 0} + mut c := []u64{len: a.digits.len + b.digits.len + 1, init: 0} karatsuba_multiply_digit_array(a.digits, b.digits, mut c) expected := integer_from_string('52348074924977237255285644820010078601114587486470740900886892189662650320988400136613780986308710610258879824881256666730655821800564143426560480113864123642197317383052431412305975584645367703594190956925565749714310612399025459615546540332117815550470167143256687163102859337019449165214274088466835988832405507818643018779158891710706073875995722420460085755') or { panic(err) @@ -93,7 +93,7 @@ fn test_multiply_karatsuba_02() { fn test_multiply_karatsuba_03() { a := integer_from_string('9729117383001812976929423642') or { panic(err) } b := integer_from_string('36889625830') or { panic(err) } - mut c := []u32{len: a.digits.len + b.digits.len + 1, init: 0} + mut c := []u64{len: a.digits.len + b.digits.len + 1, init: 0} karatsuba_multiply_digit_array(a.digits, b.digits, mut c) expected := integer_from_string('358903499915085682930564860470935872860') or { panic(err) } assert c == expected.digits @@ -106,7 +106,7 @@ fn test_multiply_karatsuba_04() { b := integer_from_string('3731153765349978524337790303019073866169615820236360525860598113101007652419849831445556596548560058884406106275477847069475858381339095110646794630677741924575724730195310383260858937531720778415641064228137265209604913544897683432011584846544175241128160109960448428720014016539849270207788214347700420341097980912654975389446512856552635037966995033809488867690903720129762695314322725472787603000346202896778207987391285985804181406019619338215824347026372158857600803429407141531830203430317068399662548922504444591385274908074684241297752512843236588392255573848392016559311506271535791162151101220135017417159573887941542933118000377620518732836102365066672711767692414034390790360321914677167248083633551615232080788875425342060374021891940380358579307475843181010943250359217102987258298') or { panic(err) } - mut c := []u32{len: a.digits.len + b.digits.len + 1, init: 0} + mut c := []u64{len: a.digits.len + b.digits.len + 1, init: 0} karatsuba_multiply_digit_array(a.digits, b.digits, mut c) expected := integer_from_string('13048451769827723841923278321193855241872058428961713256963598223823113544308869393875481721028040256602323122921447247203223525228006623227671899134332670311746731440890333942554337434694808695934023689467815943095478675997839202402226439648752044627367822829114336020280929447282163468585404769209994615771773294733509244673189891278506538595789870783779458539371487749876399583947052094519945323844069814095194450349814584625148870057340788686115533159501329566111170880650412359779760456808393370185410432176572682144365401835368082749637208673417028234063907368810782078088911618126690572051133815667647041419960704854062060909165346545952623192359395985743203202132668410599975580175587551065836786982237372286798893118604373993845408861570111075107858370728022568891638898121296878893306306949816814132427209493260832613019444633237775219695833279678133602467177003429529917049633056431504426776487440348449903443579642102845765312131382153105184606682393976592634403346380092701679847372217955349096619505179079762928343558899980175648735397570586605980720828687524561573123038494159647318551197545496201465276098212580839445369200418315000310508759594006187422953054821725336960756926970677282806506256927313164107160898391061464246294174106905930356469329214907846583235734490370187185094462047618230289621171427228984027756484384136106088551265698775255811805670865291484911854071447316265726194982034568878612017100992037975737759552371045887134952857158066446970434610533685859910727214440334704089405629838823453640146362193811088002789803040768301613468217974726406918247619138081429486763003592049227953605769429261539079312619553867710373453222399692311140251766121447367352413232353894344679702619305784387721023140535029791620') or { panic(err) @@ -121,7 +121,7 @@ fn test_multiply_karatsuba_05() { b := integer_from_string('6131999155203247409089965800322203420453654507953806420664043731807402839483789086591425043878906721774222779997525063107602234437085033473124680861190605857734959468536446941493399655926830326691521621629120060373092027010394244758024141037091879171132413714138502869706470511932123207420626930454253770927552857807363136786489394763374159969607962219873240941060072539904794844154488056908195189396404340107989455229902431972663610868653304103396580746811591826044603465964167164178855966898012392223090630560777990379102734550082243603230894035356131452642667384300201186297016843973905130369865474318465241408259975202340234896078010895740058413587215886548145454236789540259643610072467034822757986309311709163062352204171299868075798237348977707372770035406598475786685283610490871427030552271811065910679040103566128693521437169494686602257958511710409941630744336034132463653071078769355257935175422990452306401119325390221264671263162134241884024846471010012913485967133615252263790910226344828197924602249009133593805769333128797013042594573917849267639752846975907650362497446525930636106580780824323557414269051915905993675310076325449434427323319596509995607027978958874919018609350403661403557093530969901119390174552941500745876099234877620450622478016638926607899983456508753736077292767888656321808232398065656946197839116825276017609730557227380499217716159012979059487418587202959339426112335872029689748557894665457579061382636201338034208646424996597993168272103220242687875466825791773419606359011555661467526139097341579632208557586297052223542776250948355945042246570') or { panic(err) } - mut c := []u32{len: a.digits.len + b.digits.len + 1, init: 0} + mut c := []u64{len: a.digits.len + b.digits.len + 1, init: 0} karatsuba_multiply_digit_array(a.digits, b.digits, mut c) expected := integer_from_string('505848759427364274329357741986216605130929977294192546489509210457601900426556024040073016771820686868026437733614267409521396918447284490340389690485942051363279925801363126432544321078568099151287021886296438619123380168827268412907813594659842321798217885085110398790053409701983461801115224781834023082890662873169191504305955687046526953897067545715710660155497646914961948331236724535698679813769093322003853124540927645162647300288136882363116067195184081346020194278820548735254637692857032965408070350907222436479366121383668103515786877197427062836511944846826932480466793561731962606417667287322005882211404162719492993874592857953397512491054676297738305916085478022723318784939770402209256271554551983534921497675890891312043569145634049147693960076337369503134427352787659326205647773748112716745925158178083868574895622495928425018774779251441610879114247950357970220146481766134464929967308718872255190412874659776500033163244838173195170196127430202488333101394598915495516818921703991533768617658933173526938837897969225258938562812209445135532692302337200334632389376610492768735714190299971580399753641472989662439382540912225551340987306457743400281084947902809687085442657520101417063633057957437817329020615384515166388289331369947059960847133614170679450523810659652782987668106695354035493369448442611758156006527485718905514811065321814359706644699194418190701868505707278984229585722369789551619919906205950348606361748281997782424695650973254173837814677612994239768420808870410390373336826091130971864943060313050478750611821098669520946650257575417245766405176878592155758047549228495747853203792369590850631203019173404179987735732830481798878566582556613534822512685665712663431843718669560669865410543096239869112140888737377028624106518359554511507606387843187232302200517475850782456324239188729330760193745759978024087469060136973622979394674161020956014241092961397194927506142983436789500109128138664821498501821236619524363279857097783139006599000739982697284955344277927297022159974340926721938079595242213432080401390674728369851144911551677647524015332696130827540662024113950127956097630202663581097770816523549988139863573933413116244923387554856120868464947660075936449421302491693593834507504433415086543597651267221343834170846329492147815735045024399026251779389321698355035556629087431166188787768730306310422839817575075342877250911353135209339455830445625001426381558448651172971413260611230589795224559058399088555025826331277804189489960583576819035715077379287444129934501443336913174870788226087158999934636025721931831136227889916792188808832812760406918544186979726578808350730384441425398168139422687286960377555987296985431654502712799904224867124393655074866476410224402965420326358438396015233973741214619430521493940712748773125913358615541385902712067643691280755344861590017706417595031324052737978011814936301852189889555514450129648152568220313555245236446470869647730471752637583335395835914360683070164833354933631599763359799906362664251738450678551934386687346359703853569505426433767046696609137654729946153837937103113817350588833581432259043880508333957686597572923322695933631607569689529867145971190208528802200404369790036623426263214138627080') or { panic(err) @@ -130,47 +130,47 @@ fn test_multiply_karatsuba_05() { } fn test_newton_divide_03() { - a := [u32(0), 4] - b := [u32(0), 1] - mut q := []u32{cap: a.len - b.len + 1} - mut r := []u32{cap: a.len} + a := [u64(0), 4] + b := [u64(0), 1] + mut q := []u64{cap: a.len - b.len + 1} + mut r := []u64{cap: a.len} newton_divide_array_by_array(a, b, mut q, mut r) - assert q == [u32(4)] - assert r == []u32{len: 0} + assert q == [u64(4)] + assert r == []u64{len: 0} } fn test_newton_divide_04() { - a := [u32(2), 4] - b := [u32(0), 1] - mut q := []u32{cap: a.len - b.len + 1} - mut r := []u32{cap: a.len} + a := [u64(2), 4] + b := [u64(0), 1] + mut q := []u64{cap: a.len - b.len + 1} + mut r := []u64{cap: a.len} newton_divide_array_by_array(a, b, mut q, mut r) - assert q == [u32(4)] - assert r == [u32(2)] + assert q == [u64(4)] + assert r == [u64(2)] } fn test_newton_divide_05() { - a := [u32(2), 4, 5] - b := [u32(0), 1] - mut q := []u32{cap: a.len - b.len + 1} - mut r := []u32{cap: a.len} + a := [u64(2), 4, 5] + b := [u64(0), 1] + mut q := []u64{cap: a.len - b.len + 1} + mut r := []u64{cap: a.len} newton_divide_array_by_array(a, b, mut q, mut r) - assert q == [u32(4), 5] - assert r == [u32(2)] + assert q == [u64(4), 5] + assert r == [u64(2)] } fn test_newton_divide_06() { - a := [u32(2), 4, 5, 3] - b := [u32(0), 0x8000] - mut q := []u32{cap: a.len - b.len + 1} - mut r := []u32{cap: a.len} + 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} newton_divide_array_by_array(a, b, mut q, mut r) - assert q == [u32(0xa0000), 0x60000] - assert r == [u32(2), 4] + assert q == [u64(0xa000_0000_0000), 0x6000_0000_0000] + assert r == [u64(2), 4] } fn test_newton_divide_07() { @@ -180,8 +180,8 @@ fn test_newton_divide_07() { b := integer_from_string('977091076523237491790970633699383779582771973038531457285598238843271083830214915826312193418602834034688531898668229388286706296786321423078510899614439367') or { panic(err) } - mut q := []u32{cap: a.digits.len - b.digits.len + 1} - mut r := []u32{cap: a.digits.len} + mut q := []u64{cap: a.digits.len - b.digits.len + 1} + mut r := []u64{cap: a.digits.len} newton_divide_array_by_array(a.digits, b.digits, mut q, mut r) quotient := Integer{ @@ -189,5 +189,5 @@ fn test_newton_divide_07() { digits: q } assert quotient.str() == '53575430359313366047421252453000090528070240585276680372187519418517552556246806124659918940784792906379733645877657341259357264284615702179922887873492874019672838874121154927105373025311855709389770910765' - assert r == [u32(2)] + assert r == [u64(2)] } -- 2.39.5