Percy Ludgate
en.wikipedia.org
en.wikipedia.org
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 (?)
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
1/(1 + x) = 1 - x + x^2 - x^3 ...
= (1 - x) (1 + x^2) (1 + x^4) (1 + x^8) ...
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
}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.
- "...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).
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!
"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."