Usando potencias de matrices.
La formulación matricial de la ecuación de Pell te da una perspectiva muy potente desde el álgebra lineal para generar soluciones. Aunque normalmente las soluciones $(x_n, y_n)$ de la ecuación $x^2 - Dy^2 = 1$ se generan expandiendo la potencia $(x_1 + y_1\sqrt{D})^n$, puedes ver este proceso de forma equivalente como una relación de recurrencia lineal. Si guardas la solución fundamental $(x_1, y_1)$ en una matriz de $2 \times 2$, puedes usar las propiedades de la multiplicación y exponenciación de matrices para analizar la secuencia de soluciones.
Esta técnica es súper importante en las olimpiadas de matemáticas porque conecta la teoría de números con el álgebra lineal. Te permite calcular rapidísimo la $n$-ésima solución usando exponenciación binaria (elevando la matriz al cuadrado repetidamente) y hace que sea más fácil sacar identidades que involucren a $x_n$ y $y_n$. Además, el conjunto de estas matrices forma un subgrupo del grupo lineal especial $SL_2(\mathbb{Z})$, lo que te da una visión muy profunda de su estructura. El determinante de la matriz de transición corresponde a la norma del elemento en el anillo $\mathbb{Z}[\sqrt{D}]$, lo que explica por qué el determinante siempre es $1$ (o $-1$ para la ecuación de Pell negativa) en toda la secuencia.
Imagina que $D$ es un entero positivo que no es un cuadrado perfecto. Si $(x_1, y_1)$ es la solución fundamental de la ecuación de Pell $x^2 - Dy^2 = 1$, definimos la matriz fundamental como: $$ M = \begin{pmatrix} x_1 & Dy_1 \ y_1 & x_1 \end{pmatrix} $$
Puedes encontrar la $n$-ésima solución positiva $(x_n, y_n)$ usando la potencia de la matriz: $$ \begin{pmatrix} x_n \ y_n \end{pmatrix} = \begin{pmatrix} x_1 & Dy_1 \ y_1 & x_1 \end{pmatrix}^{n-1} \begin{pmatrix} x_1 \ y_1 \end{pmatrix} = \begin{pmatrix} x_1 & Dy_1 \ y_1 & x_1 \end{pmatrix}^n \begin{pmatrix} 1 \ 0 \end{pmatrix} $$
De forma equivalente, toda la estructura de la matriz se mantiene al elevarla a una potencia: $$ \begin{pmatrix} x_n & Dy_n \ y_n & x_n \end{pmatrix} = \begin{pmatrix} x_1 & Dy_1 \ y_1 & x_1 \end{pmatrix}^n $$
Para la ecuación de Pell generalizada $x^2 - Dy^2 = N$, si $(u, v)$ es una solución y $(x_n, y_n)$ son soluciones de $x^2 - Dy^2 = 1$, puedes generar nuevas soluciones $(X, Y)$ así: $$ \begin{pmatrix} X \ Y \end{pmatrix} = \begin{pmatrix} x_n & Dy_n \ y_n & x_n \end{pmatrix} \begin{pmatrix} u \ v \end{pmatrix} $$
Teorema: Sea $(x_1, y_1)$ la solución fundamental de $x^2 - Dy^2 = 1$. La solución $(x_n, y_n)$ definida por $x_n + y_n\sqrt{D} = (x_1 + y_1\sqrt{D})^n$ cumple que: $$ \begin{pmatrix} x_n \ y_n \end{pmatrix} = \begin{pmatrix} x_1 & Dy_1 \ y_1 & x_1 \end{pmatrix}^n \begin{pmatrix} 1 \ 0 \end{pmatrix} $$
Demostración: Vamos a usar inducción matemática sobre $n$.
Caso base ($n=1$): $$ \begin{pmatrix} x_1 & Dy_1 \ y_1 & x_1 \end{pmatrix}^1 \begin{pmatrix} 1 \ 0 \end{pmatrix} = \begin{pmatrix} x_1(1) + Dy_1(0) \ y_1(1) + x_1(0) \end{pmatrix} = \begin{pmatrix} x_1 \ y_1 \end{pmatrix} $$ Esto coincide con la solución fundamental.
Paso inductivo: Supón que la fórmula funciona para algún entero $k \geq 1$. Es decir, supón que: $$ \begin{pmatrix} x_k \ y_k \end{pmatrix} = M^k \begin{pmatrix} 1 \ 0 \end{pmatrix} $$ donde $M = \begin{pmatrix} x_1 & Dy_1 \ y_1 & x_1 \end{pmatrix}$.
Por la definición algebraica de las soluciones, tienes: $$ x_{k+1} + y_{k+1}\sqrt{D} = (x_k + y_k\sqrt{D})(x_1 + y_1\sqrt{D}) $$ Si expandes el lado derecho: $$ x_{k+1} + y_{k+1}\sqrt{D} = (x_k x_1 + y_k y_1 D) + (x_k y_1 + y_k x_1)\sqrt{D} $$ Al igualar las partes racionales e irracionales, obtienes estas relaciones de recurrencia:
Ahora, considera el producto de matrices $M^{k+1} \begin{pmatrix} 1 \ 0 \end{pmatrix}$. Por la hipótesis inductiva, sabes que $M^k \begin{pmatrix} 1 \ 0 \end{pmatrix} = \begin{pmatrix} x_k \ y_k \end{pmatrix}$. Por lo tanto: $$ M^{k+1} \begin{pmatrix} 1 \ 0 \end{pmatrix} = M \left( M^k \begin{pmatrix} 1 \ 0 \end{pmatrix} \right) = \begin{pmatrix} x_1 & Dy_1 \ y_1 & x_1 \end{pmatrix} \begin{pmatrix} x_k \ y_k \end{pmatrix} $$ Al hacer la multiplicación de matrices: $$ \begin{pmatrix} x_1 & Dy_1 \ y_1 & x_1 \end{pmatrix} \begin{pmatrix} x_k \ y_k \end{pmatrix} = \begin{pmatrix} x_1 x_k + Dy_1 y_k \ y_1 x_k + x_1 y_k \end{pmatrix} $$ Si comparas este resultado con las relaciones de recurrencia que sacaste arriba, verás que el componente de arriba es exactamente $x_{k+1}$ y el de abajo es exactamente $y_{k+1}$.
Así que $\begin{pmatrix} x_{k+1} \ y_{k+1} \end{pmatrix} = M^{k+1} \begin{pmatrix} 1 \ 0 \end{pmatrix}$. Por el principio de inducción matemática, la fórmula se cumple para todo $n \geq 1$. $\square$