vsnyder@jato.jpl.nasa.gov (Van Snyder) (03/16/91)
In article <1991Mar15.063023.6605@dsd.es.com> rthomson@dsd.es.com (Rich Thomson) writes: >In article <1991Mar14.113902.97@bsu-ucs.uucp> > 00lhramer@bsu-ucs.uucp (Leslie Ramer) writes: >>I would much rather use Newton's method to find the square root or any other >>root for that matter. [...] > >Well Newton's method has some interesting properties of its own (as >I'm sure many readers of comp.theory.dynamic-sys are aware). > >>Now the only problem is a starting value. A decent starting value is 1. > >If you play with Newton's method long enough you come to realize that, >in general, picking a starting value isn't always obvious. In >general, you pick a starting value that is "close" to the desired root >of the function. If you have no idea where the roots of the function >lie, then you may be in trouble. Newton's method may fail to >converge. In fact, many interesting fractal images have been computed >for simple functions like z^3 = 0, with z being a complex number. > > [stuff deleted] In the pretty inserts following page 114 of James Gleick's book "Chaos" there is a picture of the complex boundaries of the attractors in Newton's method for the zeroes of x^4-1=0. Beautiful graphics. -- vsnyder@jato.Jpl.Nasa.Gov ames!elroy!jato!vsnyder vsnyder@jato.uucp
hofbauer@csri.toronto.edu (John Hofbauer) (03/18/91)
> Now the only problem is a starting value. ...
[lots of stuff deleted]
It must be remembered that Newton's Method is only locally convergent.
If you start within the region of local convergence it will converge
quadratically (the number of correct digits doubles with each iteration).
The problem is finding the region of local convergence. This is generally
impractical and so you just rely on luck. Pick an arbitrary starting
value and let the algorithm bounce around until you fall into the
region of local convergence and then away you go. For square root
you are actually finding the positive root of x**2 - a = 0, which
has a comfortably large region of local convergence. 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. 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.vsnyder@jato.jpl.nasa.gov (Van Snyder) (03/19/91)
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
jroth@allvax.enet.dec.com (Jim Roth) (03/19/91)
In article <1991Mar19.023209.22985@jato.jpl.nasa.gov>, vsnyder@jato.jpl.nasa.gov (Van Snyder) writes... >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. A fast way of priming a sqrt routine is to do a table indexing on the exponent and the first 8 bits of the fraction - this gives about 8 bits of relative accuracy which is then doubled for each further Newton iteration. This is noticibly faster than fiddling with minimax linear approximations in a really optimized routine. An interesting aside, the I860 processor does not have a divide instruction, but it does have a couple of instructions that return an 8 bit accurate floating point approximation to the reciprocal or reciprocal of the root of a number, which are frequently needed in graphics computations. For these, Newton only needs multiplies which are so fast and pipelined that the loss of a real divide is not so bad (the divide on, say, a MIPS R2010 floating point processor, takes numerous cycles, so this isn't as crazy as it sounds.) >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. The trick of using the companion matrix is fairly neat, but a Jenkens Traub rootfinder is faster particularly for high order polynomials. Note also that Jenkens and Traub is really a form of shifted eigenvalue solving anyway, it is in fact equivalent to shifted inverse powering and Rayleigh quotient iteration. Perhaps surprisingly a complex version is slightly more accurate than one specially coded for the roots of real polynomials, though the difference is not large in practice. A Jenkens Traub finder is available from the ACM TOMS on netlib. If you can avoid it, it's best not to have to represent polynomials in polynomial form - for example in filter synthesis one needs roots of the sum of two polynomials with known roots, and it's *much* more accurate to directly iterate on the difference of the product form of those polynomials rather than multiplying them out and forming a result polynomial. - Jim
00lhramer@bsu-ucs.uucp (Leslie Ramer) (03/21/91)
In article <1991Mar19.023209.22985@jato.jpl.nasa.gov>, vsnyder@jato.jpl.nasa.gov (Van Snyder) writes: > 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. I hate to be a "Net-Picker" :-), but it's not Miller's Method, it's Muller's Method (At least according to my numerical analysis text Burden & Faires, 1989, 4th ed. section 2.6, pp 73-84) You might find it interesting to note that Muller's method uses the square root function in it's iterations. Of course, since this method finds complex roots of polynomials, it's the complex square root and negative values are "okay". I'm sure there are better techniques to use to find the square root of a complex number than the (rCIS theta)^1/2 techniques, though. :-) These genearally involve cumbersome trig functions. -- "No one runs so fast as he that is chased." ._ | .-..-. | | .-,.-.-..-..-. 00LHRAMER@bsu-ucs.bsu.edu | +-'`-, |_' .-+| | |+-'| 00LHRAMER@bsu-ucs.UUCP |__`-'`-' | \ `.|| | |`-'| 00LHRAMER@bsuvax1.bitnet