Solution (source code)

= Solution

Let $S=X+Y$ and $D=X-Y$. Since they are independent, applying the product rule for <moment-generating functions> to $S+D=2X$ gives
$$
M(2t)=M_S(t)M_D(t)=M(t)^2M(t)M(-t)=\boxed{M(t)^3M(-t)}.
$$
All factors are finite and strictly positive for every real $t$ by the stated hypothesis. Applying the same identity to $-t$ and dividing yields $\psi(2t)=\psi(t)^2$, or
$$
\boxed{\psi(t)=\psi(t/2)^2},\qquad \psi(t)=\frac{M(t)}{M(-t)}.
$$
The centered unit-variance <Taylor expansion> from part (ii) gives $M(h)=1+h^2/2+o(h^2)$ and the same expansion for $M(-h)$. Consequently $\psi(h)=1+o(h^2)$.

For the first step of <dyadic rigidity of a moment-generating function>, take logarithms, which are allowed because $\psi>0$. For each fixed $t$ and every integer $m$,
$$
\log\psi(t)=2^m\log\psi(t/2^m).
$$
Since $\log\psi(h)=o(h^2)$, the right side tends to zero: it is $2^m o(t^2/4^m)$. Thus $\psi(t)=1$ for every $t$, and the <moment-generating function> is even.

The original identity now becomes $M(2t)=M(t)^4$. Iterating this time with the correct fourth-power scaling gives
$$
\log M(t)=4^m\log M(t/2^m).
$$
Near zero, $\log M(h)=h^2/2+o(h^2)$, so the right side tends to $t^2/2$. Therefore
$$
\boxed{M(t)=e^{t^2/2}\quad\text{for all real }t}.
$$
This is the <moment-generating function of a standard normal variable>. By the <uniqueness theorem for moment-generating functions>, each of $X$ and $Y$ has the <standard normal distribution>. The proof establishes the <Gaussian characterization by independent sum and difference> directly from the two functional identities and the first two moments.