[proofplan]
We encode the complete quotients of the continued fraction of $\sqrt{D}$ by integer pairs $(m_n,d_n)$ satisfying
\begin{align*}
\alpha_n=\frac{\sqrt{D}+m_n}{d_n}.
\end{align*}
The continued fraction algorithm gives forward recurrences for these pairs and, because the conjugate complete quotient always lies between $-1$ and $0$, also gives a uniqueness rule for recovering the preceding pair. The assumed period ending in $2a_0$ identifies the terminal pair as $(m_\ell,d_\ell)=(a_0,1)$, which is the same starting data needed for the reversed cycle. Comparing the forward and backward traversals gives $m_{\ell-i}=m_{i+1}$ and $d_{\ell-i}=d_i$, and the floor formula for $a_n$ then forces $a_{\ell-i}=a_i$.
[/proofplan]
[step:Encode the complete quotients by integer pairs]
For $n \geq 0$, define the complete quotient
\begin{align*}
\alpha_n := [a_n;a_{n+1},a_{n+2},\dots],
\end{align*}
where $(a_k)_{k\geq 1}$ denotes the periodic coefficient sequence obtained from the period $a_1,\dots,a_{\ell-1},2a_0$ by the rule
\begin{align*}
a_{k+\ell}=a_k \quad \text{for all } k \geq 1.
\end{align*}
Thus $\alpha_0=\sqrt{D}$ and, for every $n\geq 0$,
\begin{align*}
a_n=\lfloor \alpha_n \rfloor,
\qquad
\alpha_{n+1}=\frac{1}{\alpha_n-a_n}.
\end{align*}
We construct integer sequences $(m_n)_{n\geq 0}$ and $(d_n)_{n\geq 0}$ inductively. Set $m_0:=0$ and $d_0:=1$, so
\begin{align*}
\alpha_0=\frac{\sqrt{D}+m_0}{d_0}.
\end{align*}
Also $d_0\mid D-m_0^2$, so the divisibility invariant holds at $n=0$. Assume for some $n\geq 0$ that $m_n\in\mathbb{Z}$, $d_n\in\mathbb{N}$, $d_n\mid D-m_n^2$, and
\begin{align*}
\alpha_n=\frac{\sqrt{D}+m_n}{d_n}.
\end{align*}
Define
\begin{align*}
m_{n+1}:=d_n a_n-m_n.
\end{align*}
Rationalising the denominator in the continued-fraction recurrence gives
\begin{align*}
\alpha_{n+1}=\frac{1}{\frac{\sqrt{D}+m_n}{d_n}-a_n}.
\end{align*}
Equivalently,
\begin{align*}
\alpha_{n+1}=\frac{d_n}{\sqrt{D}-m_{n+1}}.
\end{align*}
Multiplying numerator and denominator by $\sqrt{D}+m_{n+1}$ gives
\begin{align*}
\alpha_{n+1}=\frac{d_n(\sqrt{D}+m_{n+1})}{D-m_{n+1}^2}.
\end{align*}
Since $m_{n+1}=d_na_n-m_n$, we have
\begin{align*}
m_{n+1}\equiv -m_n \pmod{d_n}.
\end{align*}
Therefore
\begin{align*}
D-m_{n+1}^2\equiv D-m_n^2\equiv 0 \pmod{d_n},
\end{align*}
where the last congruence is the induction hypothesis $d_n\mid D-m_n^2$. Moreover $\alpha_n-a_n>0$, so the identity
\begin{align*}
\alpha_n-a_n=\frac{\sqrt{D}-m_{n+1}}{d_n}
\end{align*}
implies $m_{n+1}<\sqrt{D}$, and hence $D-m_{n+1}^2>0$. Define
\begin{align*}
d_{n+1}:=\frac{D-m_{n+1}^2}{d_n}.
\end{align*}
Then $d_{n+1}$ is a positive integer and
\begin{align*}
\alpha_{n+1}=\frac{\sqrt{D}+m_{n+1}}{d_{n+1}}.
\end{align*}
Thus, by induction, for every $n\geq 0$ the integers $m_n$ and $d_n$ are defined, $d_n>0$, $d_n\mid D-m_n^2$, and
\begin{align*}
m_{n+1}=d_n a_n-m_n.
\end{align*}
Also,
\begin{align*}
d_{n+1}=\frac{D-m_{n+1}^2}{d_n}.
\end{align*}
Finally,
\begin{align*}
a_n=\left\lfloor \frac{a_0+m_n}{d_n}\right\rfloor.
\end{align*}
Indeed, since $a_n=\lfloor \alpha_n\rfloor$ and
\begin{align*}
\alpha_n
=
\frac{a_0+m_n}{d_n}+\frac{\sqrt{D}-a_0}{d_n},
\end{align*}
it is enough to check that the second summand cannot move $\frac{a_0+m_n}{d_n}$ across the next integer. Write $a_0+m_n=q_nd_n+r_n$ with $q_n\in\mathbb{Z}$ and $r_n\in\{0,1,\dots,d_n-1\}$. Since $a_0<\sqrt{D}<a_0+1$ and $d_n>0$,
\begin{align*}
q_n+\frac{r_n}{d_n}
<
\alpha_n
<
q_n+\frac{r_n+1}{d_n}
\leq
q_n+1.
\end{align*}
Thus $\lfloor \alpha_n\rfloor=q_n=\left\lfloor \frac{a_0+m_n}{d_n}\right\rfloor$.
[/step]
[step:Prove the reduced inequalities needed for reversibility]
For every $n\geq 1$,
\begin{align*}
-1<\frac{m_n-\sqrt{D}}{d_n}<0.
\end{align*}
Indeed, for $n=1$,
\begin{align*}
m_1=a_0,
\qquad
d_1=D-a_0^2,
\end{align*}
so
\begin{align*}
\frac{m_1-\sqrt{D}}{d_1}
=
\frac{a_0-\sqrt{D}}{D-a_0^2}
=
-\frac{1}{\sqrt{D}+a_0},
\end{align*}
which lies in $(-1,0)$.
Assume the inequality holds for some $n\geq 1$. Let
\begin{align*}
\beta_n:=\frac{m_n-\sqrt{D}}{d_n}
\end{align*}
denote the algebraic conjugate of $\alpha_n$. Since $\alpha_n>1$ for $n\geq 1$, we have $a_n\geq 1$. The conjugate recurrence is
\begin{align*}
\beta_{n+1}
=
\frac{1}{\beta_n-a_n}.
\end{align*}
Because $-1<\beta_n<0$ and $a_n\geq 1$, the denominator $\beta_n-a_n$ is less than $-1$, and hence
\begin{align*}
-1<\beta_{n+1}<0.
\end{align*}
Induction proves the claim.
Consequently, for every $n\geq 1$,
\begin{align*}
0<\sqrt{D}-m_n<d_n.
\end{align*}
Since $m_n$ is an integer and $a_0<\sqrt{D}<a_0+1$, this implies
\begin{align*}
m_n\leq a_0
\qquad\text{and}\qquad
0\leq a_0-m_n<d_n.
\end{align*}
[guided]
The purpose of this step is to isolate the exact inequalities that make the continued fraction cycle reversible. We need to know that, at each stage, the integer $m_n$ sits just below $\sqrt{D}$ in a controlled way.
For $n\geq 1$, define the algebraic conjugate of the complete quotient $\alpha_n=(\sqrt{D}+m_n)/d_n$ by
\begin{align*}
\beta_n:=\frac{m_n-\sqrt{D}}{d_n}.
\end{align*}
We prove that $\beta_n$ always lies in the interval $(-1,0)$.
At the first step, the recurrence gives
\begin{align*}
m_1=d_0a_0-m_0=a_0,
\qquad
d_1=\frac{D-m_1^2}{d_0}=D-a_0^2.
\end{align*}
Therefore
\begin{align*}
\beta_1
=
\frac{a_0-\sqrt{D}}{D-a_0^2}
=
\frac{a_0-\sqrt{D}}{(\sqrt{D}-a_0)(\sqrt{D}+a_0)}
=
-\frac{1}{\sqrt{D}+a_0}.
\end{align*}
Since $D$ is not a square, $\sqrt{D}>a_0\geq 1$ unless $D=1$, which is excluded because $1$ is a square. Hence $\sqrt{D}+a_0>1$, so $-1<\beta_1<0$.
Now suppose $-1<\beta_n<0$ for some $n\geq 1$. The complete quotient recurrence is
\begin{align*}
\alpha_{n+1}=\frac{1}{\alpha_n-a_n}.
\end{align*}
Taking algebraic conjugates gives
\begin{align*}
\beta_{n+1}=\frac{1}{\beta_n-a_n}.
\end{align*}
Because $n\geq 1$, the complete quotient $\alpha_n=[a_n;a_{n+1},\dots]$ has positive integer part, so $a_n\geq 1$. Combining this with $-1<\beta_n<0$ gives
\begin{align*}
\beta_n-a_n<-1.
\end{align*}
Taking reciprocals reverses the scale on this negative interval and yields
\begin{align*}
-1<\frac{1}{\beta_n-a_n}<0.
\end{align*}
Thus $-1<\beta_{n+1}<0$.
The inequality $-1<\beta_n<0$ is exactly
\begin{align*}
-1<\frac{m_n-\sqrt{D}}{d_n}<0.
\end{align*}
Multiplying by $d_n>0$ gives
\begin{align*}
0<\sqrt{D}-m_n<d_n.
\end{align*}
Since $m_n$ is an integer and $a_0=\lfloor\sqrt{D}\rfloor$, the inequality $m_n<\sqrt{D}$ forces $m_n\leq a_0$. Subtracting $m_n$ from $a_0$ and using $\sqrt{D}>a_0$ gives
\begin{align*}
0\leq a_0-m_n<\sqrt{D}-m_n<d_n.
\end{align*}
These are the reduced inequalities used later to identify the unique predecessor of a state.
[/guided]
[/step]
[step:Identify the terminal pair of the period]
Because the period is
\begin{align*}
a_1,a_2,\dots,a_{\ell-1},2a_0,
\end{align*}
we have
\begin{align*}
\alpha_\ell=[2a_0;\overline{a_1,\dots,a_{\ell-1},2a_0}].
\end{align*}
Also
\begin{align*}
\alpha_1=[a_1;\dots,a_{\ell-1},2a_0,\overline{a_1,\dots,a_{\ell-1},2a_0}]
=\frac{1}{\sqrt{D}-a_0}.
\end{align*}
Hence
\begin{align*}
\alpha_\ell
=
2a_0+\frac{1}{\alpha_1}
=
2a_0+\sqrt{D}-a_0
=
\sqrt{D}+a_0.
\end{align*}
Since $\alpha_\ell=(\sqrt{D}+m_\ell)/d_\ell$ with integers $m_\ell$ and $d_\ell>0$, comparison of the rational and irrational parts gives
\begin{align*}
m_\ell=a_0,
\qquad
d_\ell=1.
\end{align*}
[/step]
[step:Show that reversed terminal states match shifted forward states]
Define the state domain
\begin{align*}
\mathcal{A}:=\{(M,E)\in\mathbb{Z}\times\mathbb{N}: E\mid D-M^2\}.
\end{align*}
Define the predecessor modulus map $\delta:\mathcal{A}\to\mathbb{N}$ by
\begin{align*}
\delta(M,E):=\frac{D-M^2}{E}.
\end{align*}
Define the predecessor map $P:\mathcal{A}\to\mathbb{Z}\times\mathbb{N}$ by
\begin{align*}
P(M,E):=(r,\delta(M,E)),
\end{align*}
where $r$ is the unique integer satisfying
\begin{align*}
a_0-\delta(M,E)<r\leq a_0,
\qquad
r\equiv -M \pmod{\delta(M,E)}.
\end{align*}
There is exactly one such integer because the interval $(a_0-\delta(M,E),a_0]$ has length $\delta(M,E)$ and contains exactly one representative of each residue class modulo $\delta(M,E)$.
For the actual continued-fraction states, the recurrence gives
\begin{align*}
d_{n+1}=\frac{D-m_{n+1}^2}{d_n}
\end{align*}
and
\begin{align*}
m_n\equiv -m_{n+1}\pmod{d_n}.
\end{align*}
The reduced inequalities from the previous step give
\begin{align*}
a_0-d_n<m_n\leq a_0.
\end{align*}
Therefore
\begin{align*}
P(m_{n+1},d_{n+1})=(m_n,d_n)
\end{align*}
for every $n\geq 1$.
We now compare two finite sequences of pairs. For $0\leq i\leq \ell-1$, define
\begin{align*}
R_i:=(m_{\ell-i},d_{\ell-i})
\end{align*}
and
\begin{align*}
S_i:=(m_{i+1},d_i).
\end{align*}
The terminal-pair step gives
\begin{align*}
R_0=(m_\ell,d_\ell)=(a_0,1).
\end{align*}
The initial recurrence gives $m_1=d_0a_0-m_0=a_0$, so
\begin{align*}
S_0=(m_1,d_0)=(a_0,1).
\end{align*}
Thus $R_0=S_0$.
Assume $R_i=S_i$ for some integer $i$ with $0\leq i\leq \ell-2$. Applying $P$ to $R_i$ is legitimate because $R_i=(m_{\ell-i},d_{\ell-i})$ is an actual continued-fraction state. On the reversed side,
\begin{align*}
P(R_i)=P(m_{\ell-i},d_{\ell-i})=(m_{\ell-i-1},d_{\ell-i-1})=R_{i+1}.
\end{align*}
It remains to compute $P(S_i)$ directly, because $S_i=(m_{i+1},d_i)$ is shifted relative to the usual state indexing. Since $R_i=S_i$ and $R_i\in\mathcal{A}$, the pair $S_i$ lies in $\mathcal{A}$. The recurrence gives
\begin{align*}
\delta(m_{i+1},d_i)
=
\frac{D-m_{i+1}^2}{d_i}
=
d_{i+1}.
\end{align*}
Also
\begin{align*}
m_{i+2}=d_{i+1}a_{i+1}-m_{i+1},
\end{align*}
so
\begin{align*}
m_{i+2}\equiv -m_{i+1}\pmod{d_{i+1}}.
\end{align*}
The reduced inequalities applied at index $i+2$ give
\begin{align*}
a_0-d_{i+1}<m_{i+2}\leq a_0.
\end{align*}
By the uniqueness clause in the definition of $P$, these two displayed conditions force
\begin{align*}
P(S_i)=P(m_{i+1},d_i)=(m_{i+2},d_{i+1})=S_{i+1}.
\end{align*}
Since $R_i=S_i$, applying the function $P$ to the same ordered pair gives $P(R_i)=P(S_i)$. Therefore
\begin{align*}
R_{i+1}=S_{i+1}.
\end{align*}
Induction yields
\begin{align*}
(m_{\ell-i},d_{\ell-i})=(m_{i+1},d_i)
\end{align*}
for every integer $i$ with $0\leq i\leq \ell-1$.
[guided]
The purpose of this step is to make the word "reversible" precise. A state is an ordered pair $(M,E)$ for which $E$ divides $D-M^2$; this is why we first define
\begin{align*}
\mathcal{A}:=\{(M,E)\in\mathbb{Z}\times\mathbb{N}: E\mid D-M^2\}.
\end{align*}
For a state $(M,E)\in\mathcal{A}$, the number
\begin{align*}
\delta(M,E)=\frac{D-M^2}{E}
\end{align*}
is the modulus of the predecessor state. The predecessor map $P:\mathcal{A}\to\mathbb{Z}\times\mathbb{N}$ defined by
\begin{align*}
P(M,E):=(r,\delta(M,E))
\end{align*}
chooses the unique integer $r$ satisfying
\begin{align*}
a_0-\delta(M,E)<r\leq a_0,
\qquad
r\equiv -M \pmod{\delta(M,E)}.
\end{align*}
The uniqueness is exactly the elementary residue fact that any interval of length $q\in\mathbb{N}$ contains one representative of each residue class modulo $q$ when the interval is half-open of the form $(b-q,b]$.
For an actual continued-fraction transition, the recurrence gives
\begin{align*}
d_{n+1}=\frac{D-m_{n+1}^2}{d_n},
\qquad
m_n\equiv -m_{n+1}\pmod{d_n}.
\end{align*}
The reduced inequality from the previous step says
\begin{align*}
a_0-d_n<m_n\leq a_0.
\end{align*}
Thus $m_n$ is precisely the representative selected by the definition of $P$, and therefore
\begin{align*}
P(m_{n+1},d_{n+1})=(m_n,d_n)
\end{align*}
for every $n\geq 1$.
Now define
\begin{align*}
R_i:=(m_{\ell-i},d_{\ell-i}),
\qquad
S_i:=(m_{i+1},d_i)
\end{align*}
for $0\leq i\leq \ell-1$. The terminal computation gives
\begin{align*}
R_0=(m_\ell,d_\ell)=(a_0,1),
\end{align*}
and the initial recurrence gives
\begin{align*}
S_0=(m_1,d_0)=(a_0,1).
\end{align*}
So $R_0=S_0$.
Assume $R_i=S_i$ for some $0\leq i\leq \ell-2$. On the reversed side, $R_i$ is an actual state, so the predecessor rule gives
\begin{align*}
P(R_i)=P(m_{\ell-i},d_{\ell-i})=(m_{\ell-i-1},d_{\ell-i-1})=R_{i+1}.
\end{align*}
For the shifted forward side, we must not merely assert the same formula, because $S_i=(m_{i+1},d_i)$ is indexed differently from the usual state $(m_{i+1},d_{i+1})$. We compute it. Since $R_i=S_i$ and $R_i\in\mathcal{A}$, the pair $S_i$ lies in $\mathcal{A}$. Moreover
\begin{align*}
\delta(m_{i+1},d_i)=\frac{D-m_{i+1}^2}{d_i}=d_{i+1}.
\end{align*}
The next recurrence is
\begin{align*}
m_{i+2}=d_{i+1}a_{i+1}-m_{i+1},
\end{align*}
which gives
\begin{align*}
m_{i+2}\equiv -m_{i+1}\pmod{d_{i+1}}.
\end{align*}
The reduced inequalities at index $i+2$ give
\begin{align*}
a_0-d_{i+1}<m_{i+2}\leq a_0.
\end{align*}
Thus $m_{i+2}$ is the unique representative selected by $P$, and hence
\begin{align*}
P(S_i)=P(m_{i+1},d_i)=(m_{i+2},d_{i+1})=S_{i+1}.
\end{align*}
Since $R_i=S_i$ as ordered pairs and $P$ is a function, $P(R_i)=P(S_i)$. Combining the two computations gives
\begin{align*}
R_{i+1}=S_{i+1}.
\end{align*}
Induction proves
\begin{align*}
(m_{\ell-i},d_{\ell-i})=(m_{i+1},d_i)
\end{align*}
for every $0\leq i\leq \ell-1$.
[/guided]
[/step]
[step:Compare the floor formulas in opposite directions]
Fix an integer $i$ with $1\leq i\leq \ell-1$. From the preceding step,
\begin{align*}
m_{\ell-i}=m_{i+1},
\qquad
d_{\ell-i}=d_i.
\end{align*}
Using the floor formula for $a_{\ell-i}$ gives
\begin{align*}
a_{\ell-i}
=
\left\lfloor \frac{a_0+m_{\ell-i}}{d_{\ell-i}}\right\rfloor
=
\left\lfloor \frac{a_0+m_{i+1}}{d_i}\right\rfloor.
\end{align*}
Since
\begin{align*}
m_{i+1}=d_i a_i-m_i,
\end{align*}
we have
\begin{align*}
\frac{a_0+m_{i+1}}{d_i}
=
a_i+\frac{a_0-m_i}{d_i}.
\end{align*}
The reduced inequalities give
\begin{align*}
0\leq \frac{a_0-m_i}{d_i}<1.
\end{align*}
Therefore
\begin{align*}
a_{\ell-i}
=
\left\lfloor a_i+\frac{a_0-m_i}{d_i}\right\rfloor
=
a_i.
\end{align*}
This holds for every $1\leq i\leq \ell-1$, so the period is symmetric before its final term.
[/step]