Can the polynomial $1+x+x^2+\dots +x^n$ be factored, for some $n\ge 1$, into a product of two non-constant polynomials with positive coefficients?
Thoughts
It is easy to factor it into polynomials with non-negative coefficients e.g. $$ 1+x+x^2+x^3 = (1+x)(1+0x + x^2), $$ but I have no example with positive coefficients. I believe this should be possible for large $n$, since there are so many (something like $2^{\lceil n/2\rceil-1}$) ways to factor $1+x + x^2 +\dots + x^n$ into a product of two monic polynomials with real coefficients.
Some motivation from probability theory
The question is motivated by Can the sum of two independent r.v.'s with convex support be uniformly distributed?
Namely, we can ask ourselves a discrete counterpart:
Whether a discrete uniform random variable (i.e. the one taking values $0,1,\dots,n$ with equal probabilities) can be decomposed into a sum of independent non-constant random variables, each ranging over a set of consecutive integers?
The link is provided by the probability generating function (pgf): $$ m_Y(x) = \mathbb E[x^{Y}]. $$ If the random variable $Y$ takes values $0,1,\dots,k$ with positive probabilities, then its pgf is a polynomial with positive coefficients: $$ m_Y(x) = \sum_{i=0}^k \mathbb{P}(Y=i) x^i; $$ in particular, for a random variable $U_n$, uniformly distributed on $\{0,1,\dots,n\}$, $$ m_{U_n}(x) = \frac1{n+1}\bigl(1+x+x^2+\dots + x^n\bigr). $$
Since for independent random variables, the pgf of sum if a product of pgfs: $$ m_{Y'+Y''}(x) = \mathbb E[x^{Y'+Y''}] = \mathbb E[x^{Y'}]\mathbb E[x^{Y''}] = m_{Y'}(x)m_{Y''}(x),\tag{1} $$ these two questions are equivalent$^*$.
$^*$Note that in general $(1)$ does not imply the independence of $Y'$ and $Y''$. Nevertheless, if $m_Y$ factors, say, into $m_{Y'}$ and $m_{Y''}$, then $Y$ has the same distribution as the sum of independent copies of $Y'$ and $Y''$, and, indeed, we have the desired decomposition.
$$x^n+\cdots +x+1=\frac{x^{n+1}-1}{x-1}=\frac{\prod_{1\le k\le n+1}(x - e^{2i\pi k/(n+1)})}{x-1}=\prod_{1\le k\lt n+1}(x - e^{2i\pi k/(n+1)})$$ For a factorization $x^n+\dots +x+1=g(x)h(x)$, $g$ and $h$ are products of some $x-\zeta^k$ where $\zeta$ is nth root of unity (up to units (in this case multiplying by positive real numbers) however any such factorization implies one without units). The leading term of both polynomials must be $1$ because it is the product of leading terms that are always $1$. The constant term is also one because it is the product of root of unity so its absolute value is $1$ and it must be real and positive.
$$(x^p+\cdots +1)(x^q+\cdots+1)=x^n+\cdots +x+1$$ Assume without loss of generality that $p\ge q$. Take the constant term of the first polynomial and multiply it with the leading term of the second and you get $x^q$. $$\begin{align} (x^p+\cdots+ ax^q +\cdots +1)(x^q+\cdots+1) &=\cdots +ax^q1+ \cdots +1x^q +\cdots\\ &= \cdots +(\cdots+a+1)x^q \\ \end{align}$$ All other terms must be positive and they can only add to $x^q$, the coefficient should be equal $1$. Therefore all other coefficients of $x^q$ must be zero and $a$ is zero.