diff --git a/AK/BigIntBase.h b/AK/BigIntBase.h index 9dc91333089..61a5e2da21b 100644 --- a/AK/BigIntBase.h +++ b/AK/BigIntBase.h @@ -147,6 +147,8 @@ requires(bit_size <= max_big_int_length * native_word_size) struct StaticStorage { return m_data; } + + constexpr operator StorageSpan() { return { m_data, static_size }; } }; struct IntegerWrapper { @@ -267,6 +269,12 @@ ALWAYS_INLINE constexpr WordType sub_words(WordType word1, WordType word2, bool& return output; } +template +ALWAYS_INLINE constexpr DoubleWord wide_multiply(WordType word1, WordType word2) +{ + return static_cast>(word1) * word2; +} + template constexpr DoubleWord dword(WordType low, WordType high) { @@ -584,6 +592,74 @@ struct StorageOperations { if (size2 < size && (sign1 ^ sign2)) negate(result, result); } + + template + static constexpr void div_mod_internal( + StorageSpan dividend, StorageSpan divisor, + StorageSpan quotient, StorageSpan remainder, + size_t dividend_len, size_t divisor_len) + { + // Knuth's algorithm D + // D1. Normalize + // FIXME: Investigate GCC producing bogus -Warray-bounds when dividing u128 by u32. This code + // should not be reachable at all in this case because fast paths above cover all cases + // when `operand2.size() == 1`. + AK_IGNORE_DIAGNOSTIC("-Warray-bounds", size_t shift = count_leading_zeroes(divisor[divisor_len - 1]);) + shift_left(dividend, shift, dividend); + shift_left(divisor, shift, divisor); + + auto divisor_approx = divisor[divisor_len - 1]; + + for (size_t i = dividend_len + 1; i-- > divisor_len;) { + // D3. Calculate qhat + WordType qhat; + VERIFY(dividend[i] <= divisor_approx); + if (dividend[i] == divisor_approx) { + qhat = NumericLimits::max(); + } else { + WordType rhat; + qhat = div_mod_words(dividend[i - 1], dividend[i], divisor_approx, rhat); + + auto is_qhat_too_large = [&] { + return wide_multiply(qhat, divisor[divisor_len - 2]) > dword(dividend[i - 2], rhat); + }; + if (is_qhat_too_large()) { + --qhat; + bool carry = false; + rhat = add_words(rhat, divisor_approx, carry); + if (!carry && is_qhat_too_large()) + --qhat; + } + } + + // D4. Multiply & subtract + WordType mul_carry = 0; + bool sub_carry = false; + for (size_t j = 0; j < divisor_len; ++j) { + auto mul_result = wide_multiply(qhat, divisor[j]) + mul_carry; + auto& output = dividend[i + j - divisor_len]; + output = sub_words(output, static_cast(mul_result), sub_carry); + mul_carry = mul_result >> word_size; + } + dividend[i] = sub_words(dividend[i], mul_carry, sub_carry); + + if (sub_carry) { + // D6. Add back + auto dividend_part = StorageSpan { dividend.slice(i - divisor_len, divisor_len + 1) }; + auto overflow = add(dividend_part, divisor, dividend_part); + VERIFY(overflow == 1); + } + + quotient[i - divisor_len] = qhat - sub_carry; + } + + for (size_t i = dividend_len - divisor_len + 1; i < quotient.size(); ++i) + quotient[i] = 0; + + // D8. Unnormalize + if constexpr (restore_remainder) + shift_right(StorageSpan { dividend.trim(remainder.size()) }, shift, remainder); + } }; } diff --git a/AK/UFixedBigIntDivision.h b/AK/UFixedBigIntDivision.h index 398a93ffbc0..64090d8a78b 100644 --- a/AK/UFixedBigIntDivision.h +++ b/AK/UFixedBigIntDivision.h @@ -74,64 +74,7 @@ constexpr void div_mod_internal( Ops::copy(operand1, dividend); auto divisor = operand2; - // D1. Normalize - // FIXME: Investigate GCC producing bogus -Warray-bounds when dividing u128 by u32. This code - // should not be reachable at all in this case because fast paths above cover all cases - // when `operand2.size() == 1`. - AK_IGNORE_DIAGNOSTIC("-Warray-bounds", size_t shift = count_leading_zeroes(divisor[divisor_len - 1]);) - Ops::shift_left(dividend, shift, dividend); - Ops::shift_left(divisor, shift, divisor); - - auto divisor_approx = divisor[divisor_len - 1]; - - for (size_t i = dividend_len + 1; i-- > divisor_len;) { - // D3. Calculate qhat - NativeWord qhat; - VERIFY(dividend[i] <= divisor_approx); - if (dividend[i] == divisor_approx) { - qhat = max_native_word; - } else { - NativeWord rhat; - qhat = div_mod_words(dividend[i - 1], dividend[i], divisor_approx, rhat); - - auto is_qhat_too_large = [&] { - return UFixedBigInt { qhat }.wide_multiply(divisor[divisor_len - 2]) > UFixedBigInt { dividend[i - 2], rhat }; - }; - if (is_qhat_too_large()) { - --qhat; - bool carry = false; - rhat = add_words(rhat, divisor_approx, carry); - if (!carry && is_qhat_too_large()) - --qhat; - } - } - - // D4. Multiply & subtract - NativeWord mul_carry = 0; - bool sub_carry = false; - for (size_t j = 0; j < divisor_len; ++j) { - auto mul_result = UFixedBigInt { qhat }.wide_multiply(divisor[j]) + mul_carry; - auto& output = dividend[i + j - divisor_len]; - output = sub_words(output, mul_result.low(), sub_carry); - mul_carry = mul_result.high(); - } - dividend[i] = sub_words(dividend[i], mul_carry, sub_carry); - - if (sub_carry) { - // D6. Add back - auto dividend_part = UnsignedStorageSpan { dividend.data() + i - divisor_len, divisor_len + 1 }; - VERIFY(Ops::add(dividend_part, divisor, dividend_part)); - } - - quotient[i - divisor_len] = qhat - sub_carry; - } - - for (size_t i = dividend_len - divisor_len + 1; i < quotient.size(); ++i) - quotient[i] = 0; - - // D8. Unnormalize - if constexpr (restore_remainder) - Ops::shift_right(UnsignedStorageSpan { dividend.data(), remainder.size() }, shift, remainder); + Ops::div_mod_internal(dividend, divisor, quotient, remainder, dividend_len, divisor_len); } }