Let there be a polynomial, $f(x)$ and denote $g(x) = f(x^2)$.
Find all monic polynomials, f, with complex coefficients, such that for all $r_i$ satisfying $f(r_i) = 0$, $g(r_i) =0$ as well. (source: An Old (first written in 1985) Indian Math Book Filled With Interesting Problems)
What I've basically broken the problem down to is that if all roots of $f$ are in $g$, then $f$ must be divisible by $g$.
Starting off with degree $1$, it's pretty easy to see that the possibilities are $f(x) = x$ and $f(x) = x-1$.
For degree $2$, I used polynomial long division to get the following polynomials: $f(x) = x^2, x(x-1), x^2 - 1, x^2 + x + 1,$ and $(x-1)^2.$
Afterwards, I realized that if $p(x)$ and $q(x)$ both work, then so much $p(x) \cdot q(x)$.
Furthermore, I saw that all cyclotomic polynomials work.
Other than polynomial long division (which can get quite tedious and long), are there any easier ways to determine all functions $f$ that work for higher degrees (such as $3$ or $4$ or $5$)?
"If all roots of $f$ are in $g$, then $f$ must be divisible by $g$" is not quite right; first, you mean $f$ must divide $g$, and second, it's possible that the multiplicities of the roots of $f$ are bigger than those of $g$. E.g. we could have $f(z) = z^2, g(z) = z$.
Anyway, the question is to find all monic polynomials $f$ such that the set of roots of $f$ is closed under the squaring operation $r \mapsto r^2$. $f$ can have the root $0$ with any multiplicity; we now restrict our attention to nonzero roots $r$ only. If $|r| \neq 1$ then repeated squaring $k \mapsto r^{2^k}$ rapidly sends the absolute value to either $0$ or $\infty$, and in particular would imply that $f$ has infinitely many roots, so all roots of $f$ must have absolute value exactly $1$. If $r = e^{2 \pi i t}$ then $r^{2^k} = e^{2^{k+1} \pi i t}$ which takes on infinitely many distinct values unless $t$ is rational.
So all of the roots of $f$ must be roots of unity, and moreover if a primitive $d^{th}$ root of unity $\zeta_d$ is a root then $\zeta_d^{2^k}$ must all be a root for all $k \ge 1$. So the answer is exactly the products of copies of $f(z) = z$ and polynomials of the form
$$f(z) = \prod_{k=0}^{n-1} (z - \zeta_d^{2^k})^{m_k}$$
where $\zeta_d$ is some primitive $d^{th}$ root of unity and $n$ is the smallest positive integer such that $2^n \equiv 2^i \bmod d$ for some $i < n$ and the multiplicities $m_k \ge 1$ are arbitrary positive integers. ($i = 0$ if $d$ is odd but not otherwise.)
The answer is much nicer if we ask for $f$ to have integer coefficients. Then if $f(z)$ has $\zeta_d$ as a root of multiplicity $m$ it must be divisible by the cyclotomic polynomial $\Phi_d(z)$ with multiplicity $m$. If $d$ is odd then the roots of $\Phi_d(z)$ are already closed under squaring; otherwise, if $d$ is even then squaring the roots of $\Phi_d(z)$ produces exactly the roots of $\Phi_{\frac{d}{2}}(z)$, so $f$ must also be divisible by $\Phi_{\frac{d}{2}}(z)$. Continuing in this way we conclude that $f$ is a product of copies of $z$ and cyclotomic polynomials $\Phi_{d_i}(z)$ with the property that the set of $d_i$ which occur is closed under division by $2$. This can be restated in a slightly nicer way as follows: if $d$ is odd then
$$\Phi_d(z) \Phi_{2d}(z) \Phi_{4d}(z) \dots \Phi_{2^k d}(z) = \Phi_d(z^{2^k})$$
so equivalently $f$ is a product of copies of $z$ and polynomials of the form $\Phi_d(z^{2^k})$ where $d$ is odd.
In particular, your claim that all cyclotomic polynomials work is not quite right, and $\Phi_2(z) = z + 1$ is the smallest counterexample.