Мотивация и контекст

Рассмотрим задачу: найти целые $x, y > 0$ такие, что $x^2 - 2y^2 = 1$ и $x + y$ максимально при $x + y \le 10^{18}$. Это уравнение Пелля — одна из центральных задач теории чисел. Его решения связаны с наилучшими рациональными приближениями к $\sqrt{2}$

​ — именно здесь в игру вступают непрерывные дроби.

Непрерывные дроби — это представление чисел в виде $a_0 + \cfrac{1}{a_1 + \cfrac{1}{a_2 + \cfrac{1}{\ddots}}}$. Этот инструмент возник ещё в работах Евклида (алгоритм НОД — это и есть разложение в непрерывную дробь) и сегодня используется в задачах:

• Наилучшие рациональные приближения: дроби-подходящие — лучшие приближения к иррациональному числу при заданной ограниченности знаменателя.
• Уравнение Пелля: фундаментальное решение находится через период непрерывной дроби $\sqrt{D}$.
• Криптоанализ: атака Вистра на RSA использует непрерывные дроби для нахождения малого закрытого показателя.
• Задачи на точные вычисления: представление сравнимых иррациональных через подходящие дроби.
• Задачи на цепные рекуррентности: нахождение рациональной точки на кривой, ближайшей к заданной.

Теория

Конечные непрерывные дроби

Определение. Конечной непрерывной дробью называется выражение вида: $[a_0; a_1, a_2, \ldots, a_n] = a_0 + \cfrac{1}{a_1 + \cfrac{1}{a_2 + \cfrac{1}{\ddots + \cfrac{1}{a_n}}}}$ где $a_0 \in \mathbb{Z}$, $a_i \in \mathbb{Z}_{>0}$ при $i \ge 1$.

Теорема. Каждое рациональное число имеет ровно два представления в виде конечной непрерывной дроби: $[\ldots, a_n]$ с $a_n \ge 2$ и $[\ldots, a_n - 1, 1]$.

Алгоритм разложения. Для рационального $r = p/q$ ($q > 0$):

1. $a_0 = \lfloor r \rfloor$, $r_0 = r - a_0 = \{r\}$ (дробная часть)
2. Если $r_0 = 0$, конец.
3. $r_1 = 1/r_0$, $a_1 = \lfloor r_1 \rfloor$, $r_1 = r_1 - a_1$, ...

По существу, это алгоритм Евклида: $p = a_0 q + r_0, \quad q = a_1 r_0 + r_1, \quad r_0 = a_2 r_1 + r_2, \ldots$

Пример: $\frac{355}{113} = [3; 7, 16] = 3 + \frac{1}{7 + \frac{1}{16}}$. Проверка: $7 + \frac{1}{16} = \frac{113}{16}$, $3 + \frac{16}{113} = \frac{355}{113}$. ✓

Подходящие дроби

Определение. $k$-я подходящая дробь (конвергент) числа $[a_0; a_1, a_2, \ldots]$ определяется как: $\frac{p_k}{q_k} = [a_0; a_1, \ldots, a_k]$

Рекуррентные формулы. $p_{-1} = 1, \quad p_0 = a_0, \quad p_k = a_k p_{k-1} + p_{k-2}$ $q_{-1} = 0, \quad q_0 = 1, \quad q_k = a_k q_{k-1} + q_{k-2}$

Доказательство. Индукция по $k$. База: $p_0/q_0 = a_0/1$. Шаг: $[a_0; a_1, \ldots, a_k] = [a_0; a_1, \ldots, a_{k-2}, a_{k-1} + 1/a_k]$. Применяя рекуррентное соотношение с $a_{k-1}$ заменённым на $a_{k-1} + 1/a_k$: $[a_0; \ldots, a_{k-1} + 1/a_k] = \frac{(a_{k-1} + 1/a_k) p_{k-2} + p_{k-3}}{(a_{k-1} + 1/a_k) q_{k-2} + q_{k-3}} = \frac{a_k p_{k-1} + p_{k-2}}{a_k q_{k-1} + q_{k-2}} = \frac{p_k}{q_k} \square$

