/// 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);
                }
            }
        }
    };
}