1/// Creates an unsigned division function optimized for division of integers with bitwidths
2/// larger than the largest hardware integer division supported. These functions use large radix
3/// division algorithms that require both fast division and very fast widening multiplication on the
4/// target microarchitecture. Otherwise, `impl_delegate` should be used instead.
5#[allow(unused_macros)]
6macro_rules! impl_trifecta {
7 (
8 $fn:ident, // name of the unsigned division function
9 $zero_div_fn:ident, // function called when division by zero is attempted
10 $half_division:ident, // function for division of a $uX by a $uX
11 $n_h:expr, // the number of bits in $iH or $uH
12 $uH:ident, // unsigned integer with half the bit width of $uX
13 $uX:ident, // unsigned integer with half the bit width of $uD
14 $uD:ident // unsigned integer type for the inputs and outputs of `$unsigned_name`
15 ) => {
16 /// Computes the quotient and remainder of `duo` divided by `div` and returns them as a
17 /// tuple.
18 pub fn $fn(duo: $uD, div: $uD) -> ($uD, $uD) {
19 // This is called the trifecta algorithm because it uses three main algorithms: short
20 // division for small divisors, the two possibility algorithm for large divisors, and an
21 // undersubtracting long division algorithm for intermediate cases.
22
23 // This replicates `carrying_mul` (rust-lang rfc #2417). LLVM correctly optimizes this
24 // to use a widening multiply to 128 bits on the relevant architectures.
25 fn carrying_mul(lhs: $uX, rhs: $uX) -> ($uX, $uX) {
26 let tmp = (lhs as $uD).wrapping_mul(rhs as $uD);
27 (tmp as $uX, (tmp >> ($n_h * 2)) as $uX)
28 }
29 fn carrying_mul_add(lhs: $uX, mul: $uX, add: $uX) -> ($uX, $uX) {
30 let tmp = (lhs as $uD)
31 .wrapping_mul(mul as $uD)
32 .wrapping_add(add as $uD);
33 (tmp as $uX, (tmp >> ($n_h * 2)) as $uX)
34 }
35
36 // the number of bits in a $uX
37 let n = $n_h * 2;
38
39 if div == 0 {
40 $zero_div_fn()
41 }
42
43 // Trying to use a normalization shift function will cause inelegancies in the code and
44 // inefficiencies for architectures with a native count leading zeros instruction. The
45 // undersubtracting algorithm needs both values (keeping the original `div_lz` but
46 // updating `duo_lz` multiple times), so we assume hardware support for fast
47 // `leading_zeros` calculation.
48 let div_lz = div.leading_zeros();
49 let mut duo_lz = duo.leading_zeros();
50
51 // the possible ranges of `duo` and `div` at this point:
52 // `0 <= duo < 2^n_d`
53 // `1 <= div < 2^n_d`
54
55 // quotient is 0 or 1 branch
56 if div_lz <= duo_lz {
57 // The quotient cannot be more than 1. The highest set bit of `duo` needs to be at
58 // least one place higher than `div` for the quotient to be more than 1.
59 if duo >= div {
60 return (1, duo - div);
61 } else {
62 return (0, duo);
63 }
64 }
65
66 // `_sb` is the number of significant bits (from the ones place to the highest set bit)
67 // `{2, 2^div_sb} <= duo < 2^n_d`
68 // `1 <= div < {2^duo_sb, 2^(n_d - 1)}`
69 // smaller division branch
70 if duo_lz >= n {
71 // `duo < 2^n` so it will fit in a $uX. `div` will also fit in a $uX (because of the
72 // `div_lz <= duo_lz` branch) so no numerical error.
73 let (quo, rem) = $half_division(duo as $uX, div as $uX);
74 return (quo as $uD, rem as $uD);
75 }
76
77 // `{2^n, 2^div_sb} <= duo < 2^n_d`
78 // `1 <= div < {2^duo_sb, 2^(n_d - 1)}`
79 // short division branch
80 if div_lz >= (n + $n_h) {
81 // `1 <= div < {2^duo_sb, 2^n_h}`
82
83 // It is barely possible to improve the performance of this by calculating the
84 // reciprocal and removing one `$half_division`, but only if the CPU can do fast
85 // multiplications in parallel. Other reciprocal based methods can remove two
86 // `$half_division`s, but have multiplications that cannot be done in parallel and
87 // reduce performance. I have decided to use this trivial short division method and
88 // rely on the CPU having quick divisions.
89
90 let duo_hi = (duo >> n) as $uX;
91 let div_0 = div as $uH as $uX;
92 let (quo_hi, rem_3) = $half_division(duo_hi, div_0);
93
94 let duo_mid = ((duo >> $n_h) as $uH as $uX) | (rem_3 << $n_h);
95 let (quo_1, rem_2) = $half_division(duo_mid, div_0);
96
97 let duo_lo = (duo as $uH as $uX) | (rem_2 << $n_h);
98 let (quo_0, rem_1) = $half_division(duo_lo, div_0);
99
100 return (
101 (quo_0 as $uD) | ((quo_1 as $uD) << $n_h) | ((quo_hi as $uD) << n),
102 rem_1 as $uD,
103 );
104 }
105
106 // relative leading significant bits, cannot overflow because of above branches
107 let lz_diff = div_lz - duo_lz;
108
109 // `{2^n, 2^div_sb} <= duo < 2^n_d`
110 // `2^n_h <= div < {2^duo_sb, 2^(n_d - 1)}`
111 // `mul` or `mul - 1` branch
112 if lz_diff < $n_h {
113 // Two possibility division algorithm
114
115 // The most significant bits of `duo` and `div` are within `$n_h` bits of each
116 // other. If we take the `n` most significant bits of `duo` and divide them by the
117 // corresponding bits in `div`, it produces a quotient value `quo`. It happens that
118 // `quo` or `quo - 1` will always be the correct quotient for the whole number. In
119 // other words, the bits less significant than the `n` most significant bits of
120 // `duo` and `div` can only influence the quotient to be one of two values.
121 // Because there are only two possibilities, there only needs to be one `$uH` sized
122 // division, a `$uH` by `$uD` multiplication, and only one branch with a few simple
123 // operations.
124 //
125 // Proof that the true quotient can only be `quo` or `quo - 1`.
126 // All `/` operators here are floored divisions.
127 //
128 // `shift` is the number of bits not in the higher `n` significant bits of `duo`.
129 // (definitions)
130 // 0. shift = n - duo_lz
131 // 1. duo_sig_n == duo / 2^shift
132 // 2. div_sig_n == div / 2^shift
133 // 3. quo == duo_sig_n / div_sig_n
134 //
135 //
136 // We are trying to find the true quotient, `true_quo`.
137 // 4. true_quo = duo / div. (definition)
138 //
139 // This is true because of the bits that are cut off during the bit shift.
140 // 5. duo_sig_n * 2^shift <= duo < (duo_sig_n + 1) * 2^shift.
141 // 6. div_sig_n * 2^shift <= div < (div_sig_n + 1) * 2^shift.
142 //
143 // Dividing each bound of (5) by each bound of (6) gives 4 possibilities for what
144 // `true_quo == duo / div` is bounded by:
145 // (duo_sig_n * 2^shift) / (div_sig_n * 2^shift)
146 // (duo_sig_n * 2^shift) / ((div_sig_n + 1) * 2^shift)
147 // ((duo_sig_n + 1) * 2^shift) / (div_sig_n * 2^shift)
148 // ((duo_sig_n + 1) * 2^shift) / ((div_sig_n + 1) * 2^shift)
149 //
150 // Simplifying each of these four:
151 // duo_sig_n / div_sig_n
152 // duo_sig_n / (div_sig_n + 1)
153 // (duo_sig_n + 1) / div_sig_n
154 // (duo_sig_n + 1) / (div_sig_n + 1)
155 //
156 // Taking the smallest and the largest of these as the low and high bounds
157 // and replacing `duo / div` with `true_quo`:
158 // 7. duo_sig_n / (div_sig_n + 1) <= true_quo < (duo_sig_n + 1) / div_sig_n
159 //
160 // The `lz_diff < n_h` conditional on this branch makes sure that `div_sig_n` is at
161 // least `2^n_h`, and the `div_lz <= duo_lz` branch makes sure that the highest bit
162 // of `div_sig_n` is not the `2^(n - 1)` bit.
163 // 8. `2^(n - 1) <= duo_sig_n < 2^n`
164 // 9. `2^n_h <= div_sig_n < 2^(n - 1)`
165 //
166 // We want to prove that either
167 // `(duo_sig_n + 1) / div_sig_n == duo_sig_n / (div_sig_n + 1)` or that
168 // `(duo_sig_n + 1) / div_sig_n == duo_sig_n / (div_sig_n + 1) + 1`.
169 //
170 // We also want to prove that `quo` is one of these:
171 // `duo_sig_n / div_sig_n == duo_sig_n / (div_sig_n + 1)` or
172 // `duo_sig_n / div_sig_n == (duo_sig_n + 1) / div_sig_n`.
173 //
174 // When 1 is added to the numerator of `duo_sig_n / div_sig_n` to produce
175 // `(duo_sig_n + 1) / div_sig_n`, it is not possible that the value increases by
176 // more than 1 with floored integer arithmetic and `div_sig_n != 0`. Consider
177 // `x/y + 1 < (x + 1)/y` <=> `x/y + 1 < x/y + 1/y` <=> `1 < 1/y` <=> `y < 1`.
178 // `div_sig_n` is a nonzero integer. Thus,
179 // 10. `duo_sig_n / div_sig_n == (duo_sig_n + 1) / div_sig_n` or
180 // `(duo_sig_n / div_sig_n) + 1 == (duo_sig_n + 1) / div_sig_n.
181 //
182 // When 1 is added to the denominator of `duo_sig_n / div_sig_n` to produce
183 // `duo_sig_n / (div_sig_n + 1)`, it is not possible that the value decreases by
184 // more than 1 with the bounds (8) and (9). Consider `x/y - 1 <= x/(y + 1)` <=>
185 // `(x - y)/y < x/(y + 1)` <=> `(y + 1)*(x - y) < x*y` <=> `x*y - y*y + x - y < x*y`
186 // <=> `x < y*y + y`. The smallest value of `div_sig_n` is `2^n_h` and the largest
187 // value of `duo_sig_n` is `2^n - 1`. Substituting reveals `2^n - 1 < 2^n + 2^n_h`.
188 // Thus,
189 // 11. `duo_sig_n / div_sig_n == duo_sig_n / (div_sig_n + 1)` or
190 // `(duo_sig_n / div_sig_n) - 1` == duo_sig_n / (div_sig_n + 1)`
191 //
192 // Combining both (10) and (11), we know that
193 // `quo - 1 <= duo_sig_n / (div_sig_n + 1) <= true_quo
194 // < (duo_sig_n + 1) / div_sig_n <= quo + 1` and therefore:
195 // 12. quo - 1 <= true_quo < quo + 1
196 //
197 // In a lot of division algorithms using smaller divisions to construct a larger
198 // division, we often encounter a situation where the approximate `quo` value
199 // calculated from a smaller division is multiple increments away from the true
200 // `quo` value. In those algorithms, multiple correction steps have to be applied.
201 // Those correction steps may need more multiplications to test `duo - (quo*div)`
202 // again. Because of the fact that our `quo` can only be one of two values, we can
203 // see if `duo - (quo*div)` overflows. If it did overflow, then we know that we have
204 // the larger of the two values (since the true quotient is unique, and any larger
205 // quotient will cause `duo - (quo*div)` to be negative). Also because there is only
206 // one correction needed, we can calculate the remainder `duo - (true_quo*div) ==
207 // duo - ((quo - 1)*div) == duo - (quo*div - div) == duo + div - quo*div`.
208 // If `duo - (quo*div)` did not overflow, then we have the correct answer.
209 let shift = n - duo_lz;
210 let duo_sig_n = (duo >> shift) as $uX;
211 let div_sig_n = (div >> shift) as $uX;
212 let quo = $half_division(duo_sig_n, div_sig_n).0;
213
214 // The larger `quo` value can overflow `$uD` in the right circumstances. This is a
215 // manual `carrying_mul_add` with overflow checking.
216 let div_lo = div as $uX;
217 let div_hi = (div >> n) as $uX;
218 let (tmp_lo, carry) = carrying_mul(quo, div_lo);
219 let (tmp_hi, overflow) = carrying_mul_add(quo, div_hi, carry);
220 let tmp = (tmp_lo as $uD) | ((tmp_hi as $uD) << n);
221 if (overflow != 0) || (duo < tmp) {
222 return (
223 (quo - 1) as $uD,
224 // Both the addition and subtraction can overflow, but when combined end up
225 // as a correct positive number.
226 duo.wrapping_add(div).wrapping_sub(tmp),
227 );
228 } else {
229 return (quo as $uD, duo - tmp);
230 }
231 }
232
233 // Undersubtracting long division algorithm.
234 // Instead of clearing a minimum of 1 bit from `duo` per iteration via binary long
235 // division, `n_h - 1` bits are cleared per iteration with this algorithm. It is a more
236 // complicated version of regular long division. Most integer division algorithms tend
237 // to guess a part of the quotient, and may have a larger quotient than the true
238 // quotient (which when multiplied by `div` will "oversubtract" the original dividend).
239 // They then check if the quotient was in fact too large and then have to correct it.
240 // This long division algorithm has been carefully constructed to always underguess the
241 // quotient by slim margins. This allows different subalgorithms to be blindly jumped to
242 // without needing an extra correction step.
243 //
244 // The only problem is that this subalgorithm will not work for many ranges of `duo` and
245 // `div`. Fortunately, the short division, two possibility algorithm, and other simple
246 // cases happen to exactly fill these gaps.
247 //
248 // For an example, consider the division of 76543210 by 213 and assume that `n_h` is
249 // equal to two decimal digits (note: we are working with base 10 here for readability).
250 // The first `sig_n_h` part of the divisor (21) is taken and is incremented by 1 to
251 // prevent oversubtraction. We also record the number of extra places not a part of
252 // the `sig_n` or `sig_n_h` parts.
253 //
254 // sig_n_h == 2 digits, sig_n == 4 digits
255 //
256 // vvvv <- `duo_sig_n`
257 // 76543210
258 // ^^^^ <- extra places in duo, `duo_extra == 4`
259 //
260 // vv <- `div_sig_n_h`
261 // 213
262 // ^ <- extra places in div, `div_extra == 1`
263 //
264 // The difference in extra places, `duo_extra - div_extra == extra_shl == 3`, is used
265 // for shifting partial sums in the long division.
266 //
267 // In the first step, the first `sig_n` part of duo (7654) is divided by
268 // `div_sig_n_h_add_1` (22), which results in a partial quotient of 347. This is
269 // multiplied by the whole divisor to make 73911, which is shifted left by `extra_shl`
270 // and subtracted from duo. The partial quotient is also shifted left by `extra_shl` to
271 // be added to `quo`.
272 //
273 // 347
274 // ________
275 // |76543210
276 // -73911
277 // 2632210
278 //
279 // Variables dependent on duo have to be updated:
280 //
281 // vvvv <- `duo_sig_n == 2632`
282 // 2632210
283 // ^^^ <- `duo_extra == 3`
284 //
285 // `extra_shl == 2`
286 //
287 // Two more steps are taken after this and then duo fits into `n` bits, and then a final
288 // normal long division step is made. The partial quotients are all progressively added
289 // to each other in the actual algorithm, but here I have left them all in a tower that
290 // can be added together to produce the quotient, 359357.
291 //
292 // 14
293 // 443
294 // 119
295 // 347
296 // ________
297 // |76543210
298 // -73911
299 // 2632210
300 // -25347
301 // 97510
302 // -94359
303 // 3151
304 // -2982
305 // 169 <- the remainder
306
307 let mut duo = duo;
308 let mut quo: $uD = 0;
309
310 // The number of lesser significant bits not a part of `div_sig_n_h`
311 let div_extra = (n + $n_h) - div_lz;
312
313 // The most significant `n_h` bits of div
314 let div_sig_n_h = (div >> div_extra) as $uH;
315
316 // This needs to be a `$uX` in case of overflow from the increment
317 let div_sig_n_h_add1 = (div_sig_n_h as $uX) + 1;
318
319 // `{2^n, 2^(div_sb + n_h)} <= duo < 2^n_d`
320 // `2^n_h <= div < {2^(duo_sb - n_h), 2^n}`
321 loop {
322 // The number of lesser significant bits not a part of `duo_sig_n`
323 let duo_extra = n - duo_lz;
324
325 // The most significant `n` bits of `duo`
326 let duo_sig_n = (duo >> duo_extra) as $uX;
327
328 // the two possibility algorithm requires that the difference between msbs is less
329 // than `n_h`, so the comparison is `<=` here.
330 if div_extra <= duo_extra {
331 // Undersubtracting long division step
332 let quo_part = $half_division(duo_sig_n, div_sig_n_h_add1).0 as $uD;
333 let extra_shl = duo_extra - div_extra;
334
335 // Addition to the quotient.
336 quo += (quo_part << extra_shl);
337
338 // Subtraction from `duo`. At least `n_h - 1` bits are cleared from `duo` here.
339 duo -= (div.wrapping_mul(quo_part) << extra_shl);
340 } else {
341 // Two possibility algorithm
342 let shift = n - duo_lz;
343 let duo_sig_n = (duo >> shift) as $uX;
344 let div_sig_n = (div >> shift) as $uX;
345 let quo_part = $half_division(duo_sig_n, div_sig_n).0;
346 let div_lo = div as $uX;
347 let div_hi = (div >> n) as $uX;
348
349 let (tmp_lo, carry) = carrying_mul(quo_part, div_lo);
350 // The undersubtracting long division algorithm has already run once, so
351 // overflow beyond `$uD` bits is not possible here
352 let (tmp_hi, _) = carrying_mul_add(quo_part, div_hi, carry);
353 let tmp = (tmp_lo as $uD) | ((tmp_hi as $uD) << n);
354
355 if duo < tmp {
356 return (
357 quo + ((quo_part - 1) as $uD),
358 duo.wrapping_add(div).wrapping_sub(tmp),
359 );
360 } else {
361 return (quo + (quo_part as $uD), duo - tmp);
362 }
363 }
364
365 duo_lz = duo.leading_zeros();
366
367 if div_lz <= duo_lz {
368 // quotient can have 0 or 1 added to it
369 if div <= duo {
370 return (quo + 1, duo - div);
371 } else {
372 return (quo, duo);
373 }
374 }
375
376 // This can only happen if `div_sd < n` (because of previous "quo = 0 or 1"
377 // branches), but it is not worth it to unroll further.
378 if n <= duo_lz {
379 // simple division and addition
380 let tmp = $half_division(duo as $uX, div as $uX);
381 return (quo + (tmp.0 as $uD), tmp.1 as $uD);
382 }
383 }
384 }
385 };
386}
387