Теорема (основные свойства подходящих дробей).

1. $p_k q_{k-1} - p_{k-1} q_k = (-1)^{k-1}$ (тождество Кэши).
2. $\gcd(p_k, q_k) = 1$ (числитель и знаменатель взаимно просты).
3. Подходящие дроби чётного номера монотонно возрастают, нечётного — убывают.
4. Все подходящие дроби «зажимают» значение $x$: $p_0/q_0 < p_2/q_2 < \ldots < x < \ldots < p_3/q_3 < p_1/q_1$.

Доказательство (п.1). Индукция: при $k=0$: $p_0 q_{-1} - p_{-1} q_0 = a_0 \cdot 0 - 1 \cdot 1 = -1 = (-1)^{-1}$. При переходе $k-1 \to k$: $p_k q_{k-1} - p_{k-1} q_k = (a_k p_{k-1} + p_{k-2}) q_{k-1} - p_{k-1}(a_k q_{k-1} + q_{k-2})$ $= p_{k-2} q_{k-1} - p_{k-1} q_{k-2} = -(p_{k-1} q_{k-2} - p_{k-2} q_{k-1}) = -(-1)^{k-2} = (-1)^{k-1} \square$

Теорема (наилучшее приближение). Пусть $x = [a_0; a_1, \ldots]$ — иррациональное число. Подходящая дробь $p_k/q_k$ является наилучшим приближением к $x$ из всех дробей со знаменателем $\le q_k$: для любого рационального $p/q$ с $0 < q \le q_k$ и $p/q \ne p_k/q_k$: $|x - p_k/q_k| < |x - p/q|$

Бесконечные непрерывные дроби

Теорема. Каждое иррациональное число представляется единственной бесконечной непрерывной дробью.

Алгоритм разложения иррационального числа. Для $\alpha_0 = x$: $a_k = \lfloor \alpha_k \rfloor, \quad \alpha_{k+1} = \frac{1}{\alpha_k - a_k}$

Периодические непрерывные дроби и квадратичные иррациональности

Теорема Лагранжа. Непрерывная дробь периодична тогда и только тогда, когда число является квадратичной иррациональностью, то есть имеет вид $\frac{a + b\sqrt{D}}{c}$

​​ с $a, b, c \in \mathbb{Z}$, $D > 0$ не являющимся полным квадратом.

Структура разложения $\sqrt{D}$

​. Разложение $\sqrt{D}$ имеет вид $[a_0; \overline{a_1, a_2, \ldots, a_{2l-1}, 2a_0}]$, где $a_0 = \lfloor\sqrt{D}\rfloor$ и период начинается с $a_1$.

Алгоритм разложения $\sqrt{D}$

​. Поддерживаем представление $\alpha_k = \frac{a_0 \cdot m_k + \sqrt{D}}{d_k}$ ... точнее, $\alpha_k = \frac{\sqrt{D} + m_k}{d_k}$:

• $m_0 = 0$, $d_0 = 1$, $a_0 = \lfloor\sqrt{D}\rfloor$
• $m_{k+1} = d_k \cdot a_k - m_k$
• $d_{k+1} = (D - m_{k+1}^2) / d_k$
• $a_{k+1} = \lfloor(a_0 + m_{k+1}) / d_{k+1}\rfloor$

Период заканчивается, когда $a_k = 2a_0$.

Лемма. Все $m_k, d_k$ — целые числа, $d_k \mid D - m_k^2$.

Доказательство. Индукция: при $k=0$: $d_0 = 1 \mid D - 0$. Шаг: $d_{k+1} = (D - m_{k+1}^2)/d_k = (D - (d_k a_k - m_k)^2)/d_k$. Раскрывая: $(D - m_k^2)/d_k - 2a_k m_k + a_k^2 d_k$. Первое слагаемое целое по индукционному предположению, остальные целые. $\square$

Уравнение Пелля

Уравнение Пелля: $x^2 - D y^2 = 1$, где $D > 0$ — натуральное число, не являющееся полным квадратом.

