The basic idea is to factor the polynomial over $\mathbb F_p$ for one or more suitable primes $p$, and attempt to lift these factorizations to factorizations over $\mathbb Z$. There are bounds on the size of the primes $p$ one needs to consider, making this an effective algorithm.
Over $\mathbb F_p$ one can first do a number of preliminary reductions, to reduce the problem to one where the polynomial $f$ to be factored is square-free and all irreducible factors are of fixed degree $d$. A simple way to proceed is then the (probabilistic) method of Cantor and Zassenhaus: If one takes a random monic polynomial $g$ of degree $\le 2d-1$, then the gcd of $f$ and $g$ will be a non-trivial factor of $f$ with probability about $1/2$. So, just proceed taking random $g$'s and computing the gcd with $f$, until you have a complete factorization. Other algorithms (e.g., Berlekamp) exist for this final step.
The book A Course in Computational Algebraic Number Theory by H. Cohen treats these problems in Chapter 3.4 (Factorization of Polynomials Modulo $p$) and Chapter 3.5 (Factorization of Polynomials over $\mathbb Z$ or $\mathbb Q$). These two chapters don't require any knowledge of algebraic number theory, and are easy to skim through or read.
The polynomials with $f(x)\mid f(x^2)$ are closed under multiplication. In fact, if $f$ is any such polynomial and $g(x)\mid f(x^2)/f(x)$, then $f(x)g(x)$ is also such a polynomial.
WLOG assume $x\nmid f(x)$. The relation $f(x)\mid f(x^2)$ implies
$$ \{\alpha:f(\alpha)=0\}\subseteq \{\beta:f(\beta^2)=0\}=\{\pm\sqrt{\beta}:f(\beta)=0\}.$$
Let $\alpha$ be a zero. Then $\alpha=\pm\sqrt{\beta}$ for some other zero $\beta$, or equivalently $\alpha^2=\beta$. Put another way, the square of any zero is also a zero, so the set of zeros is closed under squaring. Therefore we have a sequence of zeros $\alpha,\alpha^2,\alpha^4,\cdots$ which must eventually terminate since $f$ has finitely many zeros, in which case $\alpha^{2^n}=\alpha^{2^m}$ eventually, so $\alpha^{2^r(2^s-1)}=1$ and thus $\alpha$ is a root of unity.
We can restrict our attention to $f$ that cannot be written as a nontrivial product of other polynomials with this property. I don't think there's a very nice characterization of the possible set of zeros of $f$ beyond "start with a root of unity and keep squaring until you get a repeat." For example, over $\mathbb{C}$ we have that $f(x)=(x-i)(x+1)(x-1)$ is such a polynomial; it includes a kind of "cycle" of length two $\{-1,1\}$ in its zero set, but it also has a kind of "hangnail" at the front, namely $i$. If we think about this in terms of integers mod $n$, we can write $n=2^km$ and use the Chinese Remainder Theorem to track what $x\mapsto 2x$ does to an integer mod $n$; the sequence is eventually periodic but at the beginning the $\mathbb{Z}/2^k\mathbb{Z}$ coordinate may be nonzero.
To get the $f$ with real coefficients, just make sure the set $\{\alpha,\alpha^2,\alpha^4\cdots\}$ is closed under conjugation; if it isn't, then adjoin all their conjugates to construct an $f$ with real coefficients.
And to get $f$ with integer coefficients, if $f$ has an $n$th root of unity as a zero then $f$ is divisible by the cyclotomic polynomial $\Phi_n(x)$. If $n$ is even, then squaring primitive $2n$th roots of unity yield primitive $n$th roots of unity, meaning both $\Phi_{n}(x)$ and $\Phi_{n/2}(x)$ are factors. Writing $n=2^km$, this means it is divisible by $\Phi_{2^km}(x)\Phi_{2^{k-1}m}(x)\cdots\Phi_m(x)$. One can check these polynomials satisfy the condition.
Best Answer
The notation $[x^t]P(x)$ denotes the coefficient of $x^t$ in polynomial $P(x)$.
Suppose that $f(x)=g(x)h(x)$, where $g(x),h(x)$ are monic polynomials, $g(x),h(x)\in\Bbb Z[x]$ and $\deg g>0$, $\deg h>0$, $\deg g+\deg h=n$. Let $g(x)=\sum_k\alpha_kx^k$, $h(x)=\sum_k\beta_kx^k$, where $\alpha_k,\beta_k$ is integer whenever $k\ge0$.
First, we have $\alpha_0\beta_0=3$, so $|\alpha_0|=3$ and $|\beta_0|=1$ or $|\alpha_0|=1$ and $|\beta_0|=3$. WLOG, suppose that $|\alpha_0|=3$ and $|\beta_0|=1$. Let $m$ is the smallest positive integer such that $3\nmid\alpha_m$ (such $m$ exists because $g(x)$ is monic). Now we have $$[z^m]f(x)=\alpha_0\beta_m+\alpha_1\beta_{m-1}+\cdots+\alpha_m\beta_0\equiv\alpha_m\beta_0\not\equiv0\pmod3$$ So $m\ge n-1$, and $\deg g\ge n-1$, thus $\deg h\le 1$, hence $\deg h=1$, and $h(x)=x+\beta_0$, where $\beta_0=\pm1$, so $h(\beta_0)=0$, and $f(\beta_0)=0$, but it's obvious that $f(\pm1)\neq0$, Q.E.D.