§ Теоретическое обоснование

Сегодня я попытаюсь разобрать алгоритм Евклида в его расширенном варианте. Что это значит? Что значит "расширенный"? Это просто. Вот мы нашли НОД(a,b)=d, но нам надо узнать вот что:
ax + by = d
То есть, нам надо найти такой x и такой y, при умножении на которое у нас получится d (см. формулу выше). Вот пример. Найдем НОД(3,12) и он равен 3, потому что как 12, так и 3 делятся на 3. Хорошо, мы его нашли, выпишем:
3x + 12y = 3
Где a=3, b=12, d=НОД(3,12). И как теперь узнать, какие будут x, y? Автор классной книги, его звать Кнут (и пряник, ага), нашел этот отличный метод и назвал его "Расширенный алгоритм Евклида". Начнем по-порядку. Сначала надо вспомнить то выражение, которое мы разбирали в прошлой статье про алгоритм Евклида:
a = bq + r
Где НОД(a,b) = НОД(b,r). Это было доказано ранее. Теперь посмотрим, как происходит обычно поиск НОД по этому алгоритму. Как мы помним, чтобы перейти к следущей итерации, надо разделить "a" на "b", получить q, и остаток, и применить их как новые "a" и "b".
Шаг Поиск НОД Разложение остатков Описание
1 a = bq_1 + r_1 r_1 = ax_1 + by_1 q_1 = [a / b], \ r_1 = a - bq_1 , теперь они вместо (a, b)
2 b = r_1q_2 + r_2 r_2 = ax_2 + by_2 q_2 = [b / r_1], \ r_2 = b - r_1q_2 , теперь они вместо (b, r_1)
3 r_1 = r_2q_3 + r_3 r_3 = ax_3 + by_3 q_3 = [r_1 / r_2], \ r_3 = r_1 - r_2q_3 , теперь они вместо (r_1, r_2)
4 r_2 = r_3q_4 + r_4 r_4 = ax_4 + by_4 q_4 = [r_2 / r_3], \ r_4 = r_2 - r_3q_4 , теперь они вместо (r_2, r_3)
....
n - 1 r_{n-3} = r_{n-2}q_{n-1} + r_{n-1} r_{n-1} = ax_{n-1} + by_{n-1} q_{n-2} = [r_{n-3} / r_{n-2}], \ r_{n-1} = r_{n-3} - r_{n-2}q_{n-1} , теперь они вместо (r_{n-3}, r_{n-2})
n r_{n-2} = r_{n-1}q_{n} r_{n} = 0 Точка остановки
Все, что происходит тут - это просто магия. Обычная магия, которая происходит в фентези. Ложим в пробирку немного математических ингредиентов и начинаем смешивать. Немного выстоять... и зелье готово! Теперь немного придется включить воображение и приглядеться к формулам. Мы заметили, что:
r_{n} = r_{n-2} - r_{n-1}q_{n}
И еще то, что:
r_n = ax_n + by_n
И поэтому:
r_{n-2} = ax_{n-2} + by_{n-2}
r_{n-1} = ax_{n-1} + by_{n-1}
Давайте теперь подставим в формулу выше формулу ту, что ниже:
r_n = (ax_{n-2} + by_{n-2}) - (ax_{n-1} + by_{n-1})q_{n}
Раскроем скобки и перегруппируем:
r_n = a(x_{n-2} - x_{n-1}q_{n}) + b(y_{n-2} - by_{n-1}q_{n})
Теперь приглядимся внимательно... так вот что мы видим? Мы видим, что коэффициентами перед a и b у нас стоят
x_n и y_n :
x_n = x_{n-2} - x_{n-1}q_{n}
y_n = y_{n-2} - y_{n-1}q_{n}
Да это просто шикарно! То есть, чтобы вычислить следующий x_{n+1} и y_{n+1} , нужно всего лишь знать предыдущие 2! Вот в этом и есть сила расширенного алгоритма Евклида.
Тут встает вопрос. Вот мы хотим вычислить x_{1} и y_{1} , то есть, нам нужны коэффициенты ( x_{0} , y_{0}, x_{-1} , y_{-1} ), и где их достать? А тут все просто. Взглянем на таблицу. У нас r_3 = r_1 - r_2q_3 , то если мы попробуем найти r_2 , то в таком случае r_1=a , и r_2=b , судя по таблице.
Это значит, что:
r_{-1}=ax_{-1}+by_{-1}=a
r_{0}=ax_{0}+by_{0}=b
Отсюда:
x_{-1}=1, y_{-1}=0
x_{0}=0, y_{0}=1
Вот такой вот он, расширенный алгоритм Евклида. И результатом НОД будет r_{n-1} = ax_{n-1} + by_{n-1} .

§ Программная реализация

А теперь время программного обеспечения.
1// Расширенный НОД
2int extnod(int a, int b, int* x, int* y) {
3
4    int q, r;
5    int x2 = 1, y2 = 0, // x(i-2)=1, y(i-2)=0
6        x1 = 0, y1 = 1; // x(i-1)=0, y(i-1)=1
7
8    do {
9
10        // Алгоритм Евклида
11        q = a / b;
12        r = a % b;
13        a = b;
14        b = r;
15
16        // Расширенный алгоритм Евклида
17        *x = x2 - x1*q; // x(i) = x(i-2) - x(i-1)*q
18        *y = y2 - y1*q; // y(i) = y(i-2) - y(i-1)*q
19
20        // "Сдвиг истории" назад
21        x2 = x1; y2 = y1; // Старый x(i-1), y(i-1) перемещается в x(i-2), y(i-2)
22        x1 = *x; y1 = *y; // Новый  x(i), y(i) становится x(i-1), y(i-1)
23
24    } while (r > 0); // Перебирать, пока не станет r = 0
25
26    // И записать разультат, где x = x(i-1), y = y(i-1)
27    // т.к. ранее мы "сдвинули" в (x2, y2) старые (x1, y1)
28    *x = x2; *y = y2;
29
30    // Вернуть значение НОД
31    return a;
32}