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

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

$$gcd(a,b)=\underset{k:k|a\wedge k|b}{max}k$$

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

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

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

$$gcd(a,b)=\{\begin{array}{ll}a,& b=0\\ gcd(b,a-b),& b>0\end{array}$$

Здесь предполагается, что $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{array}{ll}a,& b=0\\ gcd(b,amodb),& b>0\end{array}$$

Можно показать, что каждые две итерации меньшее число уменьшится хотя бы в два раза, а следовательно алгоритм работает за $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^{\prime }$ и $y^{\prime }$, когда рекурсивно считали $gcd(b,amodb)$. Иными словами, у нас есть решение $(x^{\prime },y^{\prime })$ для пары $(b,amodb)$:

$$b\cdot x^{\prime }+(amodb)\cdot y^{\prime }=g$$

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

$$b\cdot x^{\prime }+(a-\lfloor \frac{a}{b}\rfloor \cdot b)\cdot y^{\prime }=g$$

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

$$a\cdot \underset{x}{\underbrace{y^{\prime }}}+b\cdot \underset{y}{\underbrace{(x^{\prime }-\lfloor \frac{a}{b}\rfloor \cdot y^{\prime })}}=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$ искомые коэффициенты.