[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*}[/step]

[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: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$.[/step]

[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: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]