Given a function $f(x)$ on [-1,1], the Remez algorithm can find the degree (at most) $n$ polynomial $P_n(x)$ that minimizes the maximum error between $P_n(x)$ and $f(x)$ on that interval. It is an iterative algorithm that finds $n+2$ nodes in the interval with extremal errors of equal magnitude and alternating signs.
Is there a variant of Remez for the special case where $f(x)$ is a polynomial of degree $n+1$ and $P_n(-1) = f(-1)$ and $P_n(1) = f(1)$, ie. endpoints are the same ?
Thanks for all answers. S.
Added: The problem doesn't really depend on the coefficients of $f(x)$ except as to the degree and the separation of the leading coefficient from lower order ones. A reference (unfortunately not easily accessed on line) is Bogacki et al, "Degree reduction of Bézier curves by uniform approximation with endpoint interpolation". From their Abstract: " The approximation of a given Bézier curve $P_n$ of degree $n$ by another of degree $m \lt n$ is called degree reduction. The paper presents two algorithms, which, for 1-degree reduction, produce an approximation which, componentwise, is the best uniform approximation to $P_n$ from the set of all polynomials of degree less than $n$ which interpolate $P_n$ at the endpoints. Algorithms are also presented for multiple degree reductions. In each case, implementation of the algorithm entails multiplying the matrix of control points by a reduction matrix that depends only on the degrees $n$ and $m$, and not on the coefficients of $P_n$. Error bounds and the results of computational experiments are presented."
I don't believe there is any essential difference in applying the Remez algorithm to the situation you describe, except that you use two degrees of freedom in specifying polynomial approximation $P(x)$ to interpolate the endpoints.
That is, if $P(x)$ is to be a polynomial approximation to $f(x)$ of degree (at most) $n$, then to minimize the absolute maximum (aka $L_{\infty}$ norm) error, one selects $n$ interior nodes and solves for $P(x)$ and $\epsilon$ to interpolate both endpoints values and values $f(x_i) \pm \epsilon$ (taking alternating signs) at these interior nodes.
Including $\epsilon$ (which may be positive or negative) and the coefficients of $P(x)$ we have $n+2$ degrees of freedom, which suffices to meet those constraints exactly.
If doing a single point exchange the Remez method continues by finding a new interior node where the absolute error $|P(x) - f(x)|$ is maximized on your interval. Obviously this will not be either endpoint since those values are exact, due to interpolation. The rules for exchanging this new node for an old one are as usual.
An older reference for imposing linear inequality restrictions is available from the AMS online papers, Bruce Chalmer's 1976 "The Remez Exchange Algorithm for Approximation with Linear Restrictions". Example 7 there is a more general case of your problem, imposing interpolation (function equality) at multiple points.
Reducing to cases $f(x) = x^{n+1} - x^{n-1}$
Finding the best minimax ($L^{\infty}$-norm on $[-1,1]$) approximation of a general degree $n+1$ polynomial $f(x)$ by a polynomial $p(x)$ of degree (at most) $n$ and interpolating endpoint values ($f(1)=p(1)$ and ($f(-1)=p(-1)$) can be reduced to the solution of just one case, $f(x) = x^{n+1} - x^{n-1}$.
That particular case has special symmetry (odd or even accordingly as $n+1$ is odd or even), with the consequence that its best minimax approximation of degree at most $n$ will actually be only degree $n-1$. This however does not imply all degree $n+1$ polynomials enjoy this characteristic (i.e. where symmetry is absent), We will illustrate this by a worked example.
First let's explain the reduction. After factoring out a leading coefficient, we may assume that $f(x)$ is a monic polynomial (has leading coefficient 1). The best minimax approximation is suitably preserved, since:
$$ ||\alpha f(x) - \alpha p(x)||_{\infty} = |\alpha| \; ||f(x) - p(x)||_{\infty} $$
Next, subtract a first degree polynomial $t(x)$ that interpolates the values $f(1)$ and $f(-1)$:
$$ t(x) = \frac{f(1)-f(-1)}{2} x + \frac{f(1)+f(-1)}{2} $$
If $t(x)$ is subtracted both from $f(x)$ and any conforming polynomial $p(x)$, then after labelling the results as $f_0(x)$ and $p_0(x)$ resp., the error of approximation is preserved:
$$ ||f_0(x) - p_0(x)||_{\infty} = ||(f(x)-t(x)) - (p(x)-t(x))||_{\infty} = ||f(x) - p(x)||_{\infty} $$
Both $f_0(x) = f(x) - t(x)$ and $p_0(x) = p(x) - t(x)$ now have roots at $x=1,-1$. and so are divisible by $(x^2 - 1)$. Since $n \gt 1$, $f_0$ inherits its leading coefficient from $f$ and so is monic.
We can write $f_0(x) = (x^2 - 1)(x^{n-1} + g(x))$ where $g(x)$ has degree at most $n-2$, and correspondingly write $p_0(x) = (x^2 - 1)q(x)$ where $q(x)$ also has degree at most $n-2$. the minimum max-norm approximation error of $f_0$ agrees with the minimum max-norm approximation error of $x^{n+1}-x^{n-1}=(x^2 - 1)x^{n-1}$:
$$ ||f_0(x) - p_0(x)||_{\infty} = ||(x^2-1)x^{n-1} - (x^2-1)(q(x)-g(x))||_{\infty} $$
because a lower order terms represented by $g(x)$ can be combined with $q(x)$ to produce a conforming approximation $(x^2-1)(q(x)-g(x))$ to $(x^2-1)x^{n-1}$.
Deriving a first few best minimax approximations, $n=2,3,4,5$
If the degrees of polynomials $f(x)$ will be known in advance. the reduced cases $x^{n+1} - x^{n-1}$ lend themselves to precomputation. Their symmetry effectively reduces the computation by half since equioscillations can be mirrored across the origin or the y-axis (depending on the odd or even nature of $n+1$).
For $n=2,3$ the best degree-one reduction (actually, a degree-two reduction) can be found analytically, but according to the last paragraph of the Introduction of this paper, "no explicit formula for the degree reduced polynomial" with endpoint constraints is known for the max-norm (motivating the author of that paper to use square-norms).
The results reported in Table I below for $n=4,5$ were derived numerically. Details should be sufficient for the interested reader to refine their accuracy as needed.
$$ \begin{array}{|c|c|r|c|c|} \hline n & f_0(x) & \text{best } p_0(x) \text{ of degree } \le n & \text{equioscillation nodes} & ||f(x)-p(x)||_\infty \\ \hline 2 & x^3 - x & 0 & \pm \sqrt{3}/3 & 2\sqrt{3}/9 \approx 0.3849 \\ \hline 3 & x^4 - x^2 & (3 - \sqrt{8})(x^2 - 1) & 0, \pm \sqrt{2 - \sqrt{2}} & 3 - \sqrt{8} \approx 0.1716 \\ \hline 4 & x^5 - x^3 & 0.381966 x(x^2-1) & \pm 0.850651,\pm 0.32492 & 0.0803246 \\ \hline 5 & x^6 - x^4 & (0.607695x^2 - 0.0384758)(x^2 - 1) & 0,\pm 0.517638,\pm 0.896575 & 0.0384758 \\ \hline \end{array} $$
Table I.
Readers interested in extending the table may find advantage in locating subsequent equioscillation nodes so that for $n+1$ they are interleaved by those for $n$. Also the equioscillation height $\epsilon$ seems to drop by a little more than half with each increment in $n$. This may help improve the Remez initialization steps.
Worked Example
Consider the (monic) sixth degree polynomial:
$$ f(x) = x^6 + x^5 - 5x^4 - 4x^3 + 6x^2 + 3x - 1 $$
Since $f(1)=1$ and $f(-1)=1$, the endpoints' interpolant is $t(x) = 1$, and:
$$ f_0(x) = f(x) - 1 = (x^2-1)(x^4 + x^3 - 4x^2 - 3x + 3) $$
Thus $f_0(x) = (x^2-1)x^4 + (x^2-1)(x^3 - 4x^2 - 3x + 3)$ has a best minimax approximation among conforming polynomials:
$$ p_0(x) = (mx^2 - \epsilon)(x^2 - 1) + (x^2-1)(x^3 - 4x^2 - 3x + 3) $$
where approximate values $m=0.607695$ and $\epsilon = 0.0384758$ are taken from Table I.
Adding $t(x)$ back in, the best polynomial preserving endpoints for $f(x)$ is:
$$ p(x) = p_0(x) + t(x) = (mx^2 - \epsilon)(x^2 - 1) + (x^2-1)(x^3 - 4x^2 - 3x + 3) + 1 $$
Substituting our approximate solution and simplifying:
$$ p(x) = x^5 - 3.39231x^4 - 4x^3 + 6.35383x^2 + 3x - 1.96152 $$
Note that the degree of $p(x)$ is exactly one less than $f(x)$, not the degree two drop from $x^6 - x^4$ to its best conforming approximation $(mx^2 - \epsilon)(x^2 - 1)$ we have in Table I.