Алгоритм Евклида

Наибольшим общим делителем (англ. greatest common divisor) целых неотрицательных чисел \(a\) и \(b\) называется наибольшее число \(x\), которое делит одновременно и \(a\), и \(b\).

\[ \gcd(a, b) = \max_{k: \; k|a \, \land \, k | b} k \]

Когда оба числа равны нулю, результат не определён — подойдёт сколько угодно большое число. Однако в этом случае мы положим в этом случае \(\gcd\) равным тоже нулю, чтобы можно было использовать следующее правило: если одно из чисел равно нулю, то их \(\gcd\) равен второму числу.

Алгоритм Евклида находит \(\gcd\) двух чисел \(a\) и \(b\) за \(O(\log \min(a, b))\). Он известен ещё с IV века до нашей эры, а возможно и ранее.

Алгоритм основывается на следующей несложной формуле:

\[ \gcd(a, b) = \begin{cases} a, & b = 0 \\ \gcd(b,\, a - b), & b > 0 \end{cases} \]

Здесь предполагается, что \(a > b\).

Докажем корректность этой формулы:

• Если \(g = \gcd(a, b)\) делит и \(a\), и \(b\), то их разность \((a-b)\) тоже будет делиться на \(g\).

• Никакой больший делитель \(d\) числа \(b\) не может делить число \((a-b)\): если \(d > g\), то \(d\) не может делить \(a\), а значит и не делит \((a - b)\).

Прямая рекурсивная реализация:

int gcd(int a, int b) {
    if (a < b)
        swap(a, b);
    if (b == 0)
        return a;
    else
        return gcd(b, a - b);
}

Этот алгоритм может работать долго — например, на паре \((10^9, 1)\) он сделает миллиард итераций. Идея дальнейшей оптимизации в том, чтобы вычитать из \(a\) не одно \(b\) за раз, а столько, чтобы в следующий раз \(a\) и \(b\) уже поменялись местами — чтобы новое \(b\) стало меньше нового \(a\).

Простой способ этого достичь — просто вычесть \(b\) из \(a\) сразу максимально возможное число раз, то есть взять вместо нового \(b\) остаток от деления \(a\) на \(b\):

\[ \gcd(a, b) = \begin{cases} a, & b = 0 \\ \gcd(b,\, a \bmod b), & b > 0 \end{cases} \]

Можно показать, что каждые две итерации меньшее число уменьшится хотя бы в два раза, а следовательно алгоритм работает за \(O(\log \min (a, b))\).

Реализация:

int gcd(int a, int b) {
    if (b == 0)
        return a;
    else
        return gcd(b, a % b);
}

Чуть более быстрая итеративная форма:

int gcd(int a, int b) {
    wihle (b > 0) {
        a %= b;
        swap(a, b);
    }
    return a;
}

В компиляторе gcc для этого уже есть встроенная функция __gcd, которая, впрочем, может непредсказуемо себя вести на отрицательных числах и \((0, 0)\).

Время работы

Оценка \(O(\log n)\) относится к худшему случаю. Асимптотика в среднем это более интересный вопрос, играющий существенную роль в некоторых алгоритмах.

Примечательно, что худшие входные данные для алгоритма — это соседние числа Фибоначчи. На графике они видны как синие точки в пропорциях золотого сечения.

Также иногда полезно знать, что нахождение \(\gcd\) группы из \(n\) чисел от \(1\) до \(A\) будет работать не за \(O(n \log A)\), а за \(O(n + \log A)\) — это несложно доказать по индукции.

Решение диофантовых уравнений

Для обычного использования \(gcd\) не нужно даже знать, как алгоритм устроен — он есть даже в компиляторе.

Расширенный алгоритм Евклида находит, помимо \(g = \gcd(a, b)\), такие целые коэффициенты \(x\) и \(y\), что

\[ a \cdot x + b \cdot y = g \]

Эта модификация алгоритма интересна, потому что с помощью неё можно искать обратный элемент по модулю: такой элемент \(a^{-1}\), что \(a \cdot a^{1} \equiv 1\), что то же самое, что найти решение в целых числах:

\[ a^{-1} \cdot a + k \cdot m = 1 \]

Заметим также, что решений бесконечно много: можно \(x\) увеличить на \(b \cdot y\), а \(y\) уменьшить на \(a \cdot x\), и равенство при этом не изменится.

Алгоритм будет тоже рекурсивный. Пусть мы посчитали нужные коэффициенты \(x'\) и \(y'\), когда рекурсивно считали \(\gcd(b, a \bmod b)\). Иными словами, у нас есть решение \((x', y')\) для пары \((b, a \bmod b)\):

\[ b \cdot x' + (a \bmod b) \cdot y' = g \]

Чтобы получить решение для исходной пары, запишем выражение \((a \bmod b)\) как \((a - \lfloor \frac{a}{b} \rfloor \cdot b)\) и подставим в приведенное выше равенство:

\[ b \cdot x' + (a - \Big \lfloor \frac{a}{b} \Big \rfloor \cdot b) \cdot y' = g \]

Теперь выполним перегруппировку слагаемых (сгруппируем по исходным \(a\) и \(b\)) и получим:

\[ a \cdot \underbrace{y'}_x + b \cdot \underbrace{(x' - \Big \lfloor \frac{a}{b} \Big \rfloor \cdot y')}_y = g \]

Сравнивая это с исходным выражением, получаем, что для иходных \(x\) и \(y\) подходят коэффициенты при \(a\) и \(b\).

Реализация:

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

Эта рекурсивная функция по прежнему возвращает значение \(\gcd(a, b)\), но помимо этого записывает в переданные по ссылке переменные \(x\) и \(y\) искомые коэффициенты.