Polynomial Functional Equation

Let $P \in \mathbb{R}[X]$ be a polynomial satisfying the equation:
$$P(X)P(X+1) = P(X^2+X+1)$$

1. Case where $P$ is constant:
   Let $\lambda \in \mathbb{R}$ such that $P = \lambda$.
   The equation becomes $\lambda^2 = \lambda$, which implies $\lambda \in \{0, 1\}$.
   Thus, the constant solutions are $P = 0$ and $P = 1$.

2. Case where $P$ is non-constant:
   Let $z \in \mathbb{C}$ be a root of $P$. Evaluating the functional equation at $X = z$ and at $X = z - 1$ yields:

$P(z)P(z+1) = 0 \implies P(z^2+z+1) = 0$ (1)
$P(z-1)P(z) = 0 \implies P((z-1)^2+(z-1)+1) = 0 \implies P((z-1)^2+z) = 0$ (2)

Therefore, $z^2+z+1$ and $(z-1)^2+z$ are also roots of $P$.

We will now study the root of maximal modulus:
Let $Z$ denote the set of (complex) roots of $P$, and let $\lambda \in Z$ be a root of maximal modulus.
Applying (1) and (2) to $z = \lambda$, the numbers $\lambda^2+\lambda+1$ and $(\lambda-1)^2+\lambda$ are also roots of $P$. By the maximality of $\vert{}\lambda\vert{}$:
$$\vert{}\lambda^2+\lambda+1\vert{} \le \vert{}\lambda\vert{} \quad \text{and} \quad \vert{}(\lambda-1)^2+\lambda\vert{} \le \vert{}\lambda\vert{}$$By the triangle inequality:
$$2\vert{}\lambda\vert{} = \vert{}(\lambda^2+\lambda+1) - ((\lambda-1)^2+\lambda)\vert{} \le \vert{}\lambda^2+\lambda+1\vert{} + \vert{}(\lambda-1)^2+\lambda\vert{} \le 2\vert{}\lambda\vert{}$$Since the bounds are equal, the case of equality in the triangle inequality implies that there exists a real number $a \ge 0$ such that:
$$\lambda^2+\lambda+1 = -a((\lambda-1)^2+\lambda)$$Furthermore, the equality case forces both $\lambda^2+\lambda+1$ and $(\lambda-1)^2+\lambda$ to have modulus exactly $\vert{}\lambda\vert{}$. This gives:
$$\vert{}\lambda\vert{} = \vert{}-a((\lambda-1)^2+\lambda)\vert{} = a\vert{}\lambda\vert{}$$If $\lambda = 0$, the initial inequality $\vert{}0^2+0+1\vert{} \le \vert{}0\vert{}$ yields $1 \le 0$, which is absurd. Thus $\lambda \neq 0$, which implies $a = 1$. The equation then becomes:
$$\lambda^2+\lambda+1 = -((\lambda-1)^2+\lambda) = -\lambda^2+\lambda-1$$Simplifying gives $2\lambda^2+2 = 0$, which yields $\lambda \in \{i, -i\}$.
Thus, necessarily $\lambda \in \{i, -i\}$, and in particular $\vert{}\lambda\vert{} = 1$.

Now let us show that $P$ has no other complex roots besides $i$ and $-i$.

