CS3343/3341 Analysis of Algorithms |
Newton's Method to
Perform Division |
Newton's Method: Newton's Method is used to find successive approximations to the roots of a function. If the function is y = f(x) and x0 is close to a root, then we usually expect the formula below to give x1 as a better approximation. Then you plug the x1 back in as x0 and iterate. See Newton's Method for a derivation of the formula below.
The Reciprocal is Enough: Suppose one wants to do division, say by a number d. It suffices to calculate the reciprocal of d, that is, 1/d. Then c/d is calculated as c*(1/d).
Background: Historically, Newton's method has been used in software and in hardware for two impressive tasks: to compute the square root, and (surprisingly) to carry out division. Early in the computer era addition (and subtraction, of course), along with multiplication, were implemented in hardware. Division was quite a bit harder to implement. Think about the multiplication and division by hand that you (surely) learned early in school. Programming multiplication would be relatively straightforward compared with division, which involved successively guessing the next digit to use for the quotient. (Many people now say that they can no longer do division by hand.) To quote from Knuth (Vol. 2, page 251):
Double-precision floating division is the most difficult routine, or at least the most frightening prospect we have encountered so far in this chapter. Actually, it is not terribly complicated once we see how to do it; ...Because of such difficulties, Newton's method was often used.
Newton's Formula for the Reciprocal of d: In order to calculate 1/d, use the function f(x) = 1/x - d, with 1/d as its root. Using Newton's method, one gets the equations:
The Normalization Step: Start with a number d, where we want 1/d. First normalize the input to some small range using powers of 2. With binary floating point numbers, multiplication or division by a power of two only needs to adjust the exponent part of the number, which can be done very quickly. For the reciprocal, normalize into a value d0 in the interval 0.5 <= d0 <= 1 by dividing by a suitable power of two. At the end, denormalize by dividing by the same power of two (since the answer inverted the previous normalization factor). Here is a diagram showing the process, similar to the diagram for square roots:
The Approximation Function:
In the graph at the right, let's use d as the horizontal
variable, and y as the vertical variable.
Then the horizontal numbers go
from 0.5 < d < 1.0,
and this is the range we have normalized d into.
The vertical numbers go
from 1.0 < y < 2.0.
The red line plots the function y = 1/d,
0.5 < d < 1.0.
This is the function we wish to calculate.
In this section we want to get a simple formula to approximate
this function. This formula will give an initial guess to use
for Newton's method. We're going to use a linear equation.
The blue line plots the function:
|
|
A Program for Reciprocal: The program below incorporates normalization and calculates the reciprocal of an arbitrary positive double. The important issue here is that this method does not use division, except for division by a power of two, which can be carried out very efficiently on a binary machine. An actual commercial software division routine would have many improvements and optimizations. The routine below only demonstrates the feasibility. For convenience the input double is given as a command line argument. In the several example outputs, there is one very small positive number and one very large one. In the data below, black digits are correct, while red digits are incorrect. The input values of d are those asked for in Recitation 5. The link Division and Square Root compares the Newton's Method algorithm for division with the one for square root, showing remarkable similarity. The program below used the final linear approximating function (the black line above). I ran this program using the earlier approximating function (the green line above) with the same four inputs. The initial values were quite different, but the final values differed only in trivial ways.
| Newton's Method, division by the double d | Outputs (red = incorrect digits) |
|---|---|
// DivD: use Newton's method for 1/d
// only divisions used are by 2
public class DivD {
private double mul; // denormalizing factor
// normalize to 0.5 <= d <= 1
public double norm(double d) {
// divide d by power of 2
mul = 1.0; // power of 2
double d0 = d;
while (d0 > 1) {
// divide by 2, will
// later divide by 2 again
d0 = d0/2.0;
mul = mul*2.0;
}
while (d0 < 0.5) {
d0 = d0*2.0;
mul = mul/2.0;
}
// d0 now normalized
return d0;
}
// denorm: reverse the normalization
public double denorm(double invD0) {
return invD0/mul;
}
// calculate 1/d0, 0.5 <= d0 <= 1
public double invD(double d0) {
// initial approx: 48/17 - 32/17*d0
double x0 = 2.82353 - 1.88235*d0;
// double x0 = 2.9166666666 - 2.0*d0;
System.out.println("1st approx to 1/d0:"+ x0);
double x1;
for ( ; ; ) {
x1 = x0*(2 - d0*x0);
System.out.println(x1);
if (Math.abs(x1 - x0) < 1.0e-7) break;
x0 = x1;
}
return x1;
}
public static void main(String[] args) {
// calculate 1/d for any d > 0
DivD divD = new DivD();
double d = Double.parseDouble(args[0]);
System.out.println("d:" + d + ", want 1/d");
// normalize to 0.5 <= d <= 1
double d0 = divD.norm(d);
// d0 now normalized
System.out.println("normalized, d0:" + d0 +
", want 1/d0");
double invD0 = divD.invD(d0);
// denoramlize: divide invD0 by mul
double invD = divD.denorm(invD0);
System.out.println("1/" + d + " = \t" + invD);
System.out.println("check on ans:\t" + 1.0/d);
}
}
| d:3.141592653589792, want 1/d (d = π) normalized, d0:0.785398163397448, want 1/d0 1st approx to 1/d0:1.3451357671288138 1.2691797691683022 1.273226599977265 1.2732395446035567 1.2732395447351632 1/3.141592653589792 = 0.3183098861837908 check on ans: 0.3183098861837908 |