From f898c84ff9cd2c9412ce6fecb00f4eaabb80ab92 Mon Sep 17 00:00:00 2001 From: Boyu Yang Date: Wed, 27 Mar 2019 02:08:28 +0800 Subject: [PATCH 1/2] Fix bug: in division, when borrow a higher unit, do not forget check if the higher unit was reduced to zero after subtraction. --- .../fixed_uint/core/internal/private_ops.rs | 61 ++++++++++++++++++- fixed-uint-tests/tests/regression_tests.rs | 16 +++++ 2 files changed, 76 insertions(+), 1 deletion(-) diff --git a/constructor/src/fixed_uint/core/internal/private_ops.rs b/constructor/src/fixed_uint/core/internal/private_ops.rs index f250e7e..ea0b7da 100644 --- a/constructor/src/fixed_uint/core/internal/private_ops.rs +++ b/constructor/src/fixed_uint/core/internal/private_ops.rs @@ -471,6 +471,53 @@ impl UintConstructor { let double_unit_suffix = &self.ts.double_unit_suffix; let inner_type = &self.ts.inner_type; let part = quote!( + // # A brief explanation of division + // + // WARNING: This method is a mess, but it is fast and it works, so don't **fuck** with it. + // + // ## Problem + // + // If $x$ and $y$ is two big unsigned integer, how to compute the result of $x / y$? + // + // ## Hypothesis + // + // $$ + // p = 2^{64} \\ + // \\ + // x = \sum_{i=0}^{n} a_{i}*p^{i} \\ + // (\forall i \in [0, n), 0 \le a_{i} \lt p; 0 \lt a_{n} \lt p) \\ + // \\ + // y = \sum_{i=0}^{m} b_{i}*p^{i} \\ + // (\forall i \in [0, m), 0 \le b_{i} \lt p; 0 \lt b_{m} \lt p) \\ + // $$ + // + // ## Solution + // + // Let $z = \sum_{i=0}^{n-m} c_{i}*p^{i}$ and $r = b_{m}$. + // + // Until $n \le 2$ or $n \lt m$, do + // + // If $a_{n} > r$, set $t = a_{n}, k = n$; + // else set $t = a_{n} * p + a_{n-1}, k = n-1$. + // + // Let $s = \left \lfloor \frac{t-1}{r} \right \rfloor$. + // Let $c_{k} = c_{k} + s$. + // Let $x = x - y * s * p^{k}$ + // + // Let $l = 0$, until $x \lt y$, do + // + // $x = x - y$ + // $l = l + 1$ + // + // Let $z = z + l$ + // $z$ is the result. + // + // ## Conclusion + // + // In fact, I use a very simple and common algorithm -- [Long division]. + // The only difference is I use base-$2^{64}$ to instead of base-10. + // + // [Long division]: https://en.wikipedia.org/wiki/Long_division #[inline] fn _div_with_rem(&self, other: &Self) -> Option<(Self, Self)> { if self < other { @@ -500,8 +547,16 @@ impl UintConstructor { // lhs_idx >= ret_idx since rhs_idx >= 0 let rhs = other.inner(); let rhs_highest = rhs[rhs_idx] as #double_unit_suffix; + let mut borrow = false; loop { + if borrow { + borrow = false; + if copy.inner()[lhs_idx + 1] != 0 { + lhs_idx += 1; + ret_idx += 1; + } + } let lhs_highest = copy.inner()[lhs_idx] as #double_unit_suffix; // if lhs highest byte is ZERO, the skip it if lhs_highest == 0 { @@ -521,6 +576,7 @@ impl UintConstructor { if ret_idx == 0 { break; } + borrow = true; lhs_idx -= 1; ret_idx -= 1; (lhs_highest << #unit_bits_size) @@ -534,10 +590,11 @@ impl UintConstructor { of }; if of { - // `ret[ret_idx+1]+1` could not overflow + // `ret[ret_idx + 1] + 1` could not overflow ret[ret_idx + 1] += 1; } let minuend = { + // could not overflow let (mut minuend_tmp, _) = other._mul_unit(quotient); // left shift let mut idx = #unit_amount - 1; @@ -553,11 +610,13 @@ impl UintConstructor { } minuend_tmp }; + // could not overflow let (copy_new, _) = copy._sub(&minuend); copy = copy_new; } let mut more: #unit_suffix = 0; while copy >= *other { + // could not overflow let (copy_new, _) = copy._sub(other); copy = copy_new; more += 1; diff --git a/fixed-uint-tests/tests/regression_tests.rs b/fixed-uint-tests/tests/regression_tests.rs index 20f350e..5fd45ea 100644 --- a/fixed-uint-tests/tests/regression_tests.rs +++ b/fixed-uint-tests/tests/regression_tests.rs @@ -19,3 +19,19 @@ fn div_throw_add_overflow() { let z = ((nfuint::U256::one() << 255) / &y) << 1; assert_eq!(x, z); } + +#[test] +fn div_too_slow() { + let x = + nfuint::U4096::from_hex_str("272184cdaf3736f0fa54c1d8529a9294bcc2ac0b180838228ab").unwrap(); + let y = nfuint::U4096::from_hex_str( + "8000000000000000000000000000000000000000000000000000000000000000", + ) + .unwrap(); + let expected = nfuint::U4096::from_hex_str( + "6b851a863a5e38d58e175cb90a7b4dd5b7bcab518f09f17ade7398cc5621e239", + ) + .unwrap(); + let result = &x * &x % y; + assert_eq!(expected, result); +} From 105b8a71944db35053ff3db83b46091ad88ea9cf Mon Sep 17 00:00:00 2001 From: Boyu Yang Date: Wed, 27 Mar 2019 10:33:12 +0800 Subject: [PATCH 2/2] Fix bug: in division, when doing estimate is too cautious, so the estimated quotient is too small each time, led to an overflow. --- .../fixed_uint/core/internal/private_ops.rs | 18 ++++++++++++--- fixed-uint-tests/tests/regression_tests.rs | 22 +++++++++++++++++-- 2 files changed, 35 insertions(+), 5 deletions(-) diff --git a/constructor/src/fixed_uint/core/internal/private_ops.rs b/constructor/src/fixed_uint/core/internal/private_ops.rs index ea0b7da..fbbeabf 100644 --- a/constructor/src/fixed_uint/core/internal/private_ops.rs +++ b/constructor/src/fixed_uint/core/internal/private_ops.rs @@ -569,6 +569,7 @@ impl UintConstructor { } // below: ret_idx > 0 // estimate highest byte of quotient + // could not overflow let divisor = rhs_highest + 1; let dividend = if lhs_highest >= divisor { lhs_highest @@ -589,9 +590,20 @@ impl UintConstructor { *ret_tmp = tmp; of }; - if of { - // `ret[ret_idx + 1] + 1` could not overflow - ret[ret_idx + 1] += 1; + { + let mut ret_idx_tmp = ret_idx + 1; + let mut ret_of = of; + loop { + if ret_of { + let ret_tmp = &mut ret[ret_idx_tmp]; + let (tmp, of) = ret_tmp.overflowing_add(1); + *ret_tmp = tmp; + ret_of = of; + ret_idx_tmp += 1; + } else { + break; + } + } } let minuend = { // could not overflow diff --git a/fixed-uint-tests/tests/regression_tests.rs b/fixed-uint-tests/tests/regression_tests.rs index 5fd45ea..4fc5d96 100644 --- a/fixed-uint-tests/tests/regression_tests.rs +++ b/fixed-uint-tests/tests/regression_tests.rs @@ -7,7 +7,7 @@ // except according to those terms. #[test] -fn div_throw_add_overflow() { +fn div_throw_add_overflow_1() { let one = nfuint::U256::one(); for i in 0..255 { let x = nfuint::U256::one() << i; @@ -20,6 +20,24 @@ fn div_throw_add_overflow() { assert_eq!(x, z); } +#[test] +fn div_throw_add_overflow_2() { + let x = nfuint::U4096::from_hex_str( + "aab1deb8c8a4ba3d000000000000000000000000000000000000000000000001", + ) + .unwrap(); + let y = nfuint::U4096::from_hex_str( + "10000000000000000000000000000000000000000000000000000000000000000", + ) + .unwrap(); + let expected = nfuint::U4096::from_hex_str( + "5563bd719149747a000000000000000000000000000000000000000000000001", + ) + .unwrap(); + let result = &x * &x % y; + assert_eq!(result, expected); +} + #[test] fn div_too_slow() { let x = @@ -33,5 +51,5 @@ fn div_too_slow() { ) .unwrap(); let result = &x * &x % y; - assert_eq!(expected, result); + assert_eq!(result, expected); }