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.
