divmod

pub fn divmod(q: ?[]u32, r: ?[]u32, u: []const u32, v: []const u32) !void

File

Code

pub fn divmod(q: ?[]u32, r: ?[]u32, u: []const u32, v: []const u32) !void {
    if (q) |q_| @memset(q_[0..], 0);
    if (r) |r_| @memset(r_[0..], 0);

    if (u.len == 0 or v.len == 0) return error.DivisionByZero;

    var m = u.len - 1;
    var n = v.len - 1;
    while (limb(u, m) == 0) : (m -= 1) {
        if (m == 0) return;
    }
    while (limb(v, n) == 0) : (n -= 1) {
        if (n == 0) return error.DivisionByZero;
    }

    if (n > m) {
        if (r) |r_| @memcpy(r_[0..u.len], u);
        return;
    }

    const s = @clz(limb(v, n));

    var vn: [max_limbs]u32 = undefined;
    var i = n;
    while (i > 0) : (i -= 1) {
        limb_set(&vn, i, shl(u32, limb(v, i), s) | shr(u32, limb(v, i - 1), 32 - s));
    }
    limb_set(&vn, 0, shl(u32, limb(v, 0), s));

    var un: [max_limbs + 1]u32 = undefined;
    limb_set(&un, m + 1, shr(u32, limb(u, m), 32 - s));
    i = m;
    while (i > 0) : (i -= 1) {
        limb_set(&un, i, shl(u32, limb(u, i), s) | shr(u32, limb(u, i - 1), 32 - s));
    }
    limb_set(&un, 0, shl(u32, limb(u, 0), s));

    var j = m - n;
    while (true) : (j -= 1) {
        const uu = (@as(u64, limb(&un, j + n + 1)) << 32) + limb(&un, j + n);
        var qhat = uu / limb(&vn, n);
        var rhat = uu % limb(&vn, n);

        while (true) {
            if (qhat >= (1 << 32) or (n > 0 and qhat * limb(&vn, n - 1) > (rhat << 32) + limb(&un, j + n - 1))) {
                qhat -= 1;
                rhat += limb(&vn, n);
                if (rhat < (1 << 32)) continue;
            }
            break;
        }
        var carry: i64 = 0;
        i = 0;
        while (i <= n) : (i += 1) {
            const p = qhat * limb(&vn, i);
            const t = limb(&un, i + j) - carry - @as(u32, @truncate(p));
            limb_set(&un, i + j, @as(u32, @truncate(@as(u64, @bitCast(t)))));
            carry = @as(i64, @intCast(p >> 32)) - @as(i64, @intCast(t >> 32));
        }
        const t = limb(&un, j + n + 1) -% carry;
        limb_set(&un, j + n + 1, @as(u32, @truncate(@as(u64, @bitCast(t)))));
        if (q) |q_| limb_set(q_, j, @as(u32, @truncate(qhat)));
        if (t < 0) {
            if (q) |q_| limb_set(q_, j, limb(q_, j) - 1);
            var carry2: u64 = 0;
            i = 0;
            while (i <= n) : (i += 1) {
                const t2 = @as(u64, limb(&un, i + j)) + @as(u64, limb(&vn, i)) + carry2;
                limb_set(&un, i + j, @as(u32, @truncate(t2)));
                carry2 = t2 >> 32;
            }
            limb_set(&un, j + n + 1, @as(u32, @truncate(limb(&un, j + n + 1) + carry2)));
        }
        if (j == 0) break;
    }
    if (r) |r_| {
        i = 0;
        while (i <= n) : (i += 1) {
            limb_set(r_, i, shr(u32, limb(&un, i), s) | shl(u32, limb(&un, i + 1), 32 - s));
        }
        limb_set(r_, n, shr(u32, limb(&un, n), s));
    }
}