In this problem we want to explore the use of Newton's method to implement division for doubles. Of course we already have division in hardware, but the method here was often used historically. Also, if one were to implement some high precision floating point routines (in hardware or software), the method given here would give us division quickly and easily once we had the other basics: addition, subtraction, and multiplication. Division by traditional means is much harder than the others.

For division, it suffices to be able to calculate the reciprocal of a positive double. For example, if you want c / d, you can calculate c * (1 / d). So given a d > 0, we want to calculate 1 / d.

In order to use Newton's method, we need an equation involving d whose root is 1 / d. One such equation is f(x) = (1 / x) - d.

1. Work out what the Newton iteration formula would be for f(x) as above. For this formula, you need the standard derivative formula:
   derivative(xⁿ) = n * x^n-1, n != 0, used with n = -1.

   [The resulting formula is very simple and only involves variables x0, d, and x1, and operators multiplication and subtraction. You should work on this part immediately, and we'll make sure you have the correct answer.]

2. Write a simple program (perhaps similar to the first program on square root of 2 or similar to some other program you've seen) that will use Newton's method to calculate 1 / d. Test this program with d = 0.6 and with an initial guess of x0 = 1. [You should get the answer in 6 iterations or so. Of course the answer should be 1 / (6 / 10) = 10 / 6 = 1.666666666666666. It's very important that you program print the successive values of the Newton approximations so you can see what's going on.]

3. In case d satisfies 0.5 ≤ d ≤ 1, we can use the following approximating formula for the initial guess x0:
   x0 = 2.82353 - 1.88235 * d;.

   Modify your program in part (ii) above to use this for the initial value of x0 instead of 1. [Convergence should be one or two steps quicker.] Try another value for d with the same modified program: d = 1/sqrt(2) = 0.7071067811865475. [This also should use the approximating formula.]

4. Now work on the general case in a way that should resemble the discussion of genearal square roots: General square roots. Here an input arbitrary positive d will be normalized so that a new variable d0 satisfies 0.5 ≤ d0 ≤ 1. You should do this by multiplying or dividing by a power of 2 (rather than a power of 4 as with square roots).

   First create a table similar to the first illustration in the link directly above, only for reciprocals rather than square roots. Fill in the various parts of the table, using d in place of r, d0 in place of r0, division by 2ⁿ in place of division by 4ⁿ. You will need to figure out what to do in order to denormalize (not the same as for square root). [You can't make a beautiful ascii table, and it's not needed -- just the information is needed.]

5. Now create a general program to find the reciprocal of an arbitrary positive d, following the pattern used in the programs in the link above, whether Java or C. Try out the following input values of d:
   1. 0.6
   2. 1 / sqrt(2)
   3. 3.141592653589793
   4. 13.0
   5. 0.003247
   6. 2.718281828459045E-8
   7. 4.812118250596034E14

   [You can easily tell what the answers should be.]