With $a_0 := 1$ we have
$$
P(x) = \sum_{k=0}^n (-1)^k a_k x^{n-k} \, .
$$
The idea is to define a polynomial whose roots are all negative or zero, because that translates more easily to a condition on its coefficients. So we define
$$
Q(y) = \sum_{k=0}^n (-1)^k a_k (y-1)^k y^{n-k} \, .
$$
$Q$ is a polynomial of degree at most $n$, $Q(1) = a_0 = 1$, and
$$
Q(y) = (y-1)^n P\left( \frac{y}{y-1}\right)
$$
for $y \ne 1$. If $Q(y) = 0$ then $\frac{y}{y-1}$ is real and in the range $[0, 1]$, and that implies $y \le 0$, so:
- All roots of $Q$ are real and negative or zero.
Also $P(x) > 0$ for $x > 1$ implies $Q(y) > 0$ for $y > 1$, therefore:
- The leading coefficient of $Q$ is positive.
It follows from 1 and 2 that $Q(0) \ge 0$. Repeated application of Rolle's theorem shows that $Q', Q'', \ldots$ also have the properties 1 and 2, so we have
- $Q(y) = \sum_{j=0}^n b_j y^{n-j}$ with nonnegative coefficients $b_0, \ldots b_n$.
Now we want to express the $a_k$ in terms of the $b_j$. Since the mapping $t \mapsto \frac{t}{t-1}$ is self-inverse, it is not difficult to see that
$$
P(x)=(x-1)^n Q\left( \frac{x}{x-1}\right)
$$
for $x \ne 1$, and therefore
$$
\begin{align}
P(x) &= \sum_{k=0}^n (-1)^k x^{n-k} \sum_{j=k}^n \binom j k b_j \\
&= \sum_{j=0}^n b_j x^{n-j} \sum_{k=0}^j \binom j k (-1)^{k} x^{j-k} \\
&= \sum_{k=0}^n (-1)^k x^{n-k} \sum_{j=k}^n \binom j k b_j \, .
\end{align}
$$
Comparing this with the original formula for $P$ gives
- $a_k = \sum_{j=k}^n \binom j k b_j$ for $k=0, \ldots, n$.
Finally we can compute
$$
\sum_{k=i}^n (-1)^{k-i}a_k = \sum_{k=i}^n (-1)^{k-i}\sum_{j=k}^n \binom j k b_j = \sum_{j=i}^n b_j \sum_{k=i}^j \binom j k (-1)^{k-i} \, ,
$$
and since all $b_j$ are nonnegative it suffices to show that
$$
c_{i, j} = \sum_{k=i}^j \binom j k (-1)^{k-i} \ge 0
$$
for $0 \le i \le j \le n$. For $j=0$ is $c_{0, 0} = 1$. For $j > 0$ we
substitute $k=j-l$ and use the formula from Truncated alternating binomial sum:
$$
c_{i, j} = \sum_{l=0}^{j-i} \binom{j}{j-l}(-1)^{j-l-i}
= (-1)^{j-i} \sum_{l=0}^{j-i} \binom j l (-1)^l
= \binom{j-1}{j-i} \ge 0 \, ,
$$
and that finishes the proof.
Since a good proof of this on MathOverflow has already been linked in the comments, I will provide an alternate and original proof.
Firstly, note that we may write $f(x)=x^r g(x)$ for some monic $g(x)\in\mathbb{Z}[x]$ with no roots at $x=0$.
Let $n=\deg(g)$, and let $\alpha_1,\cdots,\alpha_n$ be the roots of $g(x)$. Define the sequence
$$A_k=\alpha_1^k+\alpha_2^k+\cdots+\alpha_n^k$$
It is a fact that the fixed field of the Galois group of a splitting field of a polynomial in $\mathbb{Q}[x]$ is $\mathbb{Q}$, so since $A_k$ is fixed by every automorphism in $\text{Gal}(\mathbb{Q}(\alpha_1,\cdots,\alpha_n)/\mathbb{Q})$, then $A_k\in\mathbb{Q}$. Furthermore, since $\alpha_i$ are algebraic integers, then so is $A_k$, and since $A_k$ is also rational, then it must be an integer. (note that this fact can also be proven by combining of Newton's identities and Vieta's formulas)
Now, note that since $|\alpha_i|\leq 1$, then
$$|A_k|\leq |\alpha_1|^k+|\alpha_2|^k+\cdots+|\alpha_n|^k\leq n$$
Combined with the fact that $A_k$ is an integer, this means that $A_k$ only takes on finitely many values. Writing, $g(x)=x^n-b_{n-1}x^{n-1}-\cdots-b_1x-b_0$, we know that $A_k$ satisfies the linear recurrence
$$A_k=b_{n-1}A_{k-1}+\cdots+b_1A_{k-n+1}+b_0A_{k-n}$$
Since $A_k$ satisfies a linear recurrance and takes on only finitely many values, then $A_k$ is periodic with some period $m$. This means that $A_{mk}=A_0=n$ for all $k\geq 0$. Therefore, for $|x|<1$, we have that
\begin{equation}
\begin{split}
\frac{n}{1-x}&=\sum_{k=0}^\infty nx^k\\
&=\sum_{k=0}^\infty A_{mk}x^k\\
&=\sum_{k=0}^\infty \sum_{i=1}^n \alpha_i^{mk}x^k\\
&=\sum_{i=1}^n\sum_{k=0}^\infty\alpha_i^{mk}x^k\\
&=\sum_{i=1}^n\frac{1}{1-\alpha_i^mx}\\
\end{split}
\end{equation}
and therefore, multiplying both sides by $1-x$, and taking the limit as $x\rightarrow 1^-$, we have that
\begin{equation}
n=\lim_{x\rightarrow 1^-}\sum_{i=1}^n\frac{1-x}{1-\alpha_i^mx}\\
\end{equation}
Notice that if $\alpha_i$ is an $m$-th root of unity, then
$\lim\limits_{x\rightarrow 1^-}\frac{1-x}{1-\alpha_i^mx}=1$
, and otherwise $\lim\limits_{x\rightarrow 1^-}\frac{1-x}{1-\alpha_i^mx}=0$. Combining this with the formula above tells us that $g(x)$ has precisely $n$ roots $\alpha_i$ which are $m$-th roots of unity (this is every root of $g(x)$).
Equivalently, every root of $f(x)=x^r g(x)$ is either $0$ or a root of unity, as desired.
Best Answer
With the useful point of view of @dan_fulea in the comments I was able to prove the statement, I'm sharing my development.
Let $f$ be a polynomial with real coefficients whose roots $x_1, \ldots , x_m$ are all real. If we denote the multiplicity of the root $x_i$ by $r_i$, then
$$f(x)=a(x-x_1)^{r_1}\cdots (x-x_m)^{r_m}.$$
Let's suppose further that just the first $l$ roots of $f$ have multiplicity greater than 1. Then, $$\delta(f)=r_1+ \cdots +r_l +(m-l).$$
Note first that, by Role's theorem, we know about the existence of $m-1$ roots $\xi_i$ of the polynomial $f'(x)$, with $\xi_i \in ]x_i, x_i+1[$ for all $i$. Now, we assert that for each $1 \leq i \leq l$, the number $x_i$ is a root of $f'$.
Indeed, expressing $f$ as the product $f(x)=(x-x_i)^{r_i}g(x)$, by the rule of the derivative of a product we get that
$$f'(x)=r_i(x-x_i)^{r_i-1}g(x)+(x-x_i)^{r_i}g'(x)=(x-x_i)^{r_i-1}P(x),$$ where $P(x)=r_ig(x)+(x-x_i)g'(x)$. Then $x_i$ actually is a root of $f'$; moreover, if $x-x_i$ happens to divide the polynomial $P(x)$, we would get that $x-x_i$ divides $g(x)$, but this contradicts the fact of $r_i$ being the multiplicity of the root $r_i$ of the polynomial $f$.
We found then the roots $\xi_i, \ldots, \xi_{m-1}$ and $x_1, \ldots, x_l$ of the polynomial $f'$. Since $$(m-1)+(r_1-1)+ \cdots +(r_l-1)=(r_1+ \cdots +r_l)+(m-l)-1= \delta(f)-1=\delta(f'),$$
we conclude that $f'(x)=(x-\xi_1)\cdots(x-\xi_{m-1})(x-x_1)^{r_1-1}\cdots(x-x_l)^{r_l-1}$
and we see that all the roots of $f'$ are real numbers.