Division algorithm implementation

[ws-garcia]

I have been reviewing, again, the implementation of large integer division.

Several sources point out that one of the most efficient basic methods is the one proposed by Donald Knuth in his book The Art of Computers Programming (TAOCP), volume 2, Seminumerical Algorithms. In this book, the "algorithm D" is described, available at this link, with some additional comments. A more structured version of the algorithm can be found at Brent, R. P., & Zimmermann, P. (2010). "Modern Computer Arithmetic" (1st ed.). Cambridge University Press and available at this link.

Perhaps it is due to an error of interpretation, but by applying the logic outlined in Knuth's algorithm I have not been able to reach the solution of the following problem:

Suppose that u = 987654 and v = 9873, using a base b = 100 it can be seen that there is no need for "normalization" since the most significant word in the divisor is complies with vn1 > b/2; (98 > 50).

Then we have that the first estimate of the quotient, using the two most significant words of the dividend and the most significant word of the divisor, is qhat = 9876/98 = 100. Looking at the representation chosen for any base b, 0 <= wi <= b - 1 must be satisfied, being wi the ith word of the represented number. In our case, the largest value for the given base is 99, so the estimate of the quotient must be adjusted.

The problem is that after running the algorithm steps qhat is only reduced by 2 and this value is set to qi = 98, still being the wrong result.

Any help here is welcome, @sancarn!

Edit:

Here is the division function from the source you cited. You can contribute porting this to VBA and adding some test. It looks pretty good.

function divMod2(a, b) { // Implementation idea shamelessly stolen from Silent Matt's library http://silentmatt.com/biginteger/
        // Performs faster than divMod1 on larger input sizes.
        var a_l = a.length,
            b_l = b.length,
            result = [],
            part = [],
            base = BASE,
            guess, xlen, highx, highy, check;
        while (a_l) {
            part.unshift(a[--a_l]);
            trim(part);
            if (compareAbs(part, b) < 0) {
                result.push(0);
                continue;
            }
            xlen = part.length;
            highx = part[xlen - 1] * base + part[xlen - 2];
            highy = b[b_l - 1] * base + b[b_l - 2];
            if (xlen > b_l) {
                highx = (highx + 1) * base;
            }
            guess = Math.ceil(highx / highy);
            do {
                check = multiplySmall(b, guess);
                if (compareAbs(check, part) <= 0) break;
                guess--;
            } while (guess);
            result.push(guess);
            part = subtract(part, check);
        }
        result.reverse();
        return [arrayToSmall(result), arrayToSmall(part)];
    }

    function trim(v) {
        var i = v.length;
        while (v[--i] === 0);
        v.length = i + 1;
    }