Теорема (Лагранж). Уравнение Пелля всегда имеет решение в натуральных числах. Если период разложения $\sqrt{D}$

​ в непрерывную дробь имеет длину $l$, то:

• При чётном $l$: фундаментальное решение $(x_1, y_1) = (p_{l-1}, q_{l-1})$.
• При нечётном $l$: фундаментальное решение $(x_1, y_1) = (p_{2l-1}, q_{2l-1})$.

Здесь $p_k/q_k$ — подходящие дроби разложения $\sqrt{D}$

​.

Доказательство (схема). Из свойств подходящих дробей: $p_{k-1}^2 - D q_{k-1}^2 = (-1)^k d_k$ (доказывается индукцией через рекуррентные соотношения). Когда $a_k = 2a_0$ (конец периода), $d_k = 1$, и $p_{k-1}^2 - D q_{k-1}^2 = (-1)^k$. При чётном $l$ (длина периода): $k = l$, $(-1)^l = 1$, получаем первое решение. При нечётном $l$: берём $k = 2l$, $(-1)^{2l} = 1$. $\square$

Генерация всех решений. Если $(x_1, y_1)$ — фундаментальное решение, то все решения даются формулой: $x_n + y_n\sqrt{D} = (x_1 + y_1\sqrt{D})^n, \quad n = 1, 2, 3, \ldots$

​=(x1​+y1​D

​)n,n=1,2,3,…

Или рекуррентно: $x_{n+1} = x_1 x_n + D y_1 y_n, \quad y_{n+1} = x_1 y_n + y_1 x_n$

Доказательство. Если $x^2 - Dy^2 = 1$ и $u^2 - Dv^2 = 1$, то $(xu + Dvy)^2 - D(xv + yu)^2 = (x^2 - Dy^2)(u^2 - Dv^2) = 1$. Значит произведение двух решений (в кольце $\mathbb{Z}[\sqrt{D}]$

​]) тоже является решением.

Чтобы убедиться, что формула даёт все решения, используется факт: все решения образуют циклическую группу, порождённую фундаментальным решением $(x_1, y_1)$.

Связанное уравнение: $x^2 - Dy^2 = -1$ имеет решение тогда и только тогда, когда период разложения $\sqrt{D}$

​ нечётный. Минимальное решение: $(p_{l-1}, q_{l-1})$.

Обобщённое уравнение Пелля: $x^2 - Dy^2 = N$ при $|N| < \sqrt{D}$

​ решается через подходящие дроби. Решение $(x, y)$ существует тогда и только тогда, когда $N = p_k^2 - D q_k^2$ для некоторого $k$.

Ключевые формулы (сводная таблица)

Реализация на C++20 (с анализом сложности)

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef __int128 lll;

// ─── Структура для работы с большими числами через __int128 ────────────────
// Для задач Пелля числа могут быть очень большими — используем __int128
// при D до 10^6, или Python-style big int при бОльших значениях.

// ─── Разложение рационального p/q в непрерывную дробь ─────────────────────
// Возвращает коэффициенты [a0; a1, ..., an]
// Сложность: O(log(max(p, q))) — как алгоритм Евклида
vector<ll> rational_cf(ll p, ll q) {
    vector<ll> cf;
    while (q != 0) {
        cf.push_back(p / q);
        p %= q;
        swap(p, q);
    }
    return cf;
}

// ─── Вычисление подходящих дробей по коэффициентам ─────────────────────────
// Возвращает вектор пар (p_k, q_k)
// Сложность: O(n)
vector<pair<ll,ll>> convergents(const vector<ll>& cf) {
    vector<pair<ll,ll>> conv;
    ll p_prev = 1, p_curr = cf[0];
    ll q_prev = 0, q_curr = 1;
    conv.push_back({p_curr, q_curr});
    for (int i = 1; i < (int)cf.size(); i++) {
        ll p_next = cf[i] * p_curr + p_prev;
        ll q_next = cf[i] * q_curr + q_prev;
        conv.push_back({p_next, q_next});
        p_prev = p_curr; p_curr = p_next;
        q_prev = q_curr; q_curr = q_next;
    }
    return conv;
}

