Honestly, I find the "original solution" to be quite hand-waving the density argument.
For an alternative direct proof, consider that translation along the $\,x\,$ axis preserves the nature of the roots (real vs. complex) of a real polynomial, so it can be assumed WLOG that $\,a=0\,$, then:
$$p(x)=a_nx^n+a_{n-1}x^{n-1}+\ldots+a_3x^3+a_0 \quad\quad \style{font-family:inherit}{\text{with}} \;\;a_n \ne 0\;\;\style{font-family:inherit}{\text{and}}\;\; p(0)=a_0 \ne0$$
If the roots of $\,p(x)\,$ are $\,x_i \ne 0\,$ then the polynomial having as roots $\,y_i=1 / x_i\,$ is:
$$q(y)=y^np\left(\frac{1}{y}\right)=a_0y^n+a_3y^{n-3}+\ldots+a_{n-1}y+a_n$$
By Vieta's relations $\,\sum_i y_i = 0\,$ and $\,\sum_{i \lt j} y_iy_j = 0\,$, so $\,\sum_i y_i^2 = \left(\sum_i y_i\right)^2-2\sum_{i \lt j} y_iy_j=0\,$. If all $\,y_i\,$ were real, that would imply $\,y_i=0\,$, but $\,a_n \ne 0\,$ so none of the roots can be $\,0\,$. It follows that not all $\,y_i\,$ can be real, and therefore not all of $\,x_i\,$ are real.
We can go after every such polynomial all at once. The key observation is that the condition $f(x)$ divides $f(x^2)$ is equivalent to the following condition (for an algebraically closed field):
If $x$ is a root of $f$, then $x^2$ is also a root of $f$ with the same or greater multiplicity.
This is proven by completely factoring $f$ as $$f(x)=(x-r_1)^{a_1}(x-r_2)^{a_2}\ldots (x-r_k)^{a_k}$$ for distinct roots $r_i$ with multiplicities $a_i$. If we substitute in $x^2$ for $x$ and note $(x^2-k)=(x-\sqrt{k})(x+\sqrt{k})$ we get a complete factoring of $f^2$:
$$f(x^2)=(x-\sqrt{r_1})^{a_1}(x+\sqrt{r_1})^{a_1}(x-\sqrt{r_2})^{a_2}(x+\sqrt{r_2})^{a_2}\ldots (x-\sqrt{r_k})^{a_k}(x+\sqrt{r_k})^{a_k}$$
Note that since the original list of roots was distinct, this list of square roots is also distinct, except if $r_i$ was zero. Since these polynomials are completely factored, $f(x)$ divides $f(x^2)$ if and only if every term $(x-r)^a$ in the factorization of $f(x)$ appears in the factorization of $f(x^2)$ with at least the same multiplicity. Then, noting that the statement we want is trivially true if $x=0$, we are finished.*
We can then spin this into finding every possible such $f$: first, observe that if $x$ is a root, then the sequence $x,x^2,(x^2)^2, ((x^2)^2)^2,\ldots$ must be eventually periodic because all of these must be roots of $f$. This is equivalent to asking that $x$ either be $0$ or be a root of unity.
This can be used to computationally generate every possible polynomial (over $\mathbb C$ - or any field, for that matter) of each degree satisfying the given condition.
There turn out to be a lot of polynomials of this form - though note that every root must have that the sequence of squares $\{x,x^2,x^4,\ldots\}$ has size not exceeding the degree of the polynomial, which ensures that these lists are finite for each degree. For linear terms, you get
$$f(x)=x$$
$$f(x)=x-1$$
Since only $1$ and $0$ can be roots. Then, for quadratic terms, you get, letting $\gamma_{a,n}=e^{2\pi i a/n}$ be a root of unity:
$$f(x)=x^2$$
$$f(x)=x(x-1)=x^2-x$$
$$f(x)=(x-1)^2=x^2-2x+1$$
$$f(x)=(x-1)(x+1)=x^2-1$$
$$f(x)=(x-\gamma_{1,3})(x-\gamma_{2,3})=x^2+x+1$$
For cubic terms, I'll just list a few of the interesting ones, because you can start combining roots from previous "generations" in numerous uninteresting ways -- note, for instance, that we could take any of the quadratic polynomials, take a square root of any of their roots, and add that as a new root, which would give, already, quite a long list! You could also, multiply any of them by $x$ or $x-1$ to get another example. If we want to look at the "primitive" polynomials which are not divisible by any polynomial in a previous generation, you get the following conjugate pair (neither of which are polynomials over $\mathbb R$):
$$f(x)=(x-\gamma_{1,7})(x-\gamma_{2,7})(x-\gamma_{4,7})$$
$$f(x)=(x-\gamma_{3,7})(x-\gamma_{6,7})(x-\gamma_{5,7})$$
For fourth degree, you can expand the list of cubics, similarly. For degree $4$, you get a new real polynomial (which is a cyclotomic polynomial, not coincidentally) and two new complex ones:
$$f(x)=(x-\gamma_{1,5})(x-\gamma_{2,5})(x-\gamma_{4,5})(x-\gamma_{3,5})=1+x+x^2+x^3+x^4$$
$$f(x)=(x-\gamma_{1,15})(x-\gamma_{2,15})(x-\gamma_{4,15})(x-\gamma_{8,15})$$
$$f(x)=(x-\gamma_{14,15})(x-\gamma_{13,15})(x-\gamma_{11,15})(x-\gamma_{7,15})$$
I'm fairly sure that you'll get the complete list of polynomials of degree $n$ recursively as follows:
Take the product of any two polynomials already found whose degrees sum to $n$.
Take any polynomial $f$ found in the previous generation and some $r$ such that $r^2$ is a root of $f$ of greater multiplicity than the multiplicity of $r$ (which may be $0$). Multiply $f$ by $(x-r)$.
Let $r$ be a value satisfying $r^{2^n}=r$ and such that no $n'<n$ satisfies this. Take the polynomial $(x-r)(x-r^2)(x-r^4)\ldots(x-r^{2^{n-1}})$.
though I haven't formally examined this. Note that I've only listed the final case for degrees $3$ and $4$ because the first and second cases are extremely numerous.
A stronger statement
If $r$ is a root of $f(x)$ of multiplicity $a$, and a root of $g(x)-g(r)$ of multiplicity $b$ then $g(r)$ is also a root of $f(x)$ of multiplicity $c$ such that $bc \geq a$.
characterizes solutions to $f | f\circ g$, proved by similar means and gives similar results for how to list such polynomials.
Best Answer
Observe that $\{r-a, r-b, r-c, r-d \}=\{\pm 1, \pm 2\}$. Thus $$r-a + r-b + r-c + r-d=4r-(a+b+c+d)=0$$