Расширенный алгоритм Евклида

Обычный Евклид отвечает на вопрос «чему равен НОД». Расширенный отвечает на более сильный: как выразить НОД через сами числа.

Формально: даны целые $a$ и $b$, не равные нулю одновременно. Нужно найти такие целые $x$ и $y$, что

$$a x + b y = \gcd(a, b)$$

Уравнение такого вида называется диофантовым — это уравнение, решения которого ищут в целых числах.

Первое, что здесь неочевидно: решение существует всегда. Второе — что находится оно тем же алгоритмом Евклида, просто раскрученным в обратную сторону.

Откуда берутся формулы

Евклид спускается по цепочке

$$\gcd(a, b) \to \gcd(b,\; a \bmod b) \to \dots \to \gcd(d, 0)$$

Внизу цепочки задача решается глазами: для пары $(d, 0)$ подходит $x = 1$, $y = 0$, потому что $d \cdot 1 + 0 \cdot 0 = d$.

Осталось научиться подниматься на шаг вверх. Пусть для пары $(b,\; a \bmod b)$ мы уже знаем ответ:

$$b \cdot x_1 + (a \bmod b) \cdot y_1 = d$$

Подставим определение остатка $a \bmod b = a - \left\lfloor a/b \right\rfloor \cdot b$:

$$b x_1 + \left(a - \left\lfloor \tfrac{a}{b} \right\rfloor b\right) y_1 = d$$

Раскроем скобки и соберём отдельно то, что стоит при $a$, и то, что стоит при $b$:

$$a \cdot y_1 + b \cdot \left(x_1 - \left\lfloor \tfrac{a}{b} \right\rfloor y_1\right) = d$$

Слева получилось ровно то, что нужно, — комбинация $a$ и $b$. Значит,

$$x = y_1, \qquad y = x_1 - \left\lfloor \frac{a}{b} \right\rfloor \cdot y_1$$

Это весь алгоритм. Никакой отдельной идеи в нём нет — только аккуратно раскрытые скобки.

flowchart TD
    A["gcd(4, 3)"] --> B["gcd(3, 1)"]
    B --> C["gcd(1, 0)"]
    C -->|"x=1, y=0"| D["база: 1·1 + 0·0 = 1"]
    D -->|"поднимаемся"| E["gcd(3,1): x=0, y=1<br/>3·0 + 1·1 = 1"]
    E -->|"поднимаемся"| F["gcd(4,3): x=1, y=-1<br/>4·1 + 3·(-1) = 1"]

Спуск ничем не отличается от обычного Евклида. Вся работа происходит на подъёме: каждый возврат из рекурсии пересчитывает пару коэффициентов по формулам выше.

Код

long long gcdExt(long long a, long long b, long long &x, long long &y) {
    if (b == 0) { x = 1; y = 0; return a; }
    long long x1, y1;
    long long d = gcdExt(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return d;
}

Здесь x и y передаются по ссылке: функция должна вернуть три числа сразу — сам НОД и два коэффициента, — а return в C++ один. Амперсанд означает, что внутри функции мы меняем не копию, а сами переменные вызывающего.

Отдельные x1, y1 нужны обязательно. Новые значения зависят и от старого $x$, и от старого $y$ одновременно, поэтому перезаписывать их на месте нельзя — второе присваивание увидело бы уже испорченное первое.

И да: в отличие от обычного Евклида, этот алгоритм пишут именно рекурсивно. Итеративная версия существует, но читается заметно хуже, а выигрыш нулевой — глубина рекурсии здесь логарифмическая, до сорока кадров для чисел до $10^{18}$.

Как это разворачивается

Возьмём $\gcd(4, 3)$ и проследим за стеком.

вызов	что вернулось снизу	$x$	$y$	проверка
gcdExt(1, 0)	база	$1$	$0$	$1 \cdot 1 + 0 \cdot 0 = 1$
gcdExt(3, 1)	$x_1 = 1,\ y_1 = 0$	$0$	$1 - 3 \cdot 0 = 1$	$3 \cdot 0 + 1 \cdot 1 = 1$
gcdExt(4, 3)	$x_1 = 0,\ y_1 = 1$	$1$	$0 - 1 \cdot 1 = -1$	$4 \cdot 1 + 3 \cdot (-1) = 1$

Обратите внимание: равенство выполняется на каждом уровне, а не только в конце. Это удобно при отладке — можно поставить проверку внутрь функции и сразу увидеть, на каком шаге сломалось.

Заодно видно, что коэффициенты бывают отрицательными. Это нормально и неизбежно: если $a$ и $b$ положительны, а НОД меньше обоих, то без минуса комбинацию не собрать.

Общее уравнение ax + by = c

На практике справа стоит не НОД, а произвольное $c$. Тогда:

Решение существует тогда и только тогда, когда $c$ делится на $g = \gcd(a, b)$. В одну сторону очевидно: левая часть всегда кратна $g$. В другую — если $c = g \cdot t$, домножим найденную пару на $t$.

// решает a*x + b*y = c; возвращает false, если решений нет
bool diophantine(long long a, long long b, long long c, long long &x, long long &y) {
    long long g = gcdExt(a, b, x, y);
    if (c % g != 0) return false;
    x *= c / g;
    y *= c / g;
    return true;
}

Найденное решение — не единственное. Все остальные получаются сдвигом:

$$x_k = x_0 + k \cdot \frac{b}{g}, \qquad y_k = y_0 - k \cdot \frac{a}{g}, \qquad k \in \mathbb{Z}$$

Проверяется подстановкой: добавка к левой части равна $k \cdot \frac{ab}{g} - k \cdot \frac{ab}{g} = 0$. Делить нужно именно на $g$ — с меньшим шагом получатся не все решения, с бо́льшим часть потеряется.

Отсюда решаются типовые постановки: «найдите решение с наименьшим неотрицательным $x$» — берём $x_0 \bmod \frac{b}{g}$; «сколько решений в заданном отрезке» — считаем допустимые $k$ через два неравенства.

Зачем это нужно на самом деле

Обратный элемент по модулю. Если $\gcd(a, m) = 1$, то из $a x + m y = 1$ следует $a x \equiv 1 \pmod m$, то есть $x$ — обратный к $a$. В отличие от способа через малую теорему Ферма, здесь модуль не обязан быть простым — достаточно взаимной простоты. Это единственный доступный путь, когда модуль вида $2^{32}$ или просто составной.

Задачи про кувшины и гири. «Есть сосуды на 7 и 11 литров, отмерьте 5» — это буквально $7x + 11y = 5$, где знак означает направление переливания.

Линейные сравнения и китайская теорема — целиком стоят на расширенном Евклиде, про них отдельная статья.

Где ломается

Главная ловушка — переполнение. Промежуточные $x$ и $y$ могут по модулю доходить до $\max(a, b)$, а после домножения на $c/g$ — вылететь за long long. Если $a$, $b$ и $c$ порядка $10^{18}$, считать нужно в __int128 или сразу приводить по модулю.

Вторая — отрицательные входные данные. Формулы верны и для них, но % в C++ возвращает остаток со знаком делимого, и «наименьшее неотрицательное решение» придётся нормализовать вручную: ((x % n) + n) % n.