Let $z \in Z$ be any complex root of $P$. We construct the sequence $(u_n)_{n \in \mathbb{N}}$ defined by:
$$u_0 = z \quad \text{and} \quad u_{n+1} = u_n^2 + u_n + 1$$Since $P(z) = 0 \implies P(z^2+z+1) = 0$, an immediate induction shows that $u_n \in Z$ for all $n \in \mathbb{N}$.
Because $P$ is non-constant, the set of roots $Z$ is finite. By the Pigeonhole Principle, there exist indices $m < n$ such that $u_m = u_n$. Therefore, the sequence $(u_n)_{n \in \mathbb{N}}$ is eventually periodic and enters a cycle of length $p \ge 1$.
Let $N \ge 0$ be the smallest index such that $u_N$ belongs to this cycle, meaning $N = \min \{ n \in \mathbb{N} \mid u_{n+p} = u_n \}$.
For any $k \in \mathbb{N}$, write $u_k = x_k + i y_k$ with $x_k, y_k \in \mathbb{R}$. We have:
$$x_{k+1} = \text{Re}(u_{k+1}) = x_k^2 - y_k^2 + x_k + 1$$Since $\vert{}\lambda\vert{} = 1$ is the maximal modulus among all roots in $Z$, every root $u_k \in Z$ satisfies $\vert{}u_k\vert{} \le 1$. Thus, $x_k^2 + y_k^2 \le 1$, which implies $-y_k^2 \ge x_k^2 - 1$. Substituting this yields:
$$x_{k+1} \ge x_k^2 + (x_k^2 - 1) + x_k + 1 = 2x_k^2 + x_k$$Summing this inequality over one period of the cycle (from $k = N$ to $N + p - 1$):
$$\sum_{k=N}^{N+p-1} x_{k+1} \ge \sum_{k=N}^{N+p-1} (2x_k^2 + x_k)$$Since $u_{N+p} = u_N$, the sums $\sum_{k=N}^{N+p-1} x_{k+1}$ and $\sum_{k=N}^{N+p-1} x_k$ are equal, so subtracting $\sum x_k$ from both sides gives:
$$0 \ge 2 \sum_{k=N}^{N+p-1} x_k^2 \ge 0$$This forces $x_k = 0$ for all elements $u_k$ in the cycle ($k \ge N$). For any element $u_k$ in the cycle ($k \ge N$), since $x_k = 0$ and $x_{k+1} = 0$, the relation $x_{k+1} = x_k^2 - y_k^2 + x_k + 1$ becomes:
$$0 = 0 - y_k^2 + 0 + 1 \implies y_k^2 = 1 \implies y_k \in \{-1, 1\}$$Hence, $u_k \in \{i, -i\}$ for all $k \ge N$.

Finally, we prove by backwards induction that $u_k \in \{i, -i\}$ for all $k \ge 0$:
Suppose $u_{k+1} \in \{i, -i\}$.
Then, $u_{k+1} = u_k^2 + u_k + 1 = i$ or $u_{k+1} = u_k^2 + u_k + 1=-i$
We have then $u_k \in \{-1-i, -1+i, -i, i\}$, since $\vert{}u_k\vert{} \leq 1$, we have $u_k \in \{-i, i\}$

By induction from $N$ down to $0$, we conclude that $u_0 = z \in \{i, -i\}$.

We have thus shown that every complex root of $P$ must be $i$ or $-i$.

Since $P$ has real coefficients, its non-real complex roots must come in conjugate pairs. Since $i$ and $-i$ are complex conjugates of each other, $P$ must be of the form (up to a multiplicative constant):
$$P(X) = c(X^2+1)^n \quad \text{with } n \in \mathbb{N}, c \in \mathbb{R}$$The leading coefficients of $P(X)P(X+1)$ and $P(X^2+X+1)$ must match. Identifying them gives $c^2 = c \implies c \in \{0, 1\}$. Combining this with the constant case gives:
$$S \subset \{0, 1\} \cup \{(X^2+1)^n \mid n \in \mathbb{N}\}$$

Conversely, if $P = 0$ or $P = 1$, the equation is trivially satisfied.
If $P(X) = (X^2+1)^n$ for some $n \in \mathbb{N}$, we verify:
$$P(X)P(X+1) = (X^2+1)^n ((X+1)^2+1)^n = ((X^2+1)(X^2+2X+2))^n$$A direct calculation shows that $(X^2+1)(X^2+2X+2) = (X^2+X+1)^2+1$, so:
$$P(X)P(X+1) = ((X^2+X+1)^2+1)^n = P(X^2+X+1)$$Thus, these polynomials are indeed solutions.

Conclusion:
The set of solutions is
$$S = \{0, 1\} \cup \{(X^2+1)^n \mid n \in \mathbb{N}\}$$