// ─── Разложение sqrt(D) в непрерывную дробь ───────────────────────────────
// Возвращает (a0, [a1, a2, ..., period])
// Сложность: O(sqrt(D)) — длина периода O(sqrt(D))
struct SqrtCF {
    ll a0;           // целая часть sqrt(D)
    vector<ll> period; // периодическая часть

    SqrtCF(ll D) {
        a0 = (ll)sqrt((double)D);
        while ((a0+1)*(a0+1) <= D) a0++; // корректируем из-за погрешности
        while (a0*a0 > D) a0--;

        if (a0 * a0 == D) return; // D — полный квадрат

        ll m = 0, d = 1, a = a0;
        do {
            m = d * a - m;
            d = (D - m * m) / d;
            a = (a0 + m) / d;
            period.push_back(a);
        } while (a != 2 * a0);
    }

    // Получить k-й коэффициент (0-индексированный среди всех)
    // a_0 = a0, a_k = period[(k-1) % period.size()] при k >= 1
    ll get(ll k) const {
        if (k == 0) return a0;
        if (period.empty()) return -1; // D — квадрат
        return period[(k - 1) % period.size()];
    }
};

// ─── Решение уравнения Пелля x^2 - D*y^2 = 1 ─────────────────────────────
// Возвращает фундаментальное решение (x1, y1)
// При D — полный квадрат, решений нет (кроме x=1, y=0)
// Сложность: O(sqrt(D)) — поиск периода
pair<lll, lll> pell_fundamental(ll D) {
    ll a0 = (ll)sqrt((double)D);
    while ((a0+1)*(a0+1) <= D) a0++;
    while (a0*a0 > D) a0--;

    if (a0 * a0 == D) return {1, 0}; // D — квадрат

    // Разложение sqrt(D)
    ll m = 0, d = 1, a = a0;
    lll p_prev = 1, p_curr = a0;
    lll q_prev = 0, q_curr = 1;

    for (int iter = 0; iter < 1000000; iter++) {
        m = d * a - m;
        d = (D - m * m) / d;
        a = (a0 + m) / d;

        lll p_next = a * p_curr + p_prev;
        lll q_next = a * q_curr + q_prev;
        p_prev = p_curr; p_curr = p_next;
        q_prev = q_curr; q_curr = q_next;

        // Проверяем x^2 - D*y^2 = 1
        if (p_prev * p_prev - (lll)D * q_prev * q_prev == 1) {
            return {p_prev, q_prev};
        }
        // Конец периода: a == 2*a0
        // Следующая итерация продолжит до второго конца периода если нужно
    }
    return {-1, -1}; // не должно происходить
}

// ─── Структура для уравнения Пелля (все решения) ──────────────────────────
struct PellSolver {
    ll D;
    lll x1, y1; // фундаментальное решение

    PellSolver(ll D) : D(D) {
        auto [fx, fy] = pell_fundamental(D);
        x1 = fx; y1 = fy;
    }

    // Генерация n-го решения
    pair<lll, lll> nth(int n) const {
        lll xn = 1, yn = 0;
        lll base_x = x1, base_y = y1;
        // Быстрое возведение в степень в кольце Z[sqrt(D)]
        int e = n;
        while (e > 0) {
            if (e & 1) {
                lll new_x = xn * base_x + (lll)D * yn * base_y;
                lll new_y = xn * base_y + yn * base_x;
                xn = new_x; yn = new_y;
            }
            lll new_bx = base_x * base_x + (lll)D * base_y * base_y;
            lll new_by = 2 * base_x * base_y;
            base_x = new_bx; base_y = new_by;
            e >>= 1;
        }
        return {xn, yn};
    }

    // Итератор решений: генерируем (x_{n+1}, y_{n+1}) из (x_n, y_n)
    pair<lll, lll> next_solution(lll xn, lll yn) const {
        return {x1 * xn + (lll)D * y1 * yn, x1 * yn + y1 * xn};
    }
};

