function z = polyroot(c,a,b,N)
% function z = polyroot(c,a,b,N=50)
%
% finds a real root in the interval [a,b] of the polynomial with real 
% coefficients c 
%
% the polynomial is p(x) = c(1)x^(m-1) + c(2)x^(m-2) + ... + c(m)x^0
%
% returns the empty matrix if there are no roots
%
% the optional argument N controls the maximum number of iterations
% when searching for a root... searching is only used in certain cases
%
% uses mypolyval.m

% Copyright (C) 2008 Matthew T Harrison
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    See <http://www.gnu.org/licenses/> for a copy of the license.

%%%%%%%%%%%%%%%%%%%%
% INPUT CHECKING   %
%%%%%%%%%%%%%%%%%%%%

m = numel(c);
% remove leading zeros
for k = 1:m
    if c(k), break, end
end
if k > 1
    c = c(k:m);
    m = m-k+1;
end
% check interval
if a > b
    tmp = a;
    a = b;
    b = tmp;
end

%%%%%%%%%%%%%%%%%%%%
% SIMPLE CASES     %
%%%%%%%%%%%%%%%%%%%%

% empty
if m == 0
    z = [];
    return
end
% constant
if m == 1
    if c == 0
        z = a;
    else
        z = [];
    end
    return
end
% linear
if m == 2
    z = -c(2)/c(1);
    if z < a || b < z
        z = [];
    end
    return
end
% quadratic
if m == 3
    A = 2*c(1);
    B = c(2);
    D = B*B-2*A*c(3);
    if D < 0
        z = [];
        return
    end
    D = sqrt(D);
    z = (-B+D)/A;
    if z < a || b < z
        z = (-B-D)/A;
        if z < a || b < z
            z = [];
        end
    end
    return
end
% zero interior terms 
zflag = true;
for k = 2:m-1
    if c(k), zflag = false; break, end
end
if zflag
    D = -c(m)/c(1);
    if D >= 0
        z = D.^(1/(m-1));
    elseif mod(m,2)
        % m-1 is even
        z = [];
        return
    else
        % m-1 is odd
        z = -((-D).^(1/(m-1)));
    end
        if z < a || b < z
        z = [];
    end
    return
end
    
%%%%%%%%%%%%%%%%%%%%
% HIGHER ORDERS    %
%%%%%%%%%%%%%%%%%%%%

% check end points
pa = mypolyval(c,a,m);
if pa == 0
    z = a;
    return
end
pb = mypolyval(c,b,m);
if pb == 0
    z = b;
    return
end
if a == b
    z = [];
    return
end

if pa*pb > 0 
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % FIND ALL OF THE ROOTS USING A COMPANION MATRIX    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    c = c(:).';
    CM = diag(ones(m-2,1),-1);
    CM(1,1:m-1) = c(2:m)/(-c(1));
    r = eig(CM);
    
    % return any one of the valid roots

    z = [];
    for k = 1:m-1
        rk = r(k);
        if ~imag(rk) && rk >= a && rk <= b
            z = rk;
            break
        end
    end
    
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIND A ROOT BY SEARCHING    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

if nargin < 4, N = 50; end

% make p(a) negative
if pa > 0
    tmp = a;
    a = b;
    b = tmp;
    tmp = pa;
    pa = pb;
    pb = tmp;
end

% search (modified false position method)
z = a;
zold = b;
abdiff = 0;
N = -N; % this seems to speed things up a lot

% iterate to convergence
while N < 0
    N = N + 1;
    zold = z;
    
    % find the zero of the line between the two endpoints
    z = a - pa*((b-a)/(pb-pa));
    
    % exit if things aren't changing
    if z == zold, return, end
    
    % adjust repeated refinements of a single endpoint to improve the
    % chances of refining the other endpoint
    if abdiff > 2
        z = .5*(z+b);
        abdiff = 0;
    elseif abdiff < -2
        z = .5*(z+a);
        abdiff = 0;
    end
    
    % evaluate the function
    %pz = mypolyval(c,z,m);  
    % paste in the code for mypolyval since it is evaluated so frequently
    pz = c(1);
    for j = 2:m
        pz = pz.*z + c(j);
    end

    % refine the appropriate endpoint
    if pz > 0
        b = z;
        pb = pz;
        abdiff = abdiff - 1;
    elseif pz < 0
        a = z;
        pa = pz;
        abdiff = abdiff + 1;
    else
        return
    end
end

z = a - pa*((b-a)/(pb-pa));
