Subject: using Newton's method to compute inverses without division
From: rusin
Date: Oct 27 1998 10:37 (timestamp)
To: Math 229 students

RISC (Reduced Instruction Set) Computers often fail to have a Divide
routine implemented in hardware. Obviously b/a = b * (1/a), so division can
be accomplished with multiplication if reciprocals can be computed. Also,
it is simple to move the (binary) decimal point in  a  and  b  simultaneously
so that we may assume  1 < a < 2, say. But how can we effectively compute
reciprocals in this range?

Answer: Newton's method! Using Newton's method to solve
    f(x)=(1/x)-a=0
for x gives x=1/a. The procedure improves an estimate  x*  to a new estimate
	x* - f(x*)/f'(x*) 
	= x* - (1/x* - a)/(-1/x*^2)
	= x* + (1 - a x*) x*
	= (x*)( 2 - a x* )
so that the iterations can be computed using only multiplication.

How fast is the convergence? If  x* = 1/a + eps, then this new terms
works out to  1/a - a eps^2.  For  a  in  (1,2), say, this means each
iteration roughly doubles the number of digits of accuracy!

We can obtain a good starting approximation with some effort, although
for all  a  in (1,2) we can start with  x*=1  and compute 1/a correct to
within 10^(-30) or so in just a dozen iterations of Newton's method.


A similar trick applies for the computation of square roots: applying
Newton's method to the equation  x^2 - a = 0  gives a rapidly-converging
sequence constructed iteratively using  x' = (x* + a/x* ) / 2. But this
requires division. Better to compute a single reciprocal  A = 1/a  and
then solve   f(x) = 1/x^2 - A = 0  (whose solution is  x = sqrt(a))
iteratively using Newton's method: x' = (x*) (3 - A (x*)^2) / 2
(the division by  2  accomplished by a simple bit shift).

dave