// ─── Наилучшие приближения: найти все (p, q) с q ≤ Q, p/q ≈ x ─────────────
// Возвращает вектор подходящих дробей с q ≤ Q
// Для x = p0/q0 (рациональное)
vector<pair<ll,ll>> best_approximations(ll p0, ll q0, ll Q) {
    auto cf = rational_cf(p0, q0);
    auto conv = convergents(cf);
    vector<pair<ll,ll>> result;
    for (auto [p, q] : conv) {
        if (q <= Q) result.push_back({p, q});
    }
    return result;
}

// ─── Демонстрация ──────────────────────────────────────────────────────────
int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    // 1. Разложение рационального числа
    {
        ll p = 355, q = 113;
        auto cf = rational_cf(p, q);
        cout << p << "/" << q << " = [";
        for (int i = 0; i < (int)cf.size(); i++) {
            if (i) cout << "; ";
            cout << cf[i];
        }
        cout << "]\n";
        auto conv = convergents(cf);
        for (auto [pp, qq] : conv)
            cout << "  " << pp << "/" << qq << "\n";
    }

    // 2. Разложение sqrt(D)
    {
        for (ll D : {2, 3, 5, 6, 7, 13, 23}) {
            SqrtCF scf(D);
            cout << "sqrt(" << D << ") = [" << scf.a0 << "; (";
            for (int i = 0; i < (int)scf.period.size(); i++) {
                if (i) cout << ",";
                cout << scf.period[i];
            }
            cout << ")]\n";
        }
    }

    // 3. Уравнение Пелля
    {
        for (ll D : {2, 3, 5, 6, 7, 13, 23, 61}) {
            PellSolver ps(D);
            // Выводим через __int128 (приходится через строку)
            auto print128 = [](lll n) {
                if (n < 0) { cout << "-"; n = -n; }
                string s;
                if (n == 0) { cout << "0"; return; }
                while (n > 0) { s += ('0' + (int)(n%10)); n /= 10; }
                reverse(s.begin(), s.end());
                cout << s;
            };
            cout << "x^2 - " << D << "*y^2 = 1: x1=";
            print128(ps.x1); cout << ", y1="; print128(ps.y1);
            cout << ", verify: "; print128(ps.x1*ps.x1 - (lll)D*ps.y1*ps.y1);
            cout << "\n";
        }
    }

    // 4. Генерация нескольких решений уравнения Пелля для D=2
    {
        PellSolver ps(2);
        cout << "\nSolutions of x^2 - 2y^2 = 1:\n";
        lll x = 1, y = 0;
        for (int i = 0; i < 8; i++) {
            auto [nx, ny] = ps.next_solution(x, y);
            x = nx; y = ny;
            cout << "  (" << (ll)x << ", " << (ll)y << ")\n";
        }
    }

    return 0;
}

Анализ сложности:

Разбор задачи 1 (простая)

Задача. Дано натуральное число $n$. Найти дробь $p/q$ ($q \le n$), ближайшую к $\sqrt{2}$

​.

Решение. Используем подходящие дроби разложения $\sqrt{2} = [1; \overline{2}]$

​=[1;2].

$\sqrt{2} = [1; 2, 2, 2, \ldots]$

​=[1;2,2,2,…]

Подходящие дроби: $1/1, 3/2, 7/5, 17/12, 41/29, 99/70, 239/169, \ldots$

Генерируем подходящие дроби, пока знаменатель $\le n$, и берём последнюю.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

int main(){
    ll n; cin >> n;
    // sqrt(2) = [1; 2, 2, 2, ...]
    ll p_prev = 1, p_curr = 1; // p_{-1}=1, p_0=1
    ll q_prev = 0, q_curr = 1; // q_{-1}=0, q_0=1
    // Первый коэффициент a_0 = 1
    // Все последующие a_k = 2

    while (true) {
        ll p_next = 2 * p_curr + p_prev;
        ll q_next = 2 * q_curr + q_prev;
        if (q_next > n) break;
        p_prev = p_curr; p_curr = p_next;
        q_prev = q_curr; q_curr = q_next;
    }
    cout << p_curr << "/" << q_curr << "\n";
    // Погрешность: |sqrt(2) - p/q| < 1/(q * q_next)
    return 0;
}

