1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
|
/*
* Copyright (c) 2023, Dan Klishch <danilklishch@gmail.com>
*
* SPDX-License-Identifier: BSD-2-Clause
*/
#pragma once
#include <AK/Diagnostics.h>
#include <AK/UFixedBigInt.h>
namespace AK {
namespace Detail {
template<size_t dividend_bit_size, size_t divisor_bit_size, bool restore_remainder>
constexpr void div_mod_internal(
StaticStorage<false, dividend_bit_size> const& operand1,
StaticStorage<false, divisor_bit_size> const& operand2,
StaticStorage<false, dividend_bit_size>& quotient,
StaticStorage<false, divisor_bit_size>& remainder)
{
size_t dividend_len = operand1.size(), divisor_len = operand2.size();
while (divisor_len > 0 && !operand2[divisor_len - 1])
--divisor_len;
while (dividend_len > 0 && !operand1[dividend_len - 1])
--dividend_len;
// FIXME: Should raise SIGFPE instead
VERIFY(divisor_len); // VERIFY(divisor != 0)
// Fast paths
if (divisor_len == 1 && operand2[0] == 1) { // divisor == 1
quotient = operand1;
if constexpr (restore_remainder)
StorageOperations::set(0, remainder);
return;
}
if (dividend_len < divisor_len) { // dividend < divisor
StorageOperations::set(0, quotient);
if constexpr (restore_remainder)
remainder = operand1;
return;
}
if (divisor_len == 1 && dividend_len == 1) { // NativeWord / NativeWord
StorageOperations::set(operand1[0] / operand2[0], quotient);
if constexpr (restore_remainder)
StorageOperations::set(operand1[0] % operand2[0], remainder);
return;
}
if (divisor_len == 1) { // BigInt by NativeWord
auto u = (static_cast<DoubleWord>(operand1[dividend_len - 1]) << word_size) + operand1[dividend_len - 2];
auto divisor = operand2[0];
auto top = u / divisor;
quotient[dividend_len - 1] = static_cast<NativeWord>(top >> word_size);
quotient[dividend_len - 2] = static_cast<NativeWord>(top);
auto carry = static_cast<NativeWord>(u % divisor);
for (size_t i = dividend_len - 2; i--;)
quotient[i] = div_mod_words(operand1[i], carry, divisor, carry);
for (size_t i = dividend_len; i < quotient.size(); ++i)
quotient[i] = 0;
if constexpr (restore_remainder)
StorageOperations::set(carry, remainder);
return;
}
// Knuth's algorithm D
StaticStorage<false, dividend_bit_size + word_size> dividend;
StorageOperations::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]);)
StorageOperations::shift_left(dividend, shift, dividend);
StorageOperations::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_word;
} else {
NativeWord rhat;
qhat = div_mod_words(dividend[i - 1], dividend[i], divisor_approx, rhat);
auto is_qhat_too_large = [&] {
return UFixedBigInt<word_size> { qhat }.wide_multiply(divisor[divisor_len - 2]) > u128 { 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<word_size> { 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(StorageOperations::add<false>(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)
StorageOperations::shift_right(UnsignedStorageSpan { dividend.data(), remainder.size() }, shift, remainder);
}
}
}
|
