function p = laguerre(a,epsilon,N,p)
%
% p = laguerre(a,epsilon,N,p)
%
% This is Laguerre's Method applied to the polynomial q(x) with coefficients
% (with descending degree) in a and starting the iteration process at p.  If
% no p is % specified, p=0 is assumed.  The process continues until either N
% iterations are performed or |q(p)|<epsilon, and the "updated" p is returned.
%

% Troy Henderson
% Assistant Professor
% University of Mobile
% Updated 7 February, 2008

frac = [0,0.5,0.25,0.75,0.125,0.375,0.625,0.875,1];
mr = length(frac);
mt = round(N/mr);
N = mt*mr;
m = length(a)-1;
if (nargin<4), p = 0; end
converged = 0;
for n=1:N
	[b,d,f] = horner(a,p);
	if (abs(b)<epsilon), break; end
	g =d/b;
	g2 = g^2;
	h = g2 - f/b;
	sq = sqrt((m-1)*(m*h-g2));
	gp = g + sq;
	agp = abs(gp);
	gm = g - sq;
	agm = abs(gm);
	if (gp<gm)
		gp = gm;
		agp = agm;
	end
	if (agp>0)
		dx = m/gp;
	else
		dx = (1+abs(p))*(cos(n)+i*sin(n));
	end
	x1 = p - dx;
	if (mod(n,mt))
		p = x1;
	else
		p = p - frac(n/mt)*dx;
	end
end
b = horner(a,p);
if (abs(b)>=epsilon)
	disp([num2str(b)," is not within accuracy\n"]);
end

endfunction
