This is a different approach than the ones in Proof of relation between divided difference and interpolation polynomial coefficient. Possibly unnecessarily complex, but it gets the job done. We will see that the divided difference is the leading coefficient of the Lagrange polynomial, which will let us prove that the two types of polynomials are equal by induction.
Let $\ell[x_0, x_1, \cdots, x_n]$ be the leading coefficient of the Lagrange polynomial which passes through $f$ at $x_0, x_1, \cdots, x_n$. It is easily seen that
$$
\ell[x_0, x_1, \cdots, x_n] = \sum_{i=0}^n f(x_i) \prod_{\substack{j=0\\j\neq i}}^n \frac1{x_i - x_j}.
$$
We will show that $\ell[x_0, \cdots, x_n] = f[x_0, \cdots, x_n]$. Observe that
\begin{align}
& \ell[x_0, \cdots, x_n] - \ell[x_1, \cdots, x_{n+1}]\\
=\,& \sum_{i=0}^n f(x_i) \prod_{\substack{j=0\\j\neq i}}^n \frac1{x_i - x_j} - \sum_{i=1}^{n+1} f(x_i) \prod_{\substack{j=1\\j\neq i}}^{n+1} \frac1{x_i - x_j} \\
=\,& f(x_0)\prod_{\substack{j=1}}^n \frac1{x_0 - x_j} +
\sum_{i=1}^n f(x_i)
\left(\prod_{\substack{j=0\\j\neq i}}^n \frac1{x_i - x_j} - \prod_{\substack{j=1\\j\neq i}}^{n+1} \frac1{x_i - x_j}\right)
- f(x_{n+1})\prod_{\substack{j=1}}^n \frac1{x_{n+1} - x_j}\\
=\,& f(x_0)\prod_{\substack{j=1}}^n \frac1{x_0 - x_j} +
\sum_{i=1}^n f(x_i)
\left(\frac1{x_i - x_0} - \frac1{x_i - x_{n+1}}\right)\prod_{\substack{j=1\\j\neq i}}^n \frac1{x_i - x_j}
- f(x_{n+1})\prod_{\substack{j=1}}^n \frac1{x_{n+1} - x_j}\\
=\,& f(x_0)\prod_{\substack{j=1}}^n \frac1{x_0 - x_j} +
\sum_{i=1}^n f(x_i)
\frac{x_0 - x_{n+1}}{(x_i - x_0)(x_i - x_{n+1})}\prod_{\substack{j=1\\j\neq i}}^n \frac1{x_i - x_j}
- f(x_{n+1})\prod_{\substack{j=1}}^n \frac1{x_{n+1} - x_j} \\
=\,& f(x_0)\prod_{\substack{j=1}}^n \frac1{x_0 - x_j} +
(x_0 - x_{n+1})\sum_{i=1}^n f(x_i)\prod_{\substack{j=0\\j\neq i}}^{n+1} \frac1{x_i - x_j}
- f(x_{n+1})\prod_{\substack{j=1}}^n \frac1{x_{n+1} - x_j} \\
=\,& (x_0 - x_{n+1}) \left( f(x_0)\prod_{\substack{j=1}}^{n+1} \frac1{x_0 - x_j} +
\sum_{i=1}^n f(x_i)\prod_{\substack{j=0\\j\neq i}}^{n+1} \frac1{x_i - x_j}
+ f(x_{n+1})\prod_{\substack{j=0}}^n \frac1{x_{n+1} - x_j} \right) \\
=\,& (x_0 - x_{n+1}) \sum_{i=0}^{n+1} f(x_i)\prod_{\substack{j=0\\j\neq i}}^{n+1} \frac1{x_i - x_j} \\
=\,&(x_0 - x_{n+1}) \ell[x_0, \cdots, x_{n+1}].
\end{align}
Excuse the mess. The idea isn't too complex, it just takes a ton of space to write out in LaTeX. We just showed that $\ell$ satisfies the divided difference recursive formula
$$
\ell[x_0, \cdots, x_{n+1}] = \frac{\ell[x_0, \cdots, x_n] - \ell[x_1, \cdots, x_{n+1}]}{x_0 - x_{n+1}}.
$$
And since $\ell[x] = f[x] = f(x)$, it must be that $\ell[x_0, \cdots, x_n] = f[x_0, \cdots, x_n]$.
Now we can prove by induction that $q_n$ = $p_n$ for all $n$. For the base case, you can easily verify that $q_0 = p_0$. Now assume that $q_n = p_n$, and consider
$$
r(x) = q_n(x) + c \cdot (x - x_0)(x - x_1)\ldots(x - x_n),
$$
where $c$ is chosen so that $r(x_{n+1}) = f(x_{n+1})$. (Note that there is one and only one such value of $c$.) Clearly, $r$ is the unique $n+1$ degree polynomial which meets $f$ at $x_0, x_1, \cdots, x_{n+1}$. Thus $r = p_{n+1}$. But since $c$ is the leading coefficient of $r = p_{n+1}$, we have $c = \ell[x_0, \cdots x_{n+1}] = f[x_0, \cdots x_{n+1}]$. Thus
$$
p_{n+1}(x) = r(x) = q_n(x) + f[x_0, \cdots x_{n+1}] \cdot (x - x_0)(x - x_1)\ldots(x - x_n) = q_{n+1}(x).
$$
This completes the induction.
Best Answer
Perhaps this is similar to the Lagrange polynomial case.
The Chebyshev polynomials, $T_0, \cdots T_n$ are mutually orthogonal with inner product taken as a sum over the zeros of $T_{n+1}$. The zeros are easily calculated because $T_k(\cos \theta) = \cos k\theta $ (here $x = \cos \theta$). The polynomials are also easily obtained:
$$ T_0 = 1, T_1(x) = x, T_{n+1}(x) = 2x T_{n}(x) - T_{n-1}(x) $$
If interested in differentiation, the derivative of each $T_k$ is simply $k$ time the corresponding Chebyshev polynomial of the second kind (see Wikipedia).
For your problem if you precompute the first $n$ Chebyshev polynomials and save their coefficients, evaluate you $p(x)$ at the $n+1$ zeros and take inner product against the $n+1$ polynomials $ T_0, \cdots, T_{n} $ you quickly obtain coefficients for $p$ in the basis $T_0, \cdots, T_n$. Then use you pre-computed table to give you the coefficients for $p(x)$ in terms of $1, x, \cdots, x^n$. You will need to calculate the norms of each $T_0, \cdots, T_n $ under the inner product - I can't remember if there is a simple form for it or not.
Postscript: I checked. The sum inner product over the $n+1$ zeros of $T_{n+1}$, $\omega_0, \cdots \omega_n$ satisfies $$\langle T_i, T_j \rangle = \sum_{r=0}^n T_i(\omega_r)T_j(\omega_r) = \left\{ \array{0 & i\neq j \\ n+1 & i=j=0 \\ (n+1)/2 & i=j > 0} \right. $$