Пример: $n = 100$. Последняя подходящая дробь с $q \le 100$: $99/70$ (т.к. $70 \le 100 < 169$). Погрешность: $|99/70 - \sqrt{2}| \approx 5 \times 10^{-6}$

​∣≈5×10−6.

Разбор задачи 2 (средняя)

Задача. Найти наименьшее натуральное решение уравнения $x^2 - 13y^2 = 1$.

Решение. Применяем теорему Лагранжа.

Шаг 1: Разложение $\sqrt{13}$

​:

• $a_0 = 3$ (т.к. $3^2 = 9 < 13 < 16 = 4^2$)
• $m_1 = 3$, $d_1 = 13 - 9 = 4$, $a_1 = \lfloor(3+3)/4\rfloor = 1$
• $m_2 = 4 - 3 = 1$, $d_2 = (13 - 1)/4 = 3$, $a_2 = \lfloor(3+1)/3\rfloor = 1$
• $m_3 = 3 - 1 = 2$, $d_3 = (13 - 4)/3 = 3$, $a_3 = \lfloor(3+2)/3\rfloor = 1$
• $m_4 = 3 - 2 = 1$, $d_4 = (13 - 1)/3 = 4$, $a_4 = \lfloor(3+1)/4\rfloor = 1$
• $m_5 = 4 - 1 = 3$, $d_5 = (13 - 9)/4 = 1$, $a_5 = \lfloor(3+3)/1\rfloor = 6 = 2 \cdot 3$

Период: $[1, 1, 1, 1, 6]$, длина $l = 5$ (нечётная).

Подходящие дроби:

• $p_0/q_0 = 3/1$
• $p_1/q_1 = 1 \cdot 3 + 1 / 1 \cdot 1 + 0 = 4/1$
• $p_2/q_2 = 1 \cdot 4 + 3 / 1 \cdot 1 + 1 = 7/2$
• $p_3/q_3 = 1 \cdot 7 + 4 / 1 \cdot 2 + 1 = 11/3$
• $p_4/q_4 = 1 \cdot 11 + 7 / 1 \cdot 3 + 2 = 18/5$

Нечётный период: фундаментальное решение — $p_{2l-1} = p_9$... но используем прямую проверку:

При нечётном $l = 5$: $p_{l-1} = p_4 = 18$, $q_4 = 5$. Проверяем: $18^2 - 13 \cdot 25 = 324 - 325 = -1$. Значит $(18, 5)$ решает $x^2 - 13y^2 = -1$.

Для уравнения $= 1$: берём $k = 2l = 10$-й конвергент, или применяем $(x_1, y_1) = (x_0^2 + 13 y_0^2, 2 x_0 y_0) = (18^2 + 13 \cdot 25, 2 \cdot 18 \cdot 5) = (324 + 325, 180) = (649, 180)$.

Проверка: $649^2 - 13 \cdot 180^2 = 421201 - 421200 = 1$. ✓

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128 lll;

void print128(lll n) {
    if (n == 0) { cout << "0"; return; }
    string s; while (n > 0) { s += ('0' + (int)(n%10)); n /= 10; }
    reverse(s.begin(), s.end()); cout << s;
}

int main(){
    ll D; cin >> D;
    ll a0 = (ll)sqrt((double)D);
    while ((a0+1)*(a0+1) <= D) a0++;
    while (a0*a0 > D) a0--;

    if (a0*a0 == D) { cout << "D is a perfect square\n"; return 0; }

    ll m = 0, d = 1, a = a0;
    lll p_prev = 1, p_curr = a0;
    lll q_prev = 0, q_curr = 1;

    for (int iter = 0; iter < 2000000; iter++) {
        m = d * a - m;
        d = (D - m * m) / d;
        a = (a0 + m) / d;

        lll p_next = a * p_curr + p_prev;
        lll q_next = a * q_curr + q_prev;
        p_prev = p_curr; p_curr = p_next;
        q_prev = q_curr; q_curr = q_next;

        lll check = p_prev * p_prev - (lll)D * q_prev * q_prev;
        if (check == 1) {
            cout << "x1 = "; print128(p_prev);
            cout << ", y1 = "; print128(q_prev);
            cout << "\nVerify: x^2 - " << D << "*y^2 = "; print128(check);
            cout << "\n";
            return 0;
        }
        if (check == -1) {
            cout << "x^2 - " << D << "y^2 = -1 has solution: ";
            cout << "("; print128(p_prev); cout << ", "; print128(q_prev); cout << ")\n";
            // Продолжаем ещё один период для = 1
        }
    }
    return 0;
}

