Edit: According to the suggestions in comments, this only holds when f(x) is monic.
Let $p$ be a prime number and $f(x) = \sum_{i=0}^{p} a_{i}x^{i}$ be a monic polynomial with $a_{0}, a_{1}, \ldots, a_{p} \in\mathbb{Z}$. Let $\alpha_{1}, \alpha_{2}, \ldots, \alpha_{p}$ be its roots.
I am trying to show that the polynomial $$g(x) = \frac{x^{p} - \alpha_{1}^{p}}{x - \alpha_{1}} \cdot \frac{x^{p} - \alpha_{2}^{p}}{x - \alpha_{2}} \cdots \frac{x^{p} - \alpha_{p}^{p}}{x - \alpha_{p}}$$ has integer coefficients.
My attempt: I managed to show that the polynomial $$h(x) = (x^{p} - \alpha_{1}^{p})\cdot (x^{p} - \alpha_{2}^{p})\cdots (x^{p} - \alpha_{p}^{p})$$ is a polynomial with integer coefficients. This is because coefficient of $f(x)$ are values of elementary symmetric polynomial in $p$ variables, evaluated at the roots $\alpha_{1}, \alpha_{2}, \ldots, \alpha_{p}$ and coefficient of $h(x)$ are values of same elementary symmetric polynomial evaluated at $\alpha_{1}^{p}, \ldots, \alpha_{p}^{p}$. Then I used Fundamental Theorem of Symmetric Polynomials.
My strategy was to show that when I divide $h(x)$ by $(x-\alpha_{1})\cdots (x-\alpha_{p})$, I get integer polynomial as a result. I thought of doing this partially because it might help me find the coefficents of $g(x)$ too. But now, I am struggling to show that.
Are there any other strategy to this problem (preferably that can help me isolate the coefficients of $g(x)$)?
Under the assumption that $f(x)$ is a monic polynomial, here's a sketch of how one could go about determining the coefficients of this expansion. It should be evident that the expression whose expansion is requested can be written as a polynomial of degree $D=p^2-p$. We write
$$g(x;\vec{\alpha})=P(x;\alpha_1)...P(x;\alpha_p)=\sum_{n=0}^D\Gamma_n(\vec{\alpha})x^n$$
where $P(x;r)=\sum_{n=0}^{p-1}r^{p-1-n}x^n$
It is easy to see by the definition of $g$ that it is invariant under any permutation $S=\{\sigma{(1)},...,\sigma{(p)}\}$ of the set $\{\alpha_1,..., \alpha_p\}$. We readily conclude that the coefficients have the same symmetry:
$$\Gamma_n(S\vec{\alpha})=\Gamma_n(\vec{\alpha})$$
Also the coefficients can be shown to obey the very important scaling relation
$$\Gamma_n(\mu \vec{\alpha})=\mu^{D-n}\Gamma_n(\vec{\alpha})$$
Now, also by the definition of $g$ and simple power counting one should be able to show that the coefficients are polynomials in multiple variables made up only of positive powers of the variables $\alpha_1,...,\alpha_p$. This along with the scaling relation and the symmetry requirement suffices to show that the coefficients must be symmetric polynomials of total degree $D-n$.
With this proven, one can now invoke the fundamental theorem of multivariable symmetric polynomials, which asserts that no matter the form of a symmetric polynomial in many variables, it always has an expression in terms of elementary symmetric polynomials, and since those can be computed to be integer valued by the assumption that the coefficients of the polynomial with roots $\alpha_1,...,\alpha_p$ are integers. Therefore, we conclude that $\Gamma_n \in \mathbb{Z}$.
It takes quite a bit more work to find explicit expressions for the coefficients defined above. Define a basis of symmetric polynomials of total degree $N$ as follows
$$s_{N}^{q_1,..., q_p}(\alpha_1, \alpha_2,... \alpha_p)=\sum_{\text{sym}}\alpha_1^{q_1}... \alpha_p^{q_p}~~,~~ \sum_{i=1}^p q_i=N$$
I conjecture that the coefficients can be expressed as an equal weight linear combination of the polynomials in the basis defined above as follows:
$$\Gamma_n(\vec{\alpha})=\sum_{r\in R} s^{r}_{D-n}(\vec{\alpha})$$
where $r$ is a $p$-index in the set defined by
$$R=\{(r_1,..., r_p):\sum_{i=1}^p{r_i}=D-n,0\leq r_i\leq p-1 \}$$
These sets of p-indices can be generated in Mathematica by evoking the function
Interestingly enough, this symmetric polynomial basis is also good enough to solve the more general problem with
$$P(x;r)=\sum_{n=0}^{p-1}c_{p-1-n} x^n r^{p-1-n}$$
which admits a solution which is a simple generalization of the previous one
$$\Gamma_n(\vec{\alpha})=\sum_{r\in R} c_{r_1}...c_{r_p}s^{r_1,...,r_p}_{D-n}(\vec{\alpha})$$
Note that the coefficients will still be integral if the $c$'s are integral.