Уравнение
x^2 − 5y^2 = 1
— это классическое уравнение Пелля с D = 5. Поэтому естественно применять теорию уравнений Пелля (эквивалентно — факторизацию в квадратичном поле Q(√5) и использование группы единиц в кольце целых соответствующего поля). Этот подход даёт и простой способ подобрать наименьшее нетривиальное решение, и полное описание всех целых решений.

1) Нахождение фундаментального (наименьшего положительного) решения.
Проверкой малых y видим, что y = 4 даёт x^2 = 1 + 5·4^2 = 81, т.е. x = 9. Значит (x1,y1) = (9,4) — нетривиальное положительное решение. (Ещё его можно получить из непрерывной дроби √5 = [2;4,4,4,…], ближайший несокращённый дробный приблизитель 9/4 даёт x=9, y=4.)

2) Общая форма решений (стандартный вывод для уравнения Пелля).
Рассмотрим элемент 9 + 4√5 ∈ Z[√5]. Его нормa N(9 + 4√5) = (9 + 4√5)(9 − 4√5) = 1, т.е. это единица по норме. Тогда все положительные решения уравнения задаются степенями этой единицы:
x_n + y_n√5 = (9 + 4√5)^n, n = 0,1,2,….
Соответственно пары (x_n,y_n) — целые и дают x_n^2 − 5 y_n^2 = 1. Примеры: n=0 → (1,0); n=1 → (9,4); n=2 → (161,72) и т.д.

3) Рекуррентная и явная формулы.
Умножение (x+y√5)(9+4√5) даёт рекуррент:
x_{n+1} = 9 x_n + 20 yn,
y{n+1} = 4 x_n + 9 y_n,
с начальными x_0 = 1, y_0 = 0 (или x_1 = 9, y_1 = 4).

Явная формула через сопряжённые корни:
x_n = ((9+4√5)^n + (9−4√5)^n)/2,
y_n = ((9+4√5)^n − (9−4√5)^n)/(2√5).

4) Полнота решения (обоснование, что других целых решений нет).
Теорема об уравнении Пелля утверждает: если D не является полным квадратом, то все целые решения x^2 − D y^2 = 1 (в положительных x,y) получаются из фундаментальной положительной пары (x1,y1) как (x_n,y_n), где x_n + y_n√D = (x1 + y1√D)^n. Для D = 5 мы нашли фундаментальное (9,4), следовательно все целые решения с положительными x,y — именно перечисленные выше. Отрицательные значения x или y допускаются, так как уравнение содержит только квадраты: если (x,y) — решение, то любая комбинация (±x, ±y) с произвольными знаками по модулю сохраняет уравнение; поэтому полное множество целых пар — все пары (±x_n, ±y_n) (с теми же n).

Можно также пояснить с точки зрения теории чисел: в кольце целых квадратичного поля Q(√5) группа единиц бесконечна циклическая; фундаментальная единица, дающая целые a + b√5, равна 9 + 4√5, и её степени дают все единицы нормы 1, а потому все решения уравнения.

Итого: все целые решения уравнения x^2 − 5y^2 = 1 заданы формулой
x + y√5 = ±(9 + 4√5)^n, n ∈ Z,
или эквивалентно (для n ≥ 0)
(x,y) = (±x_n, ±y_n), где x_n,y_n определены рекурсией выше (или явными формулами).