Разбор задачи 3 (сложная — Div.1 D / ICPC WF)

Задача. Дано $D$ и $N$ ($1 \le D \le 10^6$, $N \le 10^{18}$). Найти $n$-е решение $(x_n, y_n)$ уравнения $x^2 - Dy^2 = 1$ по модулю $10^9 + 7$.

Источник: Подобные задачи встречаются на Codeforces Div.1 D и ICPC Regional.

Решение.

1. Находим фундаментальное решение $(x_1, y_1)$ через разложение $\sqrt{D}$.
2. Используем быстрое возведение в степень в кольце $\mathbb{Z}[\sqrt{D}]$ для вычисления $(x_1 + y_1\sqrt{D})^N \pmod{10^9+7}$.

Умножение в $\mathbb{Z}[\sqrt{D}]$

​]: $(a + b\sqrt{D})(c + d\sqrt{D}) = (ac + Dbd) + (ad + bc)\sqrt{D}$.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128 lll;

const ll MOD = 1e9 + 7;

// Кольцо Z[sqrt(D)] mod MOD
struct ZSqrtD {
    ll a, b; // a + b*sqrt(D)
    ll D;
    ZSqrtD(ll a, ll b, ll D) : a(a%MOD), b(b%MOD), D(D) {}
    ZSqrtD operator*(const ZSqrtD& o) const {
        ll na = ((__int128)a*o.a + (__int128)D*b%MOD*o.b) % MOD;
        ll nb = ((__int128)a*o.b + (__int128)b*o.a) % MOD;
        return {na, nb, D};
    }
};

ZSqrtD power(ZSqrtD base, ll exp) {
    ZSqrtD result{1, 0, base.D};
    for (; exp; exp >>= 1) {
        if (exp & 1) result = result * base;
        base = base * base;
    }
    return result;
}

// Нахождение фундаментального решения Пелля через __int128
pair<lll, lll> pell_fund(ll D) {
    ll a0 = (ll)sqrt((double)D);
    while ((a0+1)*(a0+1) <= D) a0++;
    while (a0*a0 > D) a0--;
    if (a0*a0 == D) return {1, 0};

    ll m = 0, d = 1, a = a0;
    lll pp = 1, pc = a0, qp = 0, qc = 1;

    for (int iter = 0; iter < 4000000; iter++) {
        m = d*a - m;
        d = (D - m*m)/d;
        a = (a0+m)/d;
        lll pn = a*pc+pp, qn = a*qc+qp;
        pp=pc; pc=pn; qp=qc; qc=qn;
        if (pp*pp - (lll)D*qp*qp == 1) return {pp, qp};
    }
    return {-1,-1};
}

int main(){
    ll D, N; cin >> D >> N;

    auto [x1, y1] = pell_fund(D);
    if (x1 == -1) { cout << "No solution\n"; return 0; }

    // x1 и y1 могут быть большими — берём mod MOD
    ll x1m = (ll)(x1 % MOD);
    ll y1m = (ll)(y1 % MOD);

    // (x1 + y1*sqrt(D))^N mod MOD
    ZSqrtD base{x1m, y1m, D % MOD};
    ZSqrtD res = power(base, N);

    cout << "x_" << N << " = " << res.a << " (mod " << MOD << ")\n";
    cout << "y_" << N << " = " << res.b << " (mod " << MOD << ")\n";
    return 0;
}

