- "...and was the first to exploit a multiplier-accumulator to perform division, using multiplication seeded by reciprocal, via the convergent series 1/(1 + x)".
Am I overlooking it, or is this method missing from Wikipedia's "Division Algorithm" entry?
How it works is easy enough to figure from that description. You rescale your divisor to fit in the range (0.5, 1.0], then expand it as the series 1/(1 - x) = 1 + x + x^2..., which converges, and does so linearly at >1 bit per term (since x < 1/2). What's useful about this, is you can calculate the finite sums with a circuit that's logarithmic depth, or maybe better. You could do the repeated squarings to get: x, x^2, x^4, x^8...; then you have a bunch of parallel hardware multipliers filling out all the intermediate terms, which are products of zero- or one- of each of these basic terms; then you just find the sum with a logarithmic-sized tree of adders.
Somehow, I can't find the name of this algorithm (?)
I haven't thought about it at all, but isn't it somewhat similar to Newton Raphson, in how it approximates a first reciprocal (here with a table), and then converges on better estimates in a few iterations?
Arguably the method for converging looks like it could be different at first glance. You could add it per WP:BOLD and see if anyone objects, the descriptions above are good enough to cite, such that there's no original work in describing the algorithm
late edit: I've tested and confirmed that the naive implementation does in fact work. This C snippet is within at most 1 ULP for most double floats in the range (0.5, 1.0]: I tested a few billion and found no counterexamples.
double recip(double y) {
long double x = ((long double) y) - 1.0;
long double u = 1.0 - x;
long double v = x;
for (int j=1; j<=5; ++j) {
v = v * v; // v = x^(2^j)
u = u * (1.0 + v);
}
return (double)u; // approx. 1.0/y
}
Not only similar, it's identical to Newton Raphson in ideal arithmetic!
The nth order product approximation for the reciprocal of d = 1 - x is
y_n = product_{i=1}{n} (1 + x^{2^(i-1)}) = sum_{i=0}^{2^n - 1} x^i = (1 - x^(2 n)) / (1 - x).
The Newton-Raphson iteration for the reciprocal of d is
z_{n+1} = z_n * (2 - (1-x) * z_n).
Induction with the base case z_0 = 1 = y_0 shows the sequences are equal by inserting y_n for z_n in the Newton Raphson iteration:
z_{n+1} = (1 - x^(2 n)) * (2 - (1 - x^(2 n))) / (1 - x) = ( 1 - x^(4 n) ) / ( 1 - x) = y_{n+1}.
So I guess you could call it an explicit product formula for the Newton Raphson reciprocal.
This 100% checks out! I didn't recognize this at all.
- "...in ideal arithmetic!"
Comparing the two, Newton-Raphson seems rather better conditioned. I guess it's the (1.0 + x^(2^i)) part that holds back the explicit-product version (adding a large float with a small one).
I don't know where you got the idea that there is a first guess that comes from a table. The table lookup just brings x into a suitable range, given A.
Shades of "Pray, Mr. Babbage, if you put into the machine wrong figures, will the right answers come out?"— a Taylor series can only do so much!
This is basically how division on slide rules works.
"If we want to calculate a/b we can also calculate it as a * (1/b); that is the division of two numbers is the equivalent of the first number multiplied by its reciprocal."
Multiplication by the reciprocal isn't the interesting part. Computing the reciprocal using a convergent series is. (Which is not a process that is useful or even feasible on a slide rule, owing to the addition steps required, and the fact that slide rules can trivially perform division as you point out.)
Am I overlooking it, or is this method missing from Wikipedia's "Division Algorithm" entry?
https://en.wikipedia.org/wiki/Division_algorithm?useskin=vec...
How it works is easy enough to figure from that description. You rescale your divisor to fit in the range (0.5, 1.0], then expand it as the series 1/(1 - x) = 1 + x + x^2..., which converges, and does so linearly at >1 bit per term (since x < 1/2). What's useful about this, is you can calculate the finite sums with a circuit that's logarithmic depth, or maybe better. You could do the repeated squarings to get: x, x^2, x^4, x^8...; then you have a bunch of parallel hardware multipliers filling out all the intermediate terms, which are products of zero- or one- of each of these basic terms; then you just find the sum with a logarithmic-sized tree of adders.
Somehow, I can't find the name of this algorithm (?)