[proofplan]
We use the complete quotients of the continued fraction of $\sqrt{D}$. The key identity expresses the Pell norm $p_n^2-Dq_n^2$ of a convergent in terms of the denominator $d_{n+1}$ of the next complete quotient and the determinant sign of consecutive convergents. At the end of one full period, the complete quotient has denominator $1$, so the norm is exactly $\pm 1$; the parity of the period determines the sign. If the period is odd, one more period gives denominator $1$ again and the determinant sign becomes positive.
[/proofplan]
[step:Define the complete quotients and convergents]
For the continued fraction algorithm applied to $\sqrt{D}$, define integers $m_n,d_n,a_n$ for $n \geq 0$ by
\begin{align*}
m_0 &= 0, &
d_0 &= 1, &
a_0 &= \lfloor \sqrt{D} \rfloor,
\end{align*}
and, for $n \geq 0$,
\begin{align*}
m_{n+1} &= d_n a_n - m_n,\\
d_{n+1} &= \frac{D-m_{n+1}^2}{d_n},\\
a_{n+1} &= \left\lfloor \frac{a_0+m_{n+1}}{d_{n+1}} \right\rfloor.
\end{align*}
The corresponding complete quotient $\alpha_n$ is the real number
\begin{align*}
\alpha_n := \frac{\sqrt{D}+m_n}{d_n}.
\end{align*}
Thus $\alpha_0=\sqrt{D}$, and the continued fraction tail beginning at $a_n$ is $\alpha_n$. More precisely, the recursive identities give
\begin{align*}
\alpha_n=a_n+\frac{1}{\alpha_{n+1}}
\end{align*}
for every $n \geq 0$, because
\begin{align*}
\alpha_n-a_n
&=\frac{\sqrt{D}+m_n-a_nd_n}{d_n}
=\frac{\sqrt{D}-m_{n+1}}{d_n}
=\frac{D-m_{n+1}^2}{d_n(\sqrt{D}+m_{n+1})}
=\frac{d_{n+1}}{\sqrt{D}+m_{n+1}}
=\frac{1}{\alpha_{n+1}}.
\end{align*}
Thus this recursion is exactly the standard complete quotient algorithm for the continued fraction of $\sqrt{D}$, and the recursively defined $a_n$ are the partial quotients in the stated periodic expansion.
Define the convergent numerator and denominator sequences $(p_n)_{n \geq -2}$ and $(q_n)_{n \geq -2}$ by
\begin{align*}
p_{-2} &= 0, &
p_{-1} &= 1, &
p_n &= a_n p_{n-1}+p_{n-2},\\
q_{-2} &= 1, &
q_{-1} &= 0, &
q_n &= a_n q_{n-1}+q_{n-2},
\end{align*}
for every $n \geq 0$. Then the $n$th convergent is $p_n/q_n$.
[/step]
[step:Relate the complete quotient to consecutive convergents]
For every $n \geq 0$, the continued fraction decomposition gives
\begin{align*}
\sqrt{D}
=
\frac{p_n\alpha_{n+1}+p_{n-1}}{q_n\alpha_{n+1}+q_{n-1}}.
\end{align*}
Substituting $\alpha_{n+1}=(\sqrt{D}+m_{n+1})/d_{n+1}$ gives
\begin{align*}
\sqrt{D}
=
\frac{p_n(\sqrt{D}+m_{n+1})+d_{n+1}p_{n-1}}
{q_n(\sqrt{D}+m_{n+1})+d_{n+1}q_{n-1}}.
\end{align*}
Multiplying by the denominator and comparing rational and $\sqrt{D}$ coefficients yields
\begin{align*}
Dq_n &= m_{n+1}p_n+d_{n+1}p_{n-1},\\
p_n &= m_{n+1}q_n+d_{n+1}q_{n-1}.
\end{align*}
[guided]
The complete quotient $\alpha_{n+1}$ is the tail of the continued fraction after the first $n+1$ partial quotients. Reattaching the first $n+1$ partial quotients gives the standard convergent identity
\begin{align*}
\sqrt{D}
=
\frac{p_n\alpha_{n+1}+p_{n-1}}{q_n\alpha_{n+1}+q_{n-1}}.
\end{align*}
Here $p_n/q_n$ is the $n$th convergent and $p_{n-1}/q_{n-1}$ is the preceding convergent.
Now insert the explicit complete quotient
\begin{align*}
\alpha_{n+1}=\frac{\sqrt{D}+m_{n+1}}{d_{n+1}}.
\end{align*}
This gives
\begin{align*}
\sqrt{D}
=
\frac{p_n(\sqrt{D}+m_{n+1})+d_{n+1}p_{n-1}}
{q_n(\sqrt{D}+m_{n+1})+d_{n+1}q_{n-1}}.
\end{align*}
Multiplying by the denominator gives
\begin{align*}
\sqrt{D}\bigl(q_n(\sqrt{D}+m_{n+1})+d_{n+1}q_{n-1}\bigr)
=
p_n(\sqrt{D}+m_{n+1})+d_{n+1}p_{n-1}.
\end{align*}
Expanding both sides,
\begin{align*}
Dq_n+\sqrt{D}\,(m_{n+1}q_n+d_{n+1}q_{n-1})
=
m_{n+1}p_n+d_{n+1}p_{n-1}+\sqrt{D}\,p_n.
\end{align*}
Since $D$ is not a perfect square, $\sqrt{D}$ is irrational, so equality of two expressions of the form $r+s\sqrt{D}$ forces equality of the rational coefficients and equality of the $\sqrt{D}$ coefficients. Therefore
\begin{align*}
Dq_n &= m_{n+1}p_n+d_{n+1}p_{n-1},\\
p_n &= m_{n+1}q_n+d_{n+1}q_{n-1}.
\end{align*}
[/guided]
[/step]
[step:Compute the Pell norm of each convergent]
Using the two identities from the previous step, we compute
\begin{align*}
p_n^2-Dq_n^2
&=
p_n(m_{n+1}q_n+d_{n+1}q_{n-1})
-
q_n(m_{n+1}p_n+d_{n+1}p_{n-1})\\
&=
d_{n+1}(p_nq_{n-1}-p_{n-1}q_n).
\end{align*}
The convergent determinant identity gives
\begin{align*}
p_nq_{n-1}-p_{n-1}q_n=(-1)^{n-1}.
\end{align*}
Hence, for every $n \geq 0$,
\begin{align*}
p_n^2-Dq_n^2 = (-1)^{n-1}d_{n+1}.
\end{align*}
[guided]
The goal is to convert the continued-fraction data into the Pell expression $p_n^2-Dq_n^2$. From the previous step we have
\begin{align*}
p_n &= m_{n+1}q_n+d_{n+1}q_{n-1},\\
Dq_n &= m_{n+1}p_n+d_{n+1}p_{n-1}.
\end{align*}
Substitute the first identity into the factor $p_n^2$ and the second identity into the factor $Dq_n^2$:
\begin{align*}
p_n^2-Dq_n^2
&=
p_n(m_{n+1}q_n+d_{n+1}q_{n-1})
-
q_n(m_{n+1}p_n+d_{n+1}p_{n-1})\\
&=
m_{n+1}p_nq_n+d_{n+1}p_nq_{n-1}
-
m_{n+1}p_nq_n-d_{n+1}p_{n-1}q_n\\
&=
d_{n+1}(p_nq_{n-1}-p_{n-1}q_n).
\end{align*}
The terms containing $m_{n+1}$ cancel. This cancellation is the reason complete quotients are useful here: the remaining factor is exactly the determinant of two consecutive convergents.
The determinant identity for convergents is
\begin{align*}
p_nq_{n-1}-p_{n-1}q_n=(-1)^{n-1}.
\end{align*}
It follows directly from the recurrences for $p_n$ and $q_n$ by induction, with the initial values $p_{-1}=1$, $q_{-1}=0$, $p_0=a_0$, and $q_0=1$. Therefore
\begin{align*}
p_n^2-Dq_n^2 = (-1)^{n-1}d_{n+1}.
\end{align*}
[/guided]
[/step]
[step:Evaluate the norm after one full period]
We first justify the complete quotient value at the end of one period. Since
\begin{align*}
\sqrt{D}=[a_0;\overline{a_1,\dots,a_\ell}],
\end{align*}
the tail beginning after one full period is
\begin{align*}
\alpha_\ell=[a_\ell;\overline{a_1,\dots,a_\ell}].
\end{align*}
For the standard continued fraction of a quadratic irrational $\sqrt{D}$ with $D$ positive and not a square, the last partial quotient in the minimal period is $a_\ell=2a_0$, and the complete quotient at this position is
\begin{align*}
\alpha_\ell=\sqrt{D}+a_0.
\end{align*}
Equivalently, comparing this with the defining form
\begin{align*}
\alpha_\ell=\frac{\sqrt{D}+m_\ell}{d_\ell}
\end{align*}
gives
\begin{align*}
m_\ell=a_0,
\qquad
d_\ell=1.
\end{align*}
Taking $n=\ell-1$ in the norm identity gives
\begin{align*}
p_{\ell-1}^2-Dq_{\ell-1}^2
=
(-1)^{\ell-2}d_\ell
=
(-1)^\ell.
\end{align*}
Thus, if $\ell$ is even,
\begin{align*}
p_{\ell-1}^2-Dq_{\ell-1}^2=1.
\end{align*}
[guided]
The only new input in this step is the endpoint behaviour of the complete quotients for square-root continued fractions. The hypotheses needed for that fact are exactly the present ones: $D \in \mathbb{N}$ is positive and not a square, so $\sqrt{D}$ is a quadratic irrational greater than $0$, and its continued fraction is purely periodic after the initial term:
\begin{align*}
\sqrt{D}=[a_0;\overline{a_1,\dots,a_\ell}].
\end{align*}
For such a square-root continued fraction, the final partial quotient in the minimal period is $a_\ell=2a_0$, and the complete quotient after one period is
\begin{align*}
\alpha_\ell=\sqrt{D}+a_0.
\end{align*}
This statement is the standard endpoint property of the continued fraction algorithm for nonsquare quadratic surds.
Now compare this endpoint value with the complete quotient form already defined in the proof:
\begin{align*}
\alpha_\ell=\frac{\sqrt{D}+m_\ell}{d_\ell}.
\end{align*}
Since $d_\ell$ is a positive integer and $m_\ell$ is an integer, equality with $\sqrt{D}+a_0$ forces the coefficient of $\sqrt{D}$ to be $1$ and the rational part to be $a_0$. Hence
\begin{align*}
m_\ell=a_0,
\qquad
d_\ell=1.
\end{align*}
Substituting $n=\ell-1$ into the norm identity from the previous step gives
\begin{align*}
p_{\ell-1}^2-Dq_{\ell-1}^2
=
(-1)^{\ell-2}d_\ell
=
(-1)^\ell.
\end{align*}
Therefore, when $\ell$ is even,
\begin{align*}
p_{\ell-1}^2-Dq_{\ell-1}^2=1.
\end{align*}
[/guided]
[/step]
[step:Pass to two periods when the period length is odd]
Assume now that $\ell$ is odd. The periodic expansion implies that the complete quotient tail after two periods is the same as the tail after one period:
\begin{align*}
\alpha_{2\ell}=[a_\ell;\overline{a_1,\dots,a_\ell}]=\alpha_\ell.
\end{align*}
Since the complete quotient representation $\alpha_n=(\sqrt{D}+m_n)/d_n$ is uniquely determined by the complete quotient algorithm, the equality $\alpha_{2\ell}=\alpha_\ell$ gives
\begin{align*}
m_{2\ell}=m_\ell=a_0,
\qquad
d_{2\ell}=d_\ell=1.
\end{align*}
Taking $n=2\ell-1$ in the norm identity gives
\begin{align*}
p_{2\ell-1}^2-Dq_{2\ell-1}^2
=
(-1)^{2\ell-2}d_{2\ell}
=
1.
\end{align*}
Therefore $(p_{2\ell-1},q_{2\ell-1})$ solves
\begin{align*}
x^2-Dy^2=1.
\end{align*}
Together with the even-period case, this proves the theorem.
[guided]
When $\ell$ is odd, the one-period computation gives the negative Pell norm, so we pass through a second period. The continued fraction is periodic after the initial term, so the tail after two full periods is the same continued fraction tail as the tail after one full period:
\begin{align*}
\alpha_{2\ell}=[a_\ell;\overline{a_1,\dots,a_\ell}]=\alpha_\ell.
\end{align*}
The complete quotient algorithm determines the pair $(m_n,d_n)$ from the tail $\alpha_n=(\sqrt{D}+m_n)/d_n$, with $d_n>0$. Therefore equality of the two complete quotients gives equality of their integer parameters:
\begin{align*}
m_{2\ell}=m_\ell,
\qquad
d_{2\ell}=d_\ell.
\end{align*}
From the one-period endpoint already proved, $m_\ell=a_0$ and $d_\ell=1$, so
\begin{align*}
m_{2\ell}=a_0,
\qquad
d_{2\ell}=1.
\end{align*}
Now substitute $n=2\ell-1$ into the norm identity:
\begin{align*}
p_{2\ell-1}^2-Dq_{2\ell-1}^2
=
(-1)^{2\ell-2}d_{2\ell}
=
1.
\end{align*}
Thus $(p_{2\ell-1},q_{2\ell-1})$ satisfies the Pell equation
\begin{align*}
x^2-Dy^2=1.
\end{align*}
Combining this with the even-period case proves both alternatives stated in the theorem.
[/guided]
[/step]