Test whether a given function is polynomial. You have a black box function to which you can give real number inputs and from which you can receive real number outputs. How would you test whether it is likely to be a polynomial?

Step 1
I don't know much about "likely". But $N+1$ distinct points uniquely determine a degree $\le <!--\; \le \; -->N$ polynomial by Lagrange interpolation. If the obtained points are $(x_i,y_i)$, for $i=0,\dots <!--\; \dots \; -->N$, then first define:
$$L_{N,j}(x)=\underset{k\ne <!--\; \ne \; -->j}{\overset{N}{\prod <!--\; \prod \; -->}}\frac{x-<!--\; -\; -->x_k}{x_j-<!--\; -\; -->x_k}$$
(which is clearly a polynomial of degree N). The unique polynomial passing through those points is
$$p(x)=\underset{j=0}{\overset{n}{\sum <!--\; \sum \; -->}}{y}_j{L}_{N,j}(x),$$
which is a polynomial of degree at most N. If you compute one more point from your black box, then you can compare it with this polynomial to see whether it fits.
By the same argument, there is no way to deterministically rule out polynomials given a finite set of input values, since the polynomial described above can always be constructed.
Step 2
If you choose a constant step size (like 1), your $p(x_{N+1})$ may be easy to calculate. In that case, you would need to only compute
$$\underset{j=0}{\overset{n}{\sum <!--\; \sum \; -->}}{y}_j\underset{j\ne <!--\; \ne \; -->k}{\prod <!--\; \prod \; -->}(x_{N+1}-<!--\; -\; -->x_j)=\underset{j=0}{\overset{n}{\sum <!--\; \sum \; -->}}{y}_j\underset{j\ne <!--\; \ne \; -->k}{\prod <!--\; \prod \; -->}(N+1-<!--\; -\; -->j)$$
which doesn't sound over-complex.

Step 1
I think the largest problem here is that given things like Taylor series, every function can eventually be closely approximated as a polynomial. As such, even a finite range of measured outputs of a blackbox representing a sine function will might yield a polynomial. Measuring a large range is computationally difficult as higher powers tend to yield very large values for increasing inputs. Anyway, this is what I would do: Suppose you an input set $x\in <!--\; \in \; -->{\mathbb{R}}^{n+1\times <!--\; \times \; -->1}$ and the corresponding output set $S\in <!--\; \in \; -->{\mathbb{R}}^{n+1\times <!--\; \times \; -->1}$, where n is the estimated polynomial order (or better, larger than). Then the polynomial itself can be estimated as follows:
1. create a square matrix $X=\left[\begin{array}{ccccc}{x}^0& x^1& x^2& ...& x^n\end{array}\right]\in <!--\; \in \; -->{\mathbb{R}}^{n\times <!--\; \times \; -->n}$ containing numerical powers of the input vector.
2. Suppose we construct $A={\left[\begin{array}{cccc}{a}_0& a_1& ...& a_n\end{array}\right]}^T$ representing the black box polynomial parameters such that $f(x)=a_0+a_{1}x+a_2{x}^2+...+a_n{x}^n$. Then $f(x)=XA$.
3. use the output vector S to compute A as follows: $A=X^{-<!--\; -\; -->1}S$, provided that X is invertible (which should be the case if all entries of x are unique).
Step 2
If the black box does represent a perfect polynomial, the resulting A vector should be unique and every entry with a higher order than this polynomial should be a flat zero (or approaching floating precision). The uniqueness can be verified by plugging in various different input sets.
This method can become just as troublesome for very large orders of n and can yield just an estimated Taylor series of the function (at least I tested it with a sine function and it did return the Taylor series expansion of it, where the even powers had a parameters smaller than $1e^{-<!--\; -\; -->7}$).