function [poly,levelerr,xref,iternum] = minimaxd(f,x,n,xrefinit)
% usage: [poly,levelerr] = minimaxd(f,x,n,xrefinit)
% description:  given a vector of abscissas x in increasing order
% from a = x(1) to b = x(m), polynomial degree n < m-2, vector of
% ordinates f(x) of a continuous function f on the interval [a,b],
% and optionally an ititial reference set xrefinit, this routine
% computes poly, the best minimax polynomial approximation from
% P_n to f based on this discrete set of values, final level error
% levelerr, final reference xref and number of iterations iternum.
% Default initial reference is Chebyshev. By choosing x, the user
% implicitly sets a tolerance for the accuracy of the solution.
% Smaller maximum spacing between adjacent nodes of x will yield
% better results. Abscissa x must interlace xref.  For Chebyshev,
% maximum spacing between adjacent nodes of x should be less
% than (b-a)*(1-cos(pi/(n+1)))/2. In all cases f(x) should be
% monotone between nodes of x if possible.

% programmer: t. shores    date: 10/10/09    last rev: 10/13/09

% local variables: 
% iternum, iterlim:
% m:
% a,b:
% xref, xrefndx:
% currndx,ndx:
% j:
% levelerr, newlevelerr:
% xmax, xmaxndx:
% poly, newpoly:
% levelflg:


% set iteration limit..NB: if iternum reaches iterlim, increase it!
iterlim = 4*(n+1); % will do for most, but see above
% minimal bulletproofing
m = length(x);
if (m ~= length(f))
  error('Abscissa/ordinate size mismatch');
elseif (n > m-2)
  error('Insufficient number of nodes')
end
a = x(1);
b = x(m);
% columnize everything
f = f(:);
x = x(:);
% select approximate initial reference
if (nargin == 4)
  if (length(xrefinit) ~= n+2)
    error('Initial reference has incorrect size');
  end
  xref = xrefinit(:);
else
  % no user input, so use Chebyshev reference as initial reference
  xref = (a+(b-a)*(cos((n+1:-1:0)*pi/(n+1))+1)/2)';
end
% find and index nearest x nodes <= initial reference
xrefndx = zeros(n+2,1);
currndx = 1;
for j = 1:(n+2)
  while ((currndx <= m) && (x(currndx)<xref(j)))
    currndx = currndx+1;
  end
  if (currndx > m)
    error('Abscissa vector does not interlace the reference');
  end
  xref(j) = x(currndx);
  xrefndx(j) = currndx;
  % don't reuse this index
  currndx = currndx+1;
end
% main loop of one point exchange algorithm
iternum = 0;
levelflg = 1; % set level flag true
levelerr = -1; % will always level once
while (levelflg)
  iternum = iternum+1;
  % find polynomial that levels current reference
  mtmp = vander(xref);
  polyh = [mtmp(:,2:n+2),(-1).^(1:n+2)']\f(xrefndx);
  newpoly = polyh(1:n+1);
  newlevelerr = abs(polyh(n+2));
  if (newlevelerr <= levelerr)
   % we may be cycling, so don't accept new values and quit
    break;
  else
    % accept new values
    levelerr = newlevelerr;
    poly = newpoly;
  end
  errvec = f - polyval(poly,x);
  abserrvec = abs(errvec);
  maxerr = max(abserrvec);
  % substitute in new reference point
  % locate maximum error
  xmaxndx = find(maxerr == abserrvec,1);
  xmax = x(xmaxndx);
  % locate nearest same signed reference point
  ndx = find(xmax > xref);
  if (length(ndx) == 0)
    % xmax is not larger than any reference point
    if (sign(f(xmaxndx)-polyval(poly,xmax)) == sign(f(xrefndx(1)) ...
       - polyval(poly,xref(1))))
      % replace xref(1) with xmax
      xref(1) = xmax;
      xrefndx(1) = xmaxndx;
    else
      % sign change on left, so shove last xref out
      xref = [xmax;xref(1:(n+1))];
      xrefndx = [xmaxndx;xrefndx(1:(n+1))];
    end
  elseif (length(ndx) == n+2)
    % xmax is larger than any reference point
    if (sign(f(xmaxndx)-polyval(poly,xmax)) == sign(f(xrefndx(n+2)) ...
       - polyval(poly,xref(n+2))))
      % replace xref(n+2) with xmax
      xref(n+2) = xmax;
      xrefndx(n+2) = xmaxndx;
    else
      % sign change on right, so shove first xref out
      xref = [xref(2:(n+2));xmax];
      xrefndx = [xrefndx(2:(n+2));xmaxndx];
    end
  else
    % xmax is in-between reference points
    ndx = ndx(length(ndx));
    % insert new reference point
    if (sign(f(xmaxndx)-polyval(poly,xmax)) == sign(f(xrefndx(ndx)) ...
       - polyval(poly,xref(ndx))))
      % replace xref(ndx) with xmax
      xref(ndx) = xmax;
      xrefndx(ndx) = xmaxndx;
    else
      % replace xref(ndx+1) with xmax
      xref(ndx+1) = xmax;
      xrefndx(ndx+1) = xmaxndx;
    end
  end
  if (iternum == iterlim) % iteration limit reached
    levelflg = 0;
  end
end