Nota. La menor potencia $7^n$ que da resto $1$ módulo $1000$ es $n=20$ (un divisor de $\varphi(1000)=400$. Se puede encontrar después de probar un poco si nos damos cuenta de que $7^4\equiv 401\ (\text{mod }1000)$ y que, por consiguiente, potencias de la forma $7^{4k}$ tienen por últimos dígitos $01$. Esto evitaría tener que usar el teorema de Euler-Fermat.
Nota. No es cierto en general que $(x+y)^{2n+1}-x^{2n+1}-y^{2n+1}$ sea múltiplo de $(x+y)^3-x^3-y^3=3xy(x+y)$ y el problema justamente es el $3$ final. No es difícil completar el argumento para ver que esta propiedad es cierta para todo $n$ si, y sólo si, $x\equiv y\equiv 1$ o bien $x\equiv y\equiv 2$ (mod $3$).
Comenzamos observando que \[S(X_0)=\sum_{k=1}^nk^2=\frac{n(n+1)(2n+1)}{6}\lt n(n-1)(n-2)\lt \prod_{k=1}^nk=P(X_0).\] ya que, para $n\gt 5$, se tiene que $n+1\lt 2(n-2)$ y $2n+1\lt 3(n-1)$. Esto nos dice que $D(X_0)=P(X_0)-S(X_0)\gt 0$. Ahora bien, se tiene que \begin{align*} D(X_k)&=P(X_k)-S(X_k)=P(X_{k-1})(P(X_{k-1})-1)-S(X_{k-1})-(P(X_{k-1})-1)^2\\ &=P(X_{k-1})^2-P(X_{k-1})-S(X_{k-1})-P(X_{k-1})^2+2P(X_{k-1})-1=D(X_{k-1})-1. \end{align*} Por lo tanto, la sucesión $D(X_k)$ comienza en un número positivo y decrece en una unidad en cada paso. Deducimos que será cero tras un cierto número de términos y hemos terminado.