div

    fn div(q: *Mutable, r: *Mutable, x: *Mutable, y: *Mutable) void {
        assert(!y.eqlZero()); // division by zero
        assert(q != r); // illegal aliasing

        const q_positive = (x.positive == y.positive);
        const r_positive = x.positive;

        if (x.toConst().orderAbs(y.toConst()) == .lt) {
            // q may alias x so handle r first.
            r.copy(x.toConst());
            r.positive = r_positive;

            q.set(0);
            return;
        }

        // Handle trailing zero-words of divisor/dividend. These are not handled in the following
        // algorithms.
        // Note, there must be a non-zero limb for either.
        // const x_trailing = std.mem.findScalar(Limb, x.limbs[0..x.len], 0).?;
        // const y_trailing = std.mem.findScalar(Limb, y.limbs[0..y.len], 0).?;

        const x_trailing = for (x.limbs[0..x.len], 0..) |xi, i| {
            if (xi != 0) break i;
        } else unreachable;

        const y_trailing = for (y.limbs[0..y.len], 0..) |yi, i| {
            if (yi != 0) break i;
        } else unreachable;

        const xy_trailing = @min(x_trailing, y_trailing);

        if (y.len - xy_trailing == 1) {
            const divisor = y.limbs[y.len - 1];

            // Optimization for small divisor. By using a half limb we can avoid requiring DoubleLimb
            // divisions in the hot code path. This may often require compiler_rt software-emulation.
            if (divisor < maxInt(HalfLimb)) {
                lldiv0p5(q.limbs, &r.limbs[0], x.limbs[xy_trailing..x.len], @as(HalfLimb, @intCast(divisor)));
            } else {
                lldiv1(q.limbs, &r.limbs[0], x.limbs[xy_trailing..x.len], divisor);
            }

            q.normalize(x.len - xy_trailing);
            q.positive = q_positive;

            r.len = 1;
            r.positive = r_positive;
        } else {
            // Shrink x, y such that the trailing zero limbs shared between are removed.
            var x0: Mutable = .{
                .limbs = x.limbs[xy_trailing..],
                .len = x.len - xy_trailing,
                .positive = true,
            };

            var y0: Mutable = .{
                .limbs = y.limbs[xy_trailing..],
                .len = y.len - xy_trailing,
                .positive = true,
            };

            divmod(q, r, &x0, &y0);
            q.positive = q_positive;

            r.positive = r_positive;
        }

        if (xy_trailing != 0 and r.limbs[r.len - 1] != 0) {
            // Manually shift here since we know its limb aligned.
            @memmove(r.limbs[xy_trailing..][0..r.len], r.limbs[0..r.len]);
            @memset(r.limbs[0..xy_trailing], 0);
            r.len += xy_trailing;
        }
    }