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:

Or just

As with the formula for square roots, this is an amazingly simple formula, given that it produces such good results. Notice also that the formula has only two multiplications and one subtraction. But the overwhelmingly important property of this formula is that it does not use division. This whole approach would be useless if the formula required a division.


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:

    y = −2.0 * d + 3.0.

This is just the chord joining the endpoints of the red function. It is exact at the endpoints and deviates the most in the middle. To get a better approximating function, shift the blue line halfway toward the red curve at the middle, to get the green function:

    y = −2.0 * d + 2.9166666666.

An online source suggested using the following approximation function, plotted at the right as the black line:

    y = −1.88235 * d + 2.82353 .

I don't know the criteria they used to get this formula. Both approximations work fine.

Now start regarding d as a fixed constant, and change to the notation in the first section about Newton's method. We use

    x0 = −1.88235 * d + 2.82353 .

as the initial approximation for Newton's method.

   


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

d:13.0, want 1/d normalized, d0:0.8125, want 1/d0 1st approx to 1/d0:1.294120625 1.2275083439590577 1.230760591145715 1.2307692307085834 1.2307692307692308 1/13.0 = 0.07692307692307693 check on ans: 0.07692307692307693
d:0.003247, want 1/d normalized, d0:0.831232, want 1/d0 1st approx to 1/d0:1.2588604448 1.2004429185386936 1.2030279906583141 1.2030335694228516 1.2030335694487218 1/0.003247 = 307.97659377887277 check on ans: 307.9765937788728
d:2.718281828459045E-8, want 1/d (d = e*10E-8) normalized, d0:0.912104027698647, want 1/d0 1st approx to 1/d0:1.1066309834614518 1.096270065456504 1.0963661621352847 1.0963661705596517 1/2.718281828459045E-8 = 3.6787944117144234E7 check on ans: 3.6787944117144234E7
d:4.812118250596034E14, want 1/d (d=ln(φ)*10E15) normalized, d0:0.8548039166449033, want 1/d0 1st approx to 1/d0:1.2144898475034662 1.1681562359459892 1.1698564572793013 1.1698589355094466 1.1698589355146964 1/4.812118250596034E14 = 2.0780869212350276E-15 check on ans: 2.078086921235028E-15

(Revision date: 2014-06-06. Please use ISO 8601, the International Standard.)