I don't think that $p$ being prime makes any difference. The **later thoughts** below suggest a Lagrange interpolation method which is perhaps the same as the resultant method mentioned by Abdelmalek Abdesselam.
Let $f(x)=\sum_{j=0}^nc_jx^j. $ One might require $c_n=1$ or $c_0=1$ but it is perhaps nicer not to. Then setting $u=x^n,$ $\prod_{i=0}^{n-1} f(\omega^ix)=F(u)=\sum_{j=0}^nC_j u^j$.
One can say that
- $C_j$ is a polynomial of degree $n$ in the coefficients $c_0,\cdots,c_n$ where each term has total degree $n$
- $C_j$ has a term $\pm(c_j)^n$ and no term $c_j^{n-1}$.
- $C_{n-j}$ is $C_j$ with $c_k$ replaced by $c_{n-k}$
- the roots of $F$ are the $n$th powers of the roots of $f$.
Here $c_n$ is a constant and the other $c_i$ are symmetric polynomials of the $n$ roots $\alpha_i$ of $f$. The $\alpha_i$ can be thought of as formal variables. Then $c_0,\cdots,c_{n-1}$ are also a basis for the ring of all symmetric polynomials in those variables ($\frac{1}{c_n}$ times a usual basis) . There are other bases for this ring such as $\sigma_k=\sum_{i=1}^n\alpha_i^k$ and the sum of all the terms $\alpha_1^{m_1}\cdots\alpha_n^{m_n}$ with $m_1+\cdots+m_n=k$. Transforming between these bases (more generally, expressing a given symmetric polynomial in terms of them) is a major topic of invariant theory.
In this case, we want to express the $C_i$, which are certain symmetric polynomials in the $\alpha_i^n,$ in terms of the values $c_i$. This must be a well known case. At any rate here are some results:
For $n=4,$
$$C_0=c_0^4$$
$$C_1=(4c_0^3c_4-2c_0^2c_2^2)-(4c_0^2c_3)c_1+(4c_0c_2)c_1^2-c_1^4$$
$$C_2=(6c_0^2c_4^2-8c_0c_1c_3c_4+2c_1^2c_3^2)+(4c_0c_3^2+4c_1^2c_4)c_2-(4c_0c_4+4c_1c_3)c_2^2+c_2^4$$
$$C_3=(4c_4^3c_0-2c_4^2c_2^2)-(4c_4^2c_1)c_3+4(c_4c_2)c_3^2-c_3^4$$
$$C_4=c_4^4$$
while for $n=5$ we have the following (with the others obtainable by symmetry)
$$C_0=c_0^5$$
$$C_1=5c_0^4c_5-(5c_0^3c_4-5c_0^2c_2^2)c_1+(5c_0^2c_3)c_1^2-5c_0^3c_2c_3-(5c_0c_2)c_1^3+c_1^5$$
$${\small C_2=(10c_0^3c_5^2-15c_0^2c_1c_4c_5+5c_0c_1^2c_4^2+5c_0^2c_3^2c_4+10c_0c_1^2c_3c_5-5c_0c_1c_3^3-5c_1^3c_3c_4)}$$ $${ \small+(5c_0^2c_4^2-15c_0^2c_3c_5+5c_1^2c_3^2-5c_0c_1c_3c_4-5c_1^3c_5)c_2}$$ $${\small +(5c_1^2c_4+10c_0c_1c_5+5c_0c_3^2)c_2^2-(5c_0c_4+5c_1c_3)c_2^3+c_2^5}$$
**later thoughts** In general one could consider the problem of producing from a polynomial $f(x)=\sum_{j=0}^nc_jx^j$ a polynomial $F(u)=\sum_{j=0}^nC_j u^j$ whose roots are the $m$th powers $\alpha_i^m$ of the $n$ (unknown) roots of $f$. The solution is the polynomial $\prod_{i=0}^{m-1} f(\tau^ix)$ where now $\tau$ is a primitive $m$th root of unity. This process might spread out the roots. In the case $m=2$ one has (with $ q$ repeated application) $\alpha_i^{2^q}$ and the [Dandelin–Graeffe method][1] for finding the roots of a univariate polynomial. Splitting $f$ into its even and odd parts speeds up the computation. The method was also discovered by [Nikolai Ivanovich Lobachevsky][2] and the linked article suggests that his book *Algebra ili Ichislenie Konechnykh Velichin* discussed the general product $\prod_{i=0}^{m-1} f(\omega^ix)$. Perhaps the appropriate manipulations of symmetric polynomials are discussed there.
Since a polynomial of degree $n$ is determined by its values at $n$ points (plus its leading coeffcient) one has the following method (which I doubt is new):
Let $\zeta$ be a primitive $mn$th root of unity and $\omega$ a primitive $m$th root of unity. Then the desired polynomial $F$ satisfies $F(\zeta^j)=f(\omega^j)$ for $0 \le j \le n-1$. Now Lagrange Interpolation can be used. In this case it might be particularly simple.
[1]: http://en.wikipedia.org/wiki/Graeffe's_method
[2]: http://www.jstor.org/stable/2304571