>>>>> "Richard" == Richard Fateman <fateman at cs.berkeley.edu> writes:
Richard> I thought we had better root-finding in the FORTRAN
Richard> routines. There is a Laguerre method that is supposedly
I think cpoly is it. But I could be wrong.
Richard> better than Jenkins-Traub. Converting that to lisp
Richard> bfloats might be good. (In fact, I think I have a
Richard> Laguerre method written in Lisp from first principles,
Richard> but it was done by a student who didn't totally inspire
Richard> confidence.)
Here is a 10 minute hack implementation of Laguerre's method, based on
the wikipedia page on Laguerre's method.
Barely tested, but seems to work.
Ray
laguerre1(x0, n, p, p1, p2, x) :=
block([val, g, h, a, denom1, denom2],
val : float(rectform(subst(x=x0, p))),
g : float(rectform(subst(x=x0,p1)/val)),
h : float(rectform(g^2 - subst(x=x0,p2)/val)),
denom1 : g + sqrt((n-1)*(n*h-g^2)),
denom2 : g - sqrt((n-1)*(n*h-g^2)),
if is(abs(denom1)>abs(denom2))
then
float(rectform(n/denom1))
else
float(rectform(n/denom2))
);
laguerre(x0,p,x) :=
block([n : hipow(p,x), r : x0, p1 : diff(p,x), p2 : diff(p, x, 2), dr],
do (dr : laguerre1(r,n,p,p1,p2,x),
r : r - dr,
if abs(dr) < 1e-14*abs(r) then return(r)));