🚀 go-pugleaf

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