Важный нюанс. При очень большом $D$ (до $10^{13}$) или большом периоде фундаментального решения числа $x_1, y_1$ могут быть астрономически большими, и для взятия по модулю нужно работать с ними непосредственно в алгоритме поиска (храня $p_k \bmod M$, $q_k \bmod M$ и $p_k$ в __int128 для проверки).

Задачи для самостоятельного решения

Типичные ошибки

1. Ошибка в начальных условиях рекуррентности. Правильно: $p_{-1} = 1$, $p_0 = a_0$, $q_{-1} = 0$, $q_0 = 1$. Частая ошибка: начать с $p_0 = 0$, $p_1 = a_0$.

2. Потеря точности при вычислении $a_0 = \lfloor\sqrt{D}\rfloor$

   ​⌋. Стандартный sqrt в C++ даёт погрешность для больших $D$. Всегда корректируйте:

   ll a0 = (ll)sqrt((double)D);
   while ((a0+1)*(a0+1) <= D) a0++;
   while (a0*a0 > D) a0--;

3. Бесконечный цикл при $D$ — полный квадрат. Перед алгоритмом Пелля проверьте $a_0^2 = D$.

4. Неправильное определение конца периода. Период заканчивается при $a_k = 2a_0$, но это может произойти на первом же шаге (например, $\sqrt{3} = [1; \overline{1, 2}]$

   ​=[1;1,2​] — нет, $\sqrt{2} = [1; \overline{2}]$).

5. Неучёт нечётного периода. Если длина периода нечётна, подходящая дробь $p_{l-1}/q_{l-1}$ решает $x^2 - Dy^2 = -1$, а не $= 1$. Для уравнения $= 1$ нужно взять $p_{2l-1}/q_{2l-1}$.

6. Переполнение в проверке. Числа $p_k, q_k$ растут экспоненциально. При $D \le 10^6$ фундаментальное решение может иметь сотни цифр. Используйте __int128 или BigInteger.

Совет профессионала

Атака Вистра на RSA и её применение в CTF. В задачах CTF-стиля (и иногда в олимпиадах) встречается следующий сценарий: даны $N = pq$ и $e$ (публичный ключ RSA), причём $d < N^{1/4}$ (малый секретный показатель). Тогда из соотношения $ed \equiv 1 \pmod{\phi(N)}$ следует $|e/N - k/d| < 1/(2d^2)$, то есть $k/d$ является подходящей дробью разложения $e/N$ в непрерывную дробь. Перебирая подходящие дроби и проверяя $\phi(N) = (ed - 1)/k$, $p + q = N - \phi(N) + 1$, находим факторизацию.

// Атака Вистра (концептуально)
auto cf = rational_cf(e, N);
auto conv = convergents(cf);
for (auto [k, d] : conv) {
    if (k == 0) continue;
    if ((e * d - 1) % k == 0) {
        ll phi = (e * d - 1) / k;
        // N = pq, phi = (p-1)(q-1) = N - (p+q) + 1
        // p + q = N - phi + 1
        // дискриминант: (p+q)^2 - 4N должен быть полным квадратом
        ll s = N - phi + 1;
        ll disc = s*s - 4*N;
        if (disc > 0) {
            ll sq = (ll)sqrt((double)disc);
            if (sq*sq == disc) {
                cout << "p = " << (s + sq)/2 << ", q = " << (s - sq)/2 << "\n";
            }
        }
    }
}

Итог

Непрерывные дроби — элегантный мост между числами рациональными, иррациональными и целыми диофантовыми задачами:

1. Разложение в НД вычисляется через алгоритм Евклида за $O(\log \min(p,q))$.
2. Подходящие дроби дают наилучшие рациональные приближения и вычисляются за линейное время через рекуррентность.
3. Периодические НД характеризуют квадратичные иррациональности; разложение $\sqrt{D}$ вычисляется за $O(\sqrt{D})$.
4. Уравнение Пелля всегда имеет решения (при $D$ не квадрате), которые выражаются через подходящие дроби $\sqrt{D}$.
5. Все решения уравнения Пелля порождаются фундаментальным через умножение в кольце $\mathbb{Z}[\sqrt{D}]$, а $n$-е решение вычисляется за $O(\log n)$.