divmod

fn divmod(
    q: *Mutable,
    r: *Mutable,
    x: *Mutable,
    y: *Mutable,
) void {
    // 0.
    // Normalize so that y[t] > b/2
    const lz = @clz(y.limbs[y.len - 1]);
    const norm_shift = if (lz == 0 and y.toConst().isOdd())
        limb_bits // Force an extra limb so that y is even.
    else
        lz;

    x.shiftLeft(x.toConst(), norm_shift);
    y.shiftLeft(y.toConst(), norm_shift);

    const n = x.len - 1;
    const t = y.len - 1;
    const shift = n - t;

    // 1.
    // for 0 <= j <= n - t, set q[j] to 0
    q.len = shift + 1;
    q.positive = true;
    @memset(q.limbs[0..q.len], 0);

    // 2.
    // while x >= y * b^(n - t):
    //    x -= y * b^(n - t)
    //    q[n - t] += 1
    // Note, this algorithm is performed only once if y[t] > base/2 and y is even, which we
    // enforced in step 0. This means we can replace the while with an if.
    // Note, multiplication by b^(n - t) comes down to shifting to the right by n - t limbs.
    // We can also replace x >= y * b^(n - t) by x/b^(n - t) >= y, and use shifts for that.
    {
        // x >= y * b^(n - t) can be replaced by x/b^(n - t) >= y.

        // 'divide' x by b^(n - t)
        var tmp = Mutable{
            .limbs = x.limbs[shift..],
            .len = x.len - shift,
            .positive = true,
        };

        if (tmp.toConst().order(y.toConst()) != .lt) {
            // Perform x -= y * b^(n - t)
            // Note, we can subtract y from x[n - t..] and get the result without shifting.
            // We can also re-use tmp which already contains the relevant part of x. Note that
            // this also edits x.
            // Due to the check above, this cannot underflow.
            tmp.sub(tmp.toConst(), y.toConst());

            // tmp.sub normalized tmp, but we need to normalize x now.
            x.limbs.len = tmp.limbs.len + shift;

            q.limbs[shift] += 1;
        }
    }

    // 3.
    // for i from n down to t + 1, do
    var i = n;
    while (i >= t + 1) : (i -= 1) {
        const k = i - t - 1;
        // 3.1.
        // if x_i == y_t:
        //   q[i - t - 1] = b - 1
        // else:
        //   q[i - t - 1] = (x[i] * b + x[i - 1]) / y[t]
        if (x.limbs[i] == y.limbs[t]) {
            q.limbs[k] = maxInt(Limb);
        } else {
            const q0 = (@as(DoubleLimb, x.limbs[i]) << limb_bits) | @as(DoubleLimb, x.limbs[i - 1]);
            const n0 = @as(DoubleLimb, y.limbs[t]);
            q.limbs[k] = @as(Limb, @intCast(q0 / n0));
        }

        // 3.2
        // while q[i - t - 1] * (y[t] * b + y[t - 1] > x[i] * b * b + x[i - 1] + x[i - 2]:
        //   q[i - t - 1] -= 1
        // Note, if y[t] > b / 2 this part is repeated no more than twice.

        // Extract from y.
        const y0 = if (t > 0) y.limbs[t - 1] else 0;
        const y1 = y.limbs[t];

        // Extract from x.
        // Note, big endian.
        const tmp0 = [_]Limb{
            x.limbs[i],
            if (i >= 1) x.limbs[i - 1] else 0,
            if (i >= 2) x.limbs[i - 2] else 0,
        };

        while (true) {
            // Ad-hoc 2x1 multiplication with q[i - t - 1].
            // Note, big endian.
            var tmp1 = [_]Limb{ 0, undefined, undefined };
            tmp1[2] = addMulLimbWithCarry(0, y0, q.limbs[k], &tmp1[0]);
            tmp1[1] = addMulLimbWithCarry(0, y1, q.limbs[k], &tmp1[0]);

            // Big-endian compare
            if (mem.order(Limb, &tmp1, &tmp0) != .gt)
                break;

            q.limbs[k] -= 1;
        }

        // 3.3.
        // x -= q[i - t - 1] * y * b^(i - t - 1)
        // Note, we multiply by a single limb here.
        // The shift doesn't need to be performed if we add the result of the first multiplication
        // to x[i - t - 1].
        const underflow = llmulLimb(.sub, x.limbs[k..x.len], y.limbs[0..y.len], q.limbs[k]);

        // 3.4.
        // if x < 0:
        //   x += y * b^(i - t - 1)
        //   q[i - t - 1] -= 1
        // Note, we check for x < 0 using the underflow flag from the previous operation.
        if (underflow) {
            // While we didn't properly set the signedness of x, this operation should 'flow' it back to positive.
            llaccum(.add, x.limbs[k..x.len], y.limbs[0..y.len]);
            q.limbs[k] -= 1;
        }
    }

    x.normalize(x.len);
    q.normalize(q.len);

    // De-normalize r and y.
    r.shiftRight(x.toConst(), norm_shift);
    y.shiftRight(y.toConst(), norm_shift);
}