/// Creates an unsigned division function optimized for division of integers with bitwidths
/// larger than the largest hardware integer division supported. These functions use large radix
/// division algorithms that require both fast division and very fast widening multiplication on the
/// target microarchitecture. Otherwise, `impl_delegate` should be used instead.
#[allow(unused_macros)]
macro_rules! impl_trifecta {
(
$fn:ident, // name of the unsigned division function
$zero_div_fn:ident, // function called when division by zero is attempted
$half_division:ident, // function for division of a $uX by a $uX
$n_h:expr, // the number of bits in $iH or $uH
$uH:ident, // unsigned integer with half the bit width of $uX
$uX:ident, // unsigned integer with half the bit width of $uD
$uD:ident // unsigned integer type for the inputs and outputs of `$unsigned_name`
) => {
/// Computes the quotient and remainder of `duo` divided by `div` and returns them as a
/// tuple.
pub fn $fn(duo: $uD, div: $uD) -> ($uD, $uD) {
// This is called the trifecta algorithm because it uses three main algorithms: short
// division for small divisors, the two possibility algorithm for large divisors, and an
// undersubtracting long division algorithm for intermediate cases.
// This replicates `carrying_mul` (rust-lang rfc #2417). LLVM correctly optimizes this
// to use a widening multiply to 128 bits on the relevant architectures.
fn carrying_mul(lhs: $uX, rhs: $uX) -> ($uX, $uX) {
let tmp = (lhs as $uD).wrapping_mul(rhs as $uD);
(tmp as $uX, (tmp >> ($n_h * 2)) as $uX)
}
fn carrying_mul_add(lhs: $uX, mul: $uX, add: $uX) -> ($uX, $uX) {
let tmp = (lhs as $uD)
.wrapping_mul(mul as $uD)
.wrapping_add(add as $uD);
(tmp as $uX, (tmp >> ($n_h * 2)) as $uX)
}
// the number of bits in a $uX
let n = $n_h * 2;
if div == 0 {
$zero_div_fn()
}
// Trying to use a normalization shift function will cause inelegancies in the code and
// inefficiencies for architectures with a native count leading zeros instruction. The
// undersubtracting algorithm needs both values (keeping the original `div_lz` but
// updating `duo_lz` multiple times), so we assume hardware support for fast
// `leading_zeros` calculation.
let div_lz = div.leading_zeros();
let mut duo_lz = duo.leading_zeros();
// the possible ranges of `duo` and `div` at this point:
// `0 <= duo < 2^n_d`
// `1 <= div < 2^n_d`
// quotient is 0 or 1 branch
if div_lz <= duo_lz {
// The quotient cannot be more than 1. The highest set bit of `duo` needs to be at
// least one place higher than `div` for the quotient to be more than 1.
if duo >= div {
return (1, duo - div);
} else {
return (0, duo);
}
}
// `_sb` is the number of significant bits (from the ones place to the highest set bit)
// `{2, 2^div_sb} <= duo < 2^n_d`
// `1 <= div < {2^duo_sb, 2^(n_d - 1)}`
// smaller division branch
if duo_lz >= n {
// `duo < 2^n` so it will fit in a $uX. `div` will also fit in a $uX (because of the
// `div_lz <= duo_lz` branch) so no numerical error.
let (quo, rem) = $half_division(duo as $uX, div as $uX);
return (quo as $uD, rem as $uD);
}
// `{2^n, 2^div_sb} <= duo < 2^n_d`
// `1 <= div < {2^duo_sb, 2^(n_d - 1)}`
// short division branch
if div_lz >= (n + $n_h) {
// `1 <= div < {2^duo_sb, 2^n_h}`
// It is barely possible to improve the performance of this by calculating the
// reciprocal and removing one `$half_division`, but only if the CPU can do fast
// multiplications in parallel. Other reciprocal based methods can remove two
// `$half_division`s, but have multiplications that cannot be done in parallel and
// reduce performance. I have decided to use this trivial short division method and
// rely on the CPU having quick divisions.
let duo_hi = (duo >> n) as $uX;
let div_0 = div as $uH as $uX;
let (quo_hi, rem_3) = $half_division(duo_hi, div_0);
let duo_mid = ((duo >> $n_h) as $uH as $uX) | (rem_3 << $n_h);
let (quo_1, rem_2) = $half_division(duo_mid, div_0);
let duo_lo = (duo as $uH as $uX) | (rem_2 << $n_h);
let (quo_0, rem_1) = $half_division(duo_lo, div_0);
return (
(quo_0 as $uD) | ((quo_1 as $uD) << $n_h) | ((quo_hi as $uD) << n),
rem_1 as $uD,
);
}
// relative leading significant bits, cannot overflow because of above branches
let lz_diff = div_lz - duo_lz;
// `{2^n, 2^div_sb} <= duo < 2^n_d`
// `2^n_h <= div < {2^duo_sb, 2^(n_d - 1)}`
// `mul` or `mul - 1` branch
if lz_diff < $n_h {
// Two possibility division algorithm
// The most significant bits of `duo` and `div` are within `$n_h` bits of each
// other. If we take the `n` most significant bits of `duo` and divide them by the
// corresponding bits in `div`, it produces a quotient value `quo`. It happens that
// `quo` or `quo - 1` will always be the correct quotient for the whole number. In
// other words, the bits less significant than the `n` most significant bits of
// `duo` and `div` can only influence the quotient to be one of two values.
// Because there are only two possibilities, there only needs to be one `$uH` sized
// division, a `$uH` by `$uD` multiplication, and only one branch with a few simple
// operations.
//
// Proof that the true quotient can only be `quo` or `quo - 1`.
// All `/` operators here are floored divisions.
//
// `shift` is the number of bits not in the higher `n` significant bits of `duo`.
// (definitions)
// 0. shift = n - duo_lz
// 1. duo_sig_n == duo / 2^shift
// 2. div_sig_n == div / 2^shift
// 3. quo == duo_sig_n / div_sig_n
//
//
// We are trying to find the true quotient, `true_quo`.
// 4. true_quo = duo / div. (definition)
//
// This is true because of the bits that are cut off during the bit shift.
// 5. duo_sig_n * 2^shift <= duo < (duo_sig_n + 1) * 2^shift.
// 6. div_sig_n * 2^shift <= div < (div_sig_n + 1) * 2^shift.
//
// Dividing each bound of (5) by each bound of (6) gives 4 possibilities for what
// `true_quo == duo / div` is bounded by:
// (duo_sig_n * 2^shift) / (div_sig_n * 2^shift)
// (duo_sig_n * 2^shift) / ((div_sig_n + 1) * 2^shift)
// ((duo_sig_n + 1) * 2^shift) / (div_sig_n * 2^shift)
// ((duo_sig_n + 1) * 2^shift) / ((div_sig_n + 1) * 2^shift)
//
// Simplifying each of these four:
// duo_sig_n / div_sig_n
// duo_sig_n / (div_sig_n + 1)
// (duo_sig_n + 1) / div_sig_n
// (duo_sig_n + 1) / (div_sig_n + 1)
//
// Taking the smallest and the largest of these as the low and high bounds
// and replacing `duo / div` with `true_quo`:
// 7. duo_sig_n / (div_sig_n + 1) <= true_quo < (duo_sig_n + 1) / div_sig_n
//
// The `lz_diff < n_h` conditional on this branch makes sure that `div_sig_n` is at
// least `2^n_h`, and the `div_lz <= duo_lz` branch makes sure that the highest bit
// of `div_sig_n` is not the `2^(n - 1)` bit.
// 8. `2^(n - 1) <= duo_sig_n < 2^n`
// 9. `2^n_h <= div_sig_n < 2^(n - 1)`
//
// We want to prove that either
// `(duo_sig_n + 1) / div_sig_n == duo_sig_n / (div_sig_n + 1)` or that
// `(duo_sig_n + 1) / div_sig_n == duo_sig_n / (div_sig_n + 1) + 1`.
//
// We also want to prove that `quo` is one of these:
// `duo_sig_n / div_sig_n == duo_sig_n / (div_sig_n + 1)` or
// `duo_sig_n / div_sig_n == (duo_sig_n + 1) / div_sig_n`.
//
// When 1 is added to the numerator of `duo_sig_n / div_sig_n` to produce
// `(duo_sig_n + 1) / div_sig_n`, it is not possible that the value increases by
// more than 1 with floored integer arithmetic and `div_sig_n != 0`. Consider
// `x/y + 1 < (x + 1)/y` <=> `x/y + 1 < x/y + 1/y` <=> `1 < 1/y` <=> `y < 1`.
// `div_sig_n` is a nonzero integer. Thus,
// 10. `duo_sig_n / div_sig_n == (duo_sig_n + 1) / div_sig_n` or
// `(duo_sig_n / div_sig_n) + 1 == (duo_sig_n + 1) / div_sig_n.
//
// When 1 is added to the denominator of `duo_sig_n / div_sig_n` to produce
// `duo_sig_n / (div_sig_n + 1)`, it is not possible that the value decreases by
// more than 1 with the bounds (8) and (9). Consider `x/y - 1 <= x/(y + 1)` <=>
// `(x - y)/y < x/(y + 1)` <=> `(y + 1)*(x - y) < x*y` <=> `x*y - y*y + x - y < x*y`
// <=> `x < y*y + y`. The smallest value of `div_sig_n` is `2^n_h` and the largest
// value of `duo_sig_n` is `2^n - 1`. Substituting reveals `2^n - 1 < 2^n + 2^n_h`.
// Thus,
// 11. `duo_sig_n / div_sig_n == duo_sig_n / (div_sig_n + 1)` or
// `(duo_sig_n / div_sig_n) - 1` == duo_sig_n / (div_sig_n + 1)`
//
// Combining both (10) and (11), we know that
// `quo - 1 <= duo_sig_n / (div_sig_n + 1) <= true_quo
// < (duo_sig_n + 1) / div_sig_n <= quo + 1` and therefore:
// 12. quo - 1 <= true_quo < quo + 1
//
// In a lot of division algorithms using smaller divisions to construct a larger
// division, we often encounter a situation where the approximate `quo` value
// calculated from a smaller division is multiple increments away from the true
// `quo` value. In those algorithms, multiple correction steps have to be applied.
// Those correction steps may need more multiplications to test `duo - (quo*div)`
// again. Because of the fact that our `quo` can only be one of two values, we can
// see if `duo - (quo*div)` overflows. If it did overflow, then we know that we have
// the larger of the two values (since the true quotient is unique, and any larger
// quotient will cause `duo - (quo*div)` to be negative). Also because there is only
// one correction needed, we can calculate the remainder `duo - (true_quo*div) ==
// duo - ((quo - 1)*div) == duo - (quo*div - div) == duo + div - quo*div`.
// If `duo - (quo*div)` did not overflow, then we have the correct answer.
let shift = n - duo_lz;
let duo_sig_n = (duo >> shift) as $uX;
let div_sig_n = (div >> shift) as $uX;
let quo = $half_division(duo_sig_n, div_sig_n).0;
// The larger `quo` value can overflow `$uD` in the right circumstances. This is a
// manual `carrying_mul_add` with overflow checking.
let div_lo = div as $uX;
let div_hi = (div >> n) as $uX;
let (tmp_lo, carry) = carrying_mul(quo, div_lo);
let (tmp_hi, overflow) = carrying_mul_add(quo, div_hi, carry);
let tmp = (tmp_lo as $uD) | ((tmp_hi as $uD) << n);
if (overflow != 0) || (duo < tmp) {
return (
(quo - 1) as $uD,
// Both the addition and subtraction can overflow, but when combined end up
// as a correct positive number.
duo.wrapping_add(div).wrapping_sub(tmp),
);
} else {
return (quo as $uD, duo - tmp);
}
}
// Undersubtracting long division algorithm.
// Instead of clearing a minimum of 1 bit from `duo` per iteration via binary long
// division, `n_h - 1` bits are cleared per iteration with this algorithm. It is a more
// complicated version of regular long division. Most integer division algorithms tend
// to guess a part of the quotient, and may have a larger quotient than the true
// quotient (which when multiplied by `div` will "oversubtract" the original dividend).
// They then check if the quotient was in fact too large and then have to correct it.
// This long division algorithm has been carefully constructed to always underguess the
// quotient by slim margins. This allows different subalgorithms to be blindly jumped to
// without needing an extra correction step.
//
// The only problem is that this subalgorithm will not work for many ranges of `duo` and
// `div`. Fortunately, the short division, two possibility algorithm, and other simple
// cases happen to exactly fill these gaps.
//
// For an example, consider the division of 76543210 by 213 and assume that `n_h` is
// equal to two decimal digits (note: we are working with base 10 here for readability).
// The first `sig_n_h` part of the divisor (21) is taken and is incremented by 1 to
// prevent oversubtraction. We also record the number of extra places not a part of
// the `sig_n` or `sig_n_h` parts.
//
// sig_n_h == 2 digits, sig_n == 4 digits
//
// vvvv <- `duo_sig_n`
// 76543210
// ^^^^ <- extra places in duo, `duo_extra == 4`
//
// vv <- `div_sig_n_h`
// 213
// ^ <- extra places in div, `div_extra == 1`
//
// The difference in extra places, `duo_extra - div_extra == extra_shl == 3`, is used
// for shifting partial sums in the long division.
//
// In the first step, the first `sig_n` part of duo (7654) is divided by
// `div_sig_n_h_add_1` (22), which results in a partial quotient of 347. This is
// multiplied by the whole divisor to make 73911, which is shifted left by `extra_shl`
// and subtracted from duo. The partial quotient is also shifted left by `extra_shl` to
// be added to `quo`.
//
// 347
// ________
// |76543210
// -73911
// 2632210
//
// Variables dependent on duo have to be updated:
//
// vvvv <- `duo_sig_n == 2632`
// 2632210
// ^^^ <- `duo_extra == 3`
//
// `extra_shl == 2`
//
// Two more steps are taken after this and then duo fits into `n` bits, and then a final
// normal long division step is made. The partial quotients are all progressively added
// to each other in the actual algorithm, but here I have left them all in a tower that
// can be added together to produce the quotient, 359357.
//
// 14
// 443
// 119
// 347
// ________
// |76543210
// -73911
// 2632210
// -25347
// 97510
// -94359
// 3151
// -2982
// 169 <- the remainder
let mut duo = duo;
let mut quo: $uD = 0;
// The number of lesser significant bits not a part of `div_sig_n_h`
let div_extra = (n + $n_h) - div_lz;
// The most significant `n_h` bits of div
let div_sig_n_h = (div >> div_extra) as $uH;
// This needs to be a `$uX` in case of overflow from the increment
let div_sig_n_h_add1 = (div_sig_n_h as $uX) + 1;
// `{2^n, 2^(div_sb + n_h)} <= duo < 2^n_d`
// `2^n_h <= div < {2^(duo_sb - n_h), 2^n}`
loop {
// The number of lesser significant bits not a part of `duo_sig_n`
let duo_extra = n - duo_lz;
// The most significant `n` bits of `duo`
let duo_sig_n = (duo >> duo_extra) as $uX;
// the two possibility algorithm requires that the difference between msbs is less
// than `n_h`, so the comparison is `<=` here.
if div_extra <= duo_extra {
// Undersubtracting long division step
let quo_part = $half_division(duo_sig_n, div_sig_n_h_add1).0 as $uD;
let extra_shl = duo_extra - div_extra;
// Addition to the quotient.
quo += (quo_part << extra_shl);
// Subtraction from `duo`. At least `n_h - 1` bits are cleared from `duo` here.
duo -= (div.wrapping_mul(quo_part) << extra_shl);
} else {
// Two possibility algorithm
let shift = n - duo_lz;
let duo_sig_n = (duo >> shift) as $uX;
let div_sig_n = (div >> shift) as $uX;
let quo_part = $half_division(duo_sig_n, div_sig_n).0;
let div_lo = div as $uX;
let div_hi = (div >> n) as $uX;
let (tmp_lo, carry) = carrying_mul(quo_part, div_lo);
// The undersubtracting long division algorithm has already run once, so
// overflow beyond `$uD` bits is not possible here
let (tmp_hi, _) = carrying_mul_add(quo_part, div_hi, carry);
let tmp = (tmp_lo as $uD) | ((tmp_hi as $uD) << n);
if duo < tmp {
return (
quo + ((quo_part - 1) as $uD),
duo.wrapping_add(div).wrapping_sub(tmp),
);
} else {
return (quo + (quo_part as $uD), duo - tmp);
}
}
duo_lz = duo.leading_zeros();
if div_lz <= duo_lz {
// quotient can have 0 or 1 added to it
if div <= duo {
return (quo + 1, duo - div);
} else {
return (quo, duo);
}
}
// This can only happen if `div_sd < n` (because of previous "quo = 0 or 1"
// branches), but it is not worth it to unroll further.
if n <= duo_lz {
// simple division and addition
let tmp = $half_division(duo as $uX, div as $uX);
return (quo + (tmp.0 as $uD), tmp.1 as $uD);
}
}
}
};
}