Recurrence relation for the Thue–Morse sequence I made a curious observation. Let $a_n$ be the sequence of numbers determined by a recurrence relation
$$\begin{cases}
\vphantom{\large|}a_0=0\\ 
\vphantom{\large|}a_1=1\\
\vphantom{\Large|}n\,a_n=(5-2 n)\,a_{n-1}+3{\tiny\text{ }}(n-1)\,a_{n-2}+1
\end{cases}\tag{$\small\spadesuit$}$$It also can be generated using an equivalent homogeneous recurrence relation
$$\begin{cases}
\vphantom{\large|}a_0=0\\ 
\vphantom{\large|}a_1=a_2=1\\
\vphantom{\Large|}n\,a_n = (4-n)\,a_{n-1}+(n-2)\left(5{\tiny\text{ }}a_{n-2} - 3{\tiny\text{ }}a_{n-3}\right)
\end{cases}\tag{$\small\clubsuit$}$$
0, 1, 1, 2, 1, 4, -2, 13, -23, 68, -164, 439, -1146, 3067, -8231, 22306, -60791, 166684, ...

It appears that this sequence modulo $2$ gives the Thue–Morse sequence, meaning that if we denote $t_n=(-1)^{a_n}\,$ then it satisfies
$$t_0 = 1,\quad t_n = (-1)^n \, t_{\lfloor n/2\rfloor}.\tag{$\small\diamondsuit$}$$
How can we prove this?
 A: Using (rigorously!) the language of formal power series, this is moderately easy to check.
Let $A(z)=\sum_{n=1}^{\infty}a_n z^n$, considered as a formal power series over $\mathbb{R}$ (for a while). Multiplying $$na_n=(5-2n)a_{n-1}+3(n-1)a_{n-2}+1$$ by $z^{n-1}$ and summing (formally!) over $n>0$ (yes, $n=1$ included), we obtain $$A'(z)=\big(3A(z)-2zA'(z)\big)+3\big(zA(z)+z^2A'(z)\big)+(1-z)^{-1}.$$
Considered as an ODE in the regular sense, with $A(0)=0$, it is easily solved: $$A(z)=\frac{1}{2(1-z)}\left(\sqrt{\frac{1+3z}{1-z}}-1\right),$$ which has a power series expansion convergent for $|z|<1/3$. Being unique, it is exactly our $\sum_{n=1}^{\infty}a_n z^n$.
The function $A(z)$ satisfies $(1-z)^3 A^2(z)+(1-z)^2 A(z)=z$. As an equation for (again) a formal power series $A(z)=\sum_{n=1}^{\infty}a_n z^n$ over $\mathbb{R}$, it has a unique solution (this is seen from the recurrence relation obtained after substituting $A(z)=\sum_{n=1}^{\infty}a_n z^n$ into the equation, and a by-product is that $a_n$ are indeed integers). Now the important observation is that the same is true for a formal power series over the field $\mathbb{F}_2=\mathbb{Z}/2\mathbb{Z}$ (for the same reason).
Thus, it just remains to check that the "generating function" of the Thue-Morse sequence $$B(z)=\sum_{n=0}^{\infty}b_n z^n,\qquad\color{blue}{b_0=0,}\ b_{2n}=b_n,\ b_{2n+1}=1-b_n,$$ considered as a formal power series over $\mathbb{F}_2$, satisfies the same equation. But $$B(z)=\sum_{n=0}^{\infty}b_{2n}z^{2n}+\sum_{n=0}^{\infty}b_{2n+1}z^{2n+1}=(1-z)B(z^2)+z(1-z^2)^{-1},$$ that is, $(1-z^2)B(z)-(1-z)(1-z^2)B(z^2)=z$; since $F(z)=-F(z)$ and $F(z^2)=F^2(z)$ for any formal power series $F(z)$ over $\mathbb{F}_2$, this gives $(1-z)^2 B(z)+(1-z)^3 B^2(z)=z$, as expected.
