Best polynomial approximation in uniform norm

For a given real-valued function of one real variable on an interval, the code calculates the best approximation in the uniform (max) norm by a polynomial of a given degree. Approximating in a uniform norm is much computationally harder compared to the standard least squares fit, but gives eye pleasing results. It can be viewed as an optimal polynomial interpolation, where the interpolating nodes are not known in advance, but rather determined by the algorithm. The polynomial that best approximates the data (X,Y) in the discrete uniform norm, i.e. the polynomial with the minimum value of max{ | p(x_i) - y_i | , x_i in X }, also known as min-max (or minimax) polynomial, is obtained by the exchange algorithm. Technically, the exchange algorithm requires a finite number of calculations to find the best approximation, but this finite number grows exponentially with the increase of the data points.
The screenshot example:
M = 5; N = 10000; K = 0; EPSH = 10^-12; MAXIT = 10;
X = linspace(-1,1,N); % uniformly spaced nodes on [-1,1]
k=1; Y = abs(X).^k; % the function Y to approximate
[A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT);
p = polyval(A,X); plot(X,Y,X,p) % p is the best approximation

The MATLAB/OCTAVE algorithm and the comments are based on the original FORTRAN code written by Joseph C. Simpson and available on Netlib repository: http://www.netlib.org/toms/501 See also: Communications of the ACM, V14, pp.355-356(1971)

This routine has been mostly developed and modified by students as computer assignments in Approximation Theory courses by Andrew Knyazev, University of Colorado Denver, USA: Team Fall 98 (Revision 1.0): Chanchai Aniwathananon, Crhistopher Mehl, David A. Duran, Saulo P. Oliveira; Team Spring 11 (Revision 1.1): Manuchehr Aminian.