% plot error bound and true error in polynomial evaluation
% function [xout, HornerVal, TrueVal, AbsErrorBound, TrueAbsError] = ...
%    polyplot(roots,x)
% Inputs:
%    roots = roots of polynomail (will compute coefficients from these)
%    x = points at which to evaluate polynomial
%        (program will add more points near roots to get xout)
% Outputs: 
%    xout = complete set of points at which polynomial evaluated
%           (input set x plus more clustered near roots)
%    HornerVal = polynomial evaluated using Horner: polyval(poly(roots),x)
%    TrueVal = polynomial evaluated using roots: prod(x-roots(i))
%    AbsErrorBound = computed absolute error bound for Horner
%               = eps*polyval(abs(roots),abs(x))
%    TrueAbsError = true absolute error
%                 = abs(HornerVal - TrueVal)
%    plots of RelErrorBound, TrueRelError vs x
%             AbsErrorBound, TrueAbsError vs x
%             TrueVal
function [xout, HornerVal, TrueVal, AbsErrorBound, TrueAbsError] = ...
            polyplot(roots,x)
coef = poly(roots);
numroots = length(roots);
numx = length(x);
roots_column = reshape(roots,numroots,1);
xout = reshape(x,1,numx);
xout = [xout,roots_column'];
for i=1:numroots, 
   xout=[xout,roots(i)*(1+10.^(-1:-.02:-15)), ...
         roots(i)*(1-10.^(-1:-.02:-15))]; 
end
xout = sort(xout); numx = length(xout);
HornerVal = polyval(coef,xout');
TrueVal = (prod(ones(numroots,1)*xout - roots_column*ones(1,numx)))';
AbsErrorBound = eps*polyval(abs(coef),abs(xout'));
TrueAbsError = abs(HornerVal - TrueVal);
RelErrorBound = AbsErrorBound ./ abs(TrueVal);
TrueRelError = TrueAbsError ./ abs(TrueVal);
figure(1), hold off, clf, subplot(313)
semilogy(xout,TrueRelError,'g.',xout,RelErrorBound,'b')
title('Relative Errorbound in blue, True Relative Error in green')
axis([min(xout),max(xout),1e-20,1e5])
grid
subplot(312)
semilogy(xout,TrueAbsError,'g.',xout,AbsErrorBound,'b')
Arng = [TrueAbsError;AbsErrorBound];
Arng = [1e-5*min(AbsErrorBound(find(AbsErrorBound))),max(Arng)];
Arng = round(log10(Arng));
Arng(1) = 5*floor(Arng(1)/5);
Arng(2) = 5*ceil(Arng(2)/5);
axis([min(xout),max(xout),10^Arng(1),10^Arng(2)]);grid
title(['Absolute Errorbound in blue, True Absolute Error in green', ...
      ' (max=',num2str(10^Arng(2),1),')'])
% ylabel(['error from ',num2str(10^Arng(1)),' to ',num2str(10^Arng(2))])
subplot(311)
vpT = find(TrueVal>0); vmT = find(TrueVal<0);
vpH = find(HornerVal>0); vmH = find(HornerVal<0);
rng = [abs(TrueVal([vpT;vmT]))];
semilogy(xout(vpT),TrueVal(vpT),'g.',xout(vmT),-TrueVal(vmT),'r.')
axis([min(xout),max(xout),min(rng),max(rng)])
% ylabel(['error from ',num2str(min(rng)),' to ',num2str(max(rng))])
grid
title(['True Value > 0 = green, True Value < 0 = red', ...
      ' (max=',num2str(max(rng),1),')'])