Determine the coefficients of an unknown black-box polynomial Let $p$ be a polynomial of known degree $n$:
$$p(x) = a_0 + a_1 x + \ldots + a_n x^n$$
Suppose we have a magic black box that can evaluate the polynomial for us.  How could one then determine the coefficients $a_i$?  Lagrange polynomials provide a solution, but perhaps we can do better given that we have freedom to choose where to evaluate $p$?
Of course I'm not really sure what I mean by 'better' - minimal number of operations?  Perhaps there are optimal solutions for a given degree, but is tricky in general.  I'll start you off with my solution for the linear case:
$$\begin{align*}a_0 &= p(0) \\ a_1 &= p(1) - a_0\end{align*}$$
And for a quadratic:
$$\begin{align*}
a_0 &= p(0) \\
a_2 &= \frac{1}{2}\left[p(-1) + p(1) - 2a_0\right] \\
a_1 &= p(1) - a_2 - a_0
\end{align*}$$
 A: Just evaluate at the integers from $0$ to $n$. You get the right amount of information, and what you can do with it is calculate the successive differences, $\Delta_0(m)=p(m)$, and for $i\ge0$, $\Delta_{i+1}(m)=\Delta_i(m+1)-\Delta_i(m)$. Of course $\Delta_i$ is defined only for the integers from $0$ to $n-i$. Now $p(x)=\sum_{i=0}^n\Delta_i(0)C_i(x)$, where $C_i(x)$ is the $i$-th binomial polynomial, $C_0=1$, and for $i\ge1$, $C_i(x)=C_{i-1}(x)\cdot(x-i+1)/i$. To get the original coefficients $a_m$, of course you need to expand the binomial polynomials, but if you’re doing many of these problems, this last operation needs to be done only once and the results stored.
A: To expand my comment into an answer:


*

*If you are able to evaluate your black-box polynomial at complex values, then there is an $O(n\log\,n)$ method for generating your coefficients from those values: the fast Fourier transform. As noted by copper.hat in the comments, the matrix underlying the discrete Fourier transform (the Fourier matrix) is a special case of the Vandermonde matrix, the matrix that is involved in the solution of polynomial interpolation problems. More precisely: given your black-box $p(x)$, generate the $n+1$ values $X_k=p\left(\exp\left(-\frac{2\pi ik}{n+1}\right)\right),\quad k=0,\dots,n$ and then perform the inverse transform $c_j=\frac1{n+1}\sum\limits_{k=0}^n X_k \exp\left(\frac{2\pi i jk}{n+1}\right),\quad j=0,\dots,n$. The $c_j$ are your desired polynomial coefficients.

*If you can evaluate only at real values, things are slightly more difficult. One of the more usual things to do is to evaluate your polynomial at the $n+1$ Chebyshev nodes $x_k=\cos\left(\pi\frac{k-\frac12}{n}\right), k=1,\dots,n+1$ and then use the Björck-Pereyra algorithm to solve the underlying Vandermonde linear system. (See the paper I linked to for implementation details.) This is better than the usual choice of equidistant nodes since the underlying Vandermonde linear system has better conditioning (and thus makes for more accurate solutions) when the sample points are clustered appropriately, which is the case for the Chebyshev nodes.
