% Perform Newton Interpolation on function func
% evaluated at points x(1:n)
% function [c,err,NewtonTable,NewtonTableErr,efactors] = NewtonInterpEbnd(func,x,plt)
%
% Inputs:
%  func = name of function to be interpolated
%  x = array of n points at which to perform interpolation
%      x MUST BE SORTED IN INCREASING ORDER
%  plt = 1 to produce graphics, 0 not to
%
% Outputs:
%  c = n+1 coefficients of interpolation polynomial 
%        p(x) = sum from i=0 to n c(i+1) * (x-x(1))*...*(x-x(i))
%  err = maximum  difference |p(x)-func(x)| at
%        1000 equally spaced points in interpolation interval
%  NewtonTable = divided difference table (see page 355 of text)
%  NewtonTableErr = bounds on accumumlated roundoff in NewtonTable
%  efactors = factors in error bound:
%             efactors(1) = estimate of max f^(n)(x)/n!
%             efactors(2) = estimate roundoff error in efactor(1)
%             efactors(3) = 1/n!
%  if plot = 1, 
%  Figure 1: plot of func(x) (blue), p(x) (red) at 1000 equally space points
%  Figure 2: plot of actual error |f(x)-p(x)| (blue dotted)
%                    Newton error estimate (black dashed) 
%                          = (efactor(1)+efactor(2)) * prod {i=0 to n} (x-x{i})
%                    Newton error bound factor (magenta)
%                          = 1/n! prod {i=0 to n} (x-x(i)) 
%                    roundoff bound from evaluating p(x) (red)
%                    Error bound based on estimated roundoff error 
%                       NewtonTableErr in computing NewtonTable (magenta)
%
function [c,err,NewtonTable,NewtonTableErr,efactors] = NewtonInterpEbnd(func,x,plt)
xmin = min(x);
xmax = max(x);
xx = xmin:(xmax-xmin)/999:xmax;
n = length(x);
x = reshape(x,n,1);
%
nold = n;
n = n + 10;
% Compute Divided Difference Table
% NewtonTable = zeros(n,n+1);
NewtonTable = zeros(n,nold+2);
% NewtonTableErr = zeros(n,n+1);
% NewtonTable(1:n,1) = reshape(x,n,1);
NewtonTable(1:nold,1) = x;
NewtonTable(nold+1:n,1) = rand(10,1)*(xmax-xmin) + xmin;
NewtonTable(1:n,2) = eval([func,'(NewtonTable(1:n,1))']);
NewtonTableErr(1:n,2) = eps*abs(NewtonTable(1:n,2));
% for j=1:n-1,
for j=1:nold,
   NewtonTable(1:n-j,j+2) = ...
       (NewtonTable(2:n-j+1,j+1) - NewtonTable(1:n-j,j+1)) ./ ...
       (NewtonTable(j+1:n,1)   - NewtonTable(1:n-j,1));
   NewtonTableErr(1:n-j,j+2) = ...
       ( eps*abs(NewtonTable(2:n-j+1,j+1) - NewtonTable(1:n-j,j+1)) + ...
       NewtonTableErr(2:n-j+1,j+1) + NewtonTableErr(1:n-j,j+1) ) ./ ...
       (NewtonTable(j+1:n,1)   - NewtonTable(1:n-j,1));
end
maxfn = max(abs(NewtonTable(:,nold+2)));
maxfnerr = max(abs(NewtonTableErr(:,nold+2)));
efactors = [maxfn,maxfnerr,1/gamma(nold+1)];
n = nold;
c = NewtonTable(1,2:n+1);
cbnd2f = NewtonTableErr(1,2:n+1);
cbnd2b = diag(NewtonTableErr(n:-1:1,:),1);
%
% Evaluate Newton Interpolation polynomial at xx, compute error, 
% bound error from roundoff in forming and evaluating the
% Newton polynomial
p = c(n)*ones(size(xx));
pbnd1 = eps*abs(p);
pbnd2f = cbnd2f(n);
pbnd2b = cbnd2b(n);
for i = n:-1:2,
    p = p.*(xx-x(i-1)) + c(i-1);
    pbnd1 = pbnd1.*abs(xx-x(i-1)) + eps*abs(c(i-1));
    pbnd2f = pbnd2f.*abs(xx-x(i-1)) + cbnd2f(i-1);
    pbnd2b = pbnd2b.*abs(xx-x(n+1-(i-1))) + cbnd2b(i-1);
end
fxx = eval([func,'(xx)']);
fx = eval([func,'(x)']);
err = max(abs(fxx-p));
%
% Evaluate error due to roundoff in computing divided differences
%
kk=1;
%  loop through all knots
for i = 1:n
%     look at xx just less than know if j=1
%     look at xx just greater than know if j=n
  for j = 1:2
     if ((i ~= 1 | j ~= 1) & (i ~= n | j ~= n)),
        if j==1,
           mn = (x(i-1)+x(i))/2; mx = x(i);
        else
           mn = x(i); mx = (x(i)+x(i+1))/2;
        end
%       compute range [kold,k] of xx() in range [mn,mx]
        kold = kk+1; if (kold==2), kold =1; end
        while(xx(kk)<mx), kk=kk+1; end
%%      disp('i,j,kold,kk='), [i,j,kold,kk]
%       sort indices of knots in order of increasing distance 
%       selected xx()
        [dist,loc]=sort(abs(x-(mn+mx)/2));
%%      disp('loc='),loc
%       collecte divided differences in this order
        cbndloc(1)=loc(1);
        minl = loc(1);maxl=loc(1);
        for k=2:n,
            if loc(k)<minl,
                cbndloc(k)=cbndloc(k-1)-1; minl=loc(k);
            else
                cbndloc(k)=cbndloc(k-1); maxl=loc(k);
            end
        end
%%      disp('cbndloc='),cbndloc
        for k=1:n, cbnddde(k)=NewtonTableErr(cbndloc(k),k+1); end
%       evaluate error 
        pbnd2(kold:kk) = cbnddde(n);
        for i = n:-1:2,
%           p = p.*(xx-x(i-1)) + c(i-1);
%           pbnd1 = pbnd1.*abs(xx-x(i-1)) + eps*abs(c(i-1));
%           pbnd2f = pbnd2f.*abs(xx-x(i-1)) + cbnd2f(i-1);
%           pbnd2b = pbnd2b.*abs(xx-x(n+1-(i-1))) + cbnd2b(i-1);
            pbnd2(kold:kk) = pbnd2(kold:kk).*abs(xx(kold:kk)-x(loc(i-1))) ...
            + cbnddde(i-1);
        end
%%      disp('pause'),pause
     end
  end
end
%% size(xx),size(pbnd2),size(pbnd2b),size(pbnd2f)
%
% Compute factor of standard error bound for Newton Interpolation:
%   1/n! * prod (i=1 to n) (x - x(i))
%
errfactor = ones(size(xx))/gamma(n+1);
for i=1:n, errfactor = errfactor .* (xx-x(i)); end
errfactor = abs(errfactor);
% Plot
if (plt == 1)
   figure(1)
   hold off, clf
   plot(xx,fxx,'b.',xx,p,'r',x,fx,'bx')
   grid
   title('f(x) = blue, (xi,f(xi)) = blue x, p(x) = red')
   xlabel(['Maximum error = ',num2str(err)])
   figure(2)
   hold off, clf
   semilogy(xx,abs(fxx-p),'b.',xx,pbnd1,'r',xx,errfactor,'m', ...
            xx,(maxfn+maxfnerr)*errfactor*gamma(n+1),'k--',xx,pbnd2,'green'), 
   grid
   title('|f(x) - p(x)| = blue (dotted), Newton error estimate = black dashed, Newton error factor = magenta')
   xlabel('roundoff in p(x) = red, NewtonTable error = green')
   v=axis;
   v(1) = xmin; v(2) = xmax; v(3) = 1e-20*v(4); axis(v);
end