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.indexOfScalar(Limb, x.limbs[0..x.len], 0).?;
    // const y_trailing = std.mem.indexOfScalar(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;
    }
}