Xref: utzoo sci.physics:17572 sci.math:15879 comp.theory.dynamic-sys:195
Path: utzoo!utgpu!news-server.csri.toronto.edu!rpi!zaphod.mps.ohio-state.edu!swrinde!elroy.jpl.nasa.gov!jato!vsnyder
From: vsnyder@jato.jpl.nasa.gov (Van Snyder)
Newsgroups: sci.physics,sci.math,comp.theory.dynamic-sys
Subject: Re: Newton's method (was square roots by hand or computer)
Message-ID: <1991Mar19.023209.22985@jato.jpl.nasa.gov>
Date: 19 Mar 91 02:32:09 GMT
References: <1991Mar17.120702.9205@jarvis.csri.toronto.edu>
Reply-To: vsnyder@jato.Jpl.Nasa.Gov (Van Snyder)
Organization: Jet Propulsion Laboratory, Pasadena, CA
Lines: 35

In article <1991Mar17.120702.9205@jarvis.csri.toronto.edu> hofbauer@csri.toronto.edu (John Hofbauer) writes:
>    [lots of stuff deleted]
>... The SQRT routine
>on your favourite computer uses Newton's Method. The starting value
>is determined from a simple linear min-max approximation which guarantees
>that the first 3 bits are correct. Three iterations gives full 24-bit
>single precision accuracy.

Argument reduction helps a lot: by fiddling with the exponent, you can
get the fraction between 1/4 and 1 (on a machine with base-2 arithmetic).
Then the result has half the exponent of the adjusted argument, and the
fraction of the result is the square root of the fraction of the argument.
But starting with a number in [1/4,1) allows one to get a very good initial
iterate using approximations in "Computer Approximations" by Hart, et. al.

>Finding roots of arbitrary polynomials gets
>trickier. Some roots lie arbitrarily close to one end of the region of
>local convergence, so even if you started with a value very close to
>the root, the iteration may converge to some other root. This can
>be vary frustrating if you are trying to find all the roots. I speak
>from first hand experience.

Most codes to find all roots of polynomials, e.g. Jenkins and Traub, use a
process called "deflation", equivalent to synthetic division, to get rid of
roots that have been already discovered.  I.e., after finding a1, divide the
given polynomial by (x-a1).  They also use a variation of Newton's method
called Miller's method, that finds complex roots.  To avoid these problems,
at some expense in space, simply form the companion matrix and compute it's
eigenvalues, using a routine from EISPACK.  This is easy since the companion
matrix is already in Hessenberg form.

-- 
vsnyder@jato.Jpl.Nasa.Gov
ames!elroy!jato!vsnyder
vsnyder@jato.uucp