1 // Copyright 2021 the V8 project authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style license that can be
3 // found in the LICENSE file.
4
5 // Karatsuba multiplication. This is loosely based on Go's implementation
6 // found at https://golang.org/src/math/big/nat.go, licensed as follows:
7 //
8 // Copyright 2009 The Go Authors. All rights reserved.
9 // Use of this source code is governed by a BSD-style
10 // license that can be found in the LICENSE file [1].
11 //
12 // [1] https://golang.org/LICENSE
13
14 #include <algorithm>
15 #include <utility>
16
17 #include "src/bigint/bigint-internal.h"
18 #include "src/bigint/digit-arithmetic.h"
19 #include "src/bigint/util.h"
20 #include "src/bigint/vector-arithmetic.h"
21
22 namespace v8 {
23 namespace bigint {
24
25 // If Karatsuba is the best supported algorithm, then it must check for
26 // termination requests. If there are more advanced algorithms available
27 // for larger inputs, then Karatsuba will only be used for sufficiently
28 // small chunks that checking for termination requests is not necessary.
29 #if V8_ADVANCED_BIGINT_ALGORITHMS
30 #define MAYBE_TERMINATE
31 #else
32 #define MAYBE_TERMINATE \
33 if (should_terminate()) return;
34 #endif
35
36 namespace {
37
38 // The Karatsuba algorithm sometimes finishes more quickly when the
39 // input length is rounded up a bit. This method encodes some heuristics
40 // to accomplish this. The details have been determined experimentally.
RoundUpLen(int len)41 int RoundUpLen(int len) {
42 if (len <= 36) return RoundUp(len, 2);
43 // Keep the 4 or 5 most significant non-zero bits.
44 int shift = BitLength(len) - 5;
45 if ((len >> shift) >= 0x18) {
46 shift++;
47 }
48 // Round up, unless we're only just above the threshold. This smoothes
49 // the steps by which time goes up as input size increases.
50 int additive = ((1 << shift) - 1);
51 if (shift >= 2 && (len & additive) < (1 << (shift - 2))) {
52 return len;
53 }
54 return ((len + additive) >> shift) << shift;
55 }
56
57 // This method makes the final decision how much to bump up the input size.
KaratsubaLength(int n)58 int KaratsubaLength(int n) {
59 n = RoundUpLen(n);
60 int i = 0;
61 while (n > kKaratsubaThreshold) {
62 n >>= 1;
63 i++;
64 }
65 return n << i;
66 }
67
68 // Performs the specific subtraction required by {KaratsubaMain} below.
KaratsubaSubtractionHelper(RWDigits result,Digits X,Digits Y,int * sign)69 void KaratsubaSubtractionHelper(RWDigits result, Digits X, Digits Y,
70 int* sign) {
71 X.Normalize();
72 Y.Normalize();
73 digit_t borrow = 0;
74 int i = 0;
75 if (!GreaterThanOrEqual(X, Y)) {
76 *sign = -(*sign);
77 std::swap(X, Y);
78 }
79 for (; i < Y.len(); i++) {
80 result[i] = digit_sub2(X[i], Y[i], borrow, &borrow);
81 }
82 for (; i < X.len(); i++) {
83 result[i] = digit_sub(X[i], borrow, &borrow);
84 }
85 DCHECK(borrow == 0); // NOLINT(readability/check)
86 for (; i < result.len(); i++) result[i] = 0;
87 }
88
89 } // namespace
90
MultiplyKaratsuba(RWDigits Z,Digits X,Digits Y)91 void ProcessorImpl::MultiplyKaratsuba(RWDigits Z, Digits X, Digits Y) {
92 DCHECK(X.len() >= Y.len());
93 DCHECK(Y.len() >= kKaratsubaThreshold);
94 DCHECK(Z.len() >= X.len() + Y.len());
95 int k = KaratsubaLength(Y.len());
96 int scratch_len = 4 * k;
97 ScratchDigits scratch(scratch_len);
98 KaratsubaStart(Z, X, Y, scratch, k);
99 }
100
101 // Entry point for Karatsuba-based multiplication, takes care of inputs
102 // with unequal lengths by chopping the larger into chunks.
KaratsubaStart(RWDigits Z,Digits X,Digits Y,RWDigits scratch,int k)103 void ProcessorImpl::KaratsubaStart(RWDigits Z, Digits X, Digits Y,
104 RWDigits scratch, int k) {
105 KaratsubaMain(Z, X, Y, scratch, k);
106 MAYBE_TERMINATE
107 for (int i = 2 * k; i < Z.len(); i++) Z[i] = 0;
108 if (k < Y.len() || X.len() != Y.len()) {
109 ScratchDigits T(2 * k);
110 // Add X0 * Y1 * b.
111 Digits X0(X, 0, k);
112 Digits Y1 = Y + std::min(k, Y.len());
113 if (Y1.len() > 0) {
114 KaratsubaChunk(T, X0, Y1, scratch);
115 MAYBE_TERMINATE
116 AddAndReturnOverflow(Z + k, T); // Can't overflow.
117 }
118
119 // Add Xi * Y0 << i and Xi * Y1 * b << (i + k).
120 Digits Y0(Y, 0, k);
121 for (int i = k; i < X.len(); i += k) {
122 Digits Xi(X, i, k);
123 KaratsubaChunk(T, Xi, Y0, scratch);
124 MAYBE_TERMINATE
125 AddAndReturnOverflow(Z + i, T); // Can't overflow.
126 if (Y1.len() > 0) {
127 KaratsubaChunk(T, Xi, Y1, scratch);
128 MAYBE_TERMINATE
129 AddAndReturnOverflow(Z + (i + k), T); // Can't overflow.
130 }
131 }
132 }
133 }
134
135 // Entry point for chunk-wise multiplications, selects an appropriate
136 // algorithm for the inputs based on their sizes.
KaratsubaChunk(RWDigits Z,Digits X,Digits Y,RWDigits scratch)137 void ProcessorImpl::KaratsubaChunk(RWDigits Z, Digits X, Digits Y,
138 RWDigits scratch) {
139 X.Normalize();
140 Y.Normalize();
141 if (X.len() == 0 || Y.len() == 0) return Z.Clear();
142 if (X.len() < Y.len()) std::swap(X, Y);
143 if (Y.len() == 1) return MultiplySingle(Z, X, Y[0]);
144 if (Y.len() < kKaratsubaThreshold) return MultiplySchoolbook(Z, X, Y);
145 int k = KaratsubaLength(Y.len());
146 DCHECK(scratch.len() >= 4 * k);
147 return KaratsubaStart(Z, X, Y, scratch, k);
148 }
149
150 // The main recursive Karatsuba method.
KaratsubaMain(RWDigits Z,Digits X,Digits Y,RWDigits scratch,int n)151 void ProcessorImpl::KaratsubaMain(RWDigits Z, Digits X, Digits Y,
152 RWDigits scratch, int n) {
153 if (n < kKaratsubaThreshold) {
154 X.Normalize();
155 Y.Normalize();
156 if (X.len() >= Y.len()) {
157 return MultiplySchoolbook(RWDigits(Z, 0, 2 * n), X, Y);
158 } else {
159 return MultiplySchoolbook(RWDigits(Z, 0, 2 * n), Y, X);
160 }
161 }
162 DCHECK(scratch.len() >= 4 * n);
163 DCHECK((n & 1) == 0); // NOLINT(readability/check)
164 int n2 = n >> 1;
165 Digits X0(X, 0, n2);
166 Digits X1(X, n2, n2);
167 Digits Y0(Y, 0, n2);
168 Digits Y1(Y, n2, n2);
169 RWDigits scratch_for_recursion(scratch, 2 * n, 2 * n);
170 RWDigits P0(scratch, 0, n);
171 KaratsubaMain(P0, X0, Y0, scratch_for_recursion, n2);
172 MAYBE_TERMINATE
173 for (int i = 0; i < n; i++) Z[i] = P0[i];
174 RWDigits P2(scratch, n, n);
175 KaratsubaMain(P2, X1, Y1, scratch_for_recursion, n2);
176 MAYBE_TERMINATE
177 RWDigits Z2 = Z + n;
178 int end = std::min(Z2.len(), P2.len());
179 for (int i = 0; i < end; i++) Z2[i] = P2[i];
180 for (int i = end; i < n; i++) {
181 DCHECK(P2[i] == 0); // NOLINT(readability/check)
182 }
183 // The intermediate result can be one digit too large; the subtraction
184 // below will fix this.
185 digit_t overflow = AddAndReturnOverflow(Z + n2, P0);
186 overflow += AddAndReturnOverflow(Z + n2, P2);
187 RWDigits X_diff(scratch, 0, n2);
188 RWDigits Y_diff(scratch, n2, n2);
189 int sign = 1;
190 KaratsubaSubtractionHelper(X_diff, X1, X0, &sign);
191 KaratsubaSubtractionHelper(Y_diff, Y0, Y1, &sign);
192 RWDigits P1(scratch, n, n);
193 KaratsubaMain(P1, X_diff, Y_diff, scratch_for_recursion, n2);
194 if (sign > 0) {
195 overflow += AddAndReturnOverflow(Z + n2, P1);
196 } else {
197 overflow -= SubAndReturnBorrow(Z + n2, P1);
198 }
199 // The intermediate result may have been bigger, but the final result fits.
200 DCHECK(overflow == 0); // NOLINT(readability/check)
201 USE(overflow);
202 }
203
204 #undef MAYBE_TERMINATE
205
206 } // namespace bigint
207 } // namespace v8
208