Question: Let $f(X)\in\Bbb{Z}[X]$ be a non-constant polynomial such that $$\left(\frac{f(n)}{p}\right)=\left(\frac{n}{p}\right)$$ where $\left(\frac{\cdot}{p}\right)$ is the Legendre symbol and $p$ is a fixed odd prime. Find all such polynomials $f$.
Basic Observations: We know that $f(n+p)\equiv f(n)\pmod{p}$ for all $n\in\Bbb{Z}$. Then $$\left(\frac{f(n)}{p}\right)=\left(\frac{n}{p}\right)$$ if and only if $$\left(\frac{f(j)}{p}\right)=\left(\frac{j}{p}\right)\tag{1}$$ for all $j\in\{0,1,\ldots,p-1\}$. By Euler's criterion for quadratic residues, we get that $$\boxed{\left(\frac{f(j)}{p}\right)\equiv f(j)^{\frac{p-1}{2}}\pmod{p}}\qquad\boxed{\left(\frac{j}{p}\right)\equiv j^{\frac{p-1}{2}}\pmod{p}}$$ Therefore by condition $(1)$, we get that $$f(j)^{\frac{p-1}{2}}-j^{\frac{p-1}{2}}\equiv0\pmod{p}\;(\,\forall\,0\leq j\leq p-1)$$ Therefore we get that the polynomial $f(X)^{\frac{p-1}{2}}-X^{\frac{p-1}{2}}$ is in the ideal $(X^p-X)$ of $\Bbb{Z}/p\Bbb{Z}[X]$. In fact, this is an iff condition. Any $f$ satisfies the condition of the question if and only if $f(X)^{\frac{p-1}{2}}-X^{\frac{p-1}{2}}$ is in the ideal $(X^p-X)$ of $\Bbb{Z}/p\Bbb{Z}[X]$. I don't know how to proceed further to characterize all such polynomials. It is obvious that any polynomial of the form $aX^{2k+1}$ works, where $a$ is a quadratic residue modulo $p$. But what are all such polynomials? Any help would be highly appreciated. Thanks in advance!
It makes more sense to first work in $\mathbb F_p[x]$. Observe that for $f\in \mathbb F_p[x]$, we have $\left(\frac{f(n)}{p}\right) = \left(\frac{n}{p}\right)$ for all $n\in \mathbb F_p$ if and only if $f(n)=0$ and $\left(\frac{nf(n)}{p}\right)=1$ for all $n\in \mathbb F_p^{\times}$.
We claim that $g\in \mathbb F_p[x]$ satisfies $g(0)=0$ and $\left(\frac{g(n)}{p}\right)=1$ for all $n\in \mathbb F_p^{\times}$ if and only if $g\equiv h^2\mod (x^p-x)$ for some $h\in \mathbb F_p[x]$ with $h(0)=0$ and no zeroes in $\mathbb F_p^\times$. One direction is obvious. For the other, suppose $g$ satisfies $g(0)=0$ and $g(n)$ is a nonzero square for all $n\in \mathbb F_p^{\times}$. Let $\tau : \mathbb F_p^{\times}\to \mathbb F_p^{\times}$ be a function (not necessarily polynomial yet) with $g(n)=\tau(n)^2$. Via Lagrange interpolation, we can construct a polynomial $h$ with $h(0)=0$ and $h(n)=\tau(n)$ for all $n$. Then $g - h^2$ has zeroes everywhere in $\mathbb F_p$, so $g\equiv h^2\mod (x^p-x)$
Going back to $f$, this means that we need $f(0)=0$ and $xf(x)\equiv h(x)^2\mod (x^p-x)$ for some polynomial $h(x)\in \mathbb F_p[x]$ whose only root is $0$. This implies that $x\mid h(x)$, so we can write $h(x) = x r(x)$ for some $r$ which has no roots in $\mathbb F_p^\times$. Then $xf(x)\equiv x^2r(x)^2\mod (x^p-x)$. Since $f(0)=0$, this implies $f(x)\equiv xr(x)^2\mod (x^p-x)$. This classifies all possible $f\in \mathbb F_p[x]$.
So, for polynomials in $\mathbb Z[x]$, the answer is that all possible $f$ are those of the form $$f(x) = xr(x)^2 + p r_1(x) + (x^p-x)r_2(x)$$ for polynomials $r,r_1,r_2$, with the additional restriction that $p\mid r(n)\implies p\mid n$.
The example you gave of $ax^{2k+1}$ is when $r(x) = bx^k$, where $b^2\equiv a\mod p$, and $r_1$ is some suitable constant, and $r_2=0$.