Прелюдия

1. Понятие степени с целым положительным показателем может быть определено на любой алгебраической структуре с ассоциативной бинарной операцией. Можно возводить в степень числа, перестановки, матрицы, многочлены… Любой из приведённых ниже алгоритмов пригоден для всех указанных случаев, если соответствующим образом реализовать оператор умножения.
2. Если в структуре есть нейтральный элемент [latex] e [/latex], то можно определить возведение в нулевую степень как [latex] a^{0} = e [/latex]. На практике мы обычно имеем дело именно с этим случаем
3. Для каждого обратимого элемента структуры  [latex] \left( M, * \right) [/latex] можно определить степени с отрицательным показателем, положив  [latex] a^{n} = a^{-(-n)} [/latex]  для всех  [latex] n < 0 [/latex]. Для обработки этого случая достаточно в начале каждого из приведённых ниже алгоритмов реализовать следующий псевдокод:
   
   Python

   if (n < 0): if (is_invertible(a)): return (inverse(a), -n) else: error("Элемент необратим")

   1 2 3 4 5

   if (n < 0): if (is_invertible(a)): return (inverse(a), -n) else: error("Элемент необратим")

   где функция is_invertible проверяет является ли элемент обратимым, а функция inverse возвращает для обратимого элемента его обратный.

Итак, пусть есть непустое множество  [latex] M [/latex], на котором задана ассоциативная бинарная операция  [latex] * [/latex], относительно которой в  [latex] M [/latex] есть нейтральный элемент  [latex] e [/latex].  Будем считать, что у нас перегружен оператор умножения  [latex] * [/latex], а также что есть алгоритм  [latex] e() [/latex], возвращающий нейтральный элемент структуры  [latex] \left(M,*\right) [/latex]. Кроме того, деление у нас целочисленное (как в С/С++ и Java), т.е. результатом вычисления выражения  [latex] 5 / 2 [/latex]  является число  [latex] 2 [/latex].

Самой «дорогой» по времени выполнения операцией, встречающейся в этих алгоритмах является умножение, поэтому под сложностью алгоритма возведения в степень будем понимать функцию  [latex] T(n) [/latex], сопоставляющую показателю степени  [latex] n [/latex]  количество умножений, необходимое рассматриваемому алгоритму для возведения в степень  [latex] n [/latex]  произвольного элемента из  [latex] M [/latex].

NB! Алгоритмы приведены на языке Python. Я видел, что этот подход используется во многих онлайн-курсах по Computer Science и многих вводных курсах в программирование. На выбор Python в качестве языка-иллюстратора есть, по крайней мере, три причины:

• синтаксис минималистичный и, как говорят в английском, self-explanatory, поэтому код Python похож на псевдокод, легко воспринимается и очень легко переводится на любой другой язык программирования
• код получается очень компактным: вложенность конструкций регулируется отступами, поэтому нет большого количество скобок; отсутствие необходимости указывать тип переменных и возвращаемый тип функций существенно сокращает используемое количество ключевых слов
• наличие множественного возврата из функций и параллельного присваивания делает код в большей мере соответствующим интуитивному восприятию алгоритмов

Наивный  алгоритм

Наивный алгоритм основан на том, что  [latex] a^{n} = [/latex] [latex]\begin{cases}e & \text{} n=0 \\a * a^{n-1} & \text{} n>0\end{cases} [/latex].

Доказательство корректности обеих версий алгоритма легко проводится индукцией с учётом отмеченного выше рекуррентного соотношения. В итеративной реализации очевидно, что для завершения алгоритма при  [latex] n > 0 [/latex]  необходимо  [latex] n [/latex]  умножений.

Для рекурсивной реализации это также нетрудно доказать. Отметим, что  [latex] T(0) = 0 [/latex]  и что при  [latex] n > 0 [/latex]

[latex] T(n) = 1 + T(n-1) [/latex]

(нам нужно столько умножений, сколько нужно для вычисления  [latex] a ^ {n-1} [/latex] , и ещё одно умножение для вычисления  [latex] a * pow(a, n-1) [/latex]). А поэтому

[latex] T(1) = 1 + T(n-1) = 1 + 0 = 1 [/latex]

если  [latex] T(n) = n [/latex], то  [latex] T(n+1) = 1 + T(n) = 1 + n [/latex]

и, по принципу математической индукции, [latex] T(n) = n [/latex].

Время линейно зависит от [latex] n [/latex], поэтому для возведения в степень с показателем, например, [latex] 10^{18} [/latex] понадобится… ну очень много времени. А ведь в криптографии приходится иметь дело с числами длиной в сотни и иногда даже тысячи битов, поэтому данный алгоритм для возведения в столь большие степени абсолютно неприменим. В олимпиадных задачах тысячебитная степень вряд ли встретится, но вот степень с восемнадцатью нулями вполне возможна. Поэтому у нас есть все причины перейти к следующему алгоритму.

Быстрый рекурсивный алгоритм

Этот алгоритм основан на том, что  [latex] a^{n} = [/latex][latex]\begin{cases}\left(a^{2}\right)^{\frac{n}{2}} & \text{} n\equiv_{2} 0 \\a * \left(a^{2}\right)^{\frac{n}{2}} & \text{} n \not \equiv_{2} 0 \end{cases} [/latex].

Доказательство корректности

Докажем следующую теорему: для любых  [latex] a \in M [/latex]  и  [latex] n \in \mathbb{N}_{0}[/latex], поданных на вход алгоритма pow, этот алгоритм завершается и возвращает [latex] a ^ {n} [/latex]. Доказательство будем проводить индукцией по  [latex] n [/latex]. При  [latex] n = 0 [/latex]  или  [latex] n = 1 [/latex]  алгоритм корректно завершается при любом  [latex] a \in M [/latex]. Пусть  [latex] n > 1 [/latex]  и для всех  [latex] 0 \leq k < n [/latex]  алгоритм корректно завершается при любом  [latex] a \in M [/latex].

Подадим на вход алгоритма произвольное  [latex] a \in M [/latex]  и рассмотрим два случая, в зависимости от чётности/нечётности числа  [latex] n [/latex]. В обоих случаях алгоритм делает рекурсивный вызов  [latex] pow \left( a \cdot a, n/2 \right) [/latex]. По предположению индукции, на входной паре  [latex] \left(a \cdot a, n/2 \right) [/latex]  алгоритм завершается. Значит, на входной паре  [latex] \left( a, n \right) [/latex]  алгоритм также завершится. Осталось показать, что он возвратит правильный результат.

Если  [latex] n [/latex]  чётное, то алгоритм возвращает

[latex] pow \left( a \cdot a, n/2 \right) = {\left( a \cdot a \right)} ^{\frac{n}{2}} = a^{n} [/latex]

Если  [latex] n [/latex]  нечётное, то алгоритм возвращает

[latex] a \cdot pow \left( a \cdot a, n/2 \right) = a \cdot {\left( a \cdot a \right)} ^{\frac{n-1}{2}} = a^{n} [/latex]

Первая оценка времени

Проанализируем время работы. [latex] T(0) = T(1) = 0 [/latex], а при  [latex] n > 1 [/latex]

[latex] T\left( n \right) = [/latex] [latex]\begin{cases}1 + T(n/2) & \text{} n \equiv_{2} 0 \\2 + T(n/2) & \text{} n \not \equiv_{2} 0\end{cases} [/latex]

Из приведённого соотношения следует, что для любого  [latex] n > 1 [/latex]

[latex] T(n) \leq 2 + T(n/2) [/latex]

Проведём неформальный анализ функции  [latex] T [/latex]:

[latex] T(n) \leq 2 + T(n/2) \leq 4 + T(n/4) \leq 6 + T(n/8) \leq \dots[/latex]

Делить число  [latex] n [/latex]  на  [latex] 2 [/latex]  до получения нуля мы можем примерно столько раз, сколько цифр в двоичной записи числа, а их примерно  [latex] \log_2 n [/latex]. При каждом делении мы увеличиваем верхнюю границу на  [latex] 2 [/latex], поэтому естественно предположить следующую оценку:  [latex] T(n) \leq 2 \log_2 n [/latex].

Проведём формальное доказательство этой оценки индукцией по  [latex] n > 1 [/latex]. При  [latex] n = 2 [/latex]  оценка справедлива, поскольку

[latex] T(2) = 1 + T(1) = 1 + 0 = 1 \leq 2 = 2 \cdot 1 = 2 \cdot log_2 2 [/latex]

Пусть  [latex] n > 2 [/latex]  и для всех  [latex] 1 < k < n [/latex]  доказываемая оценка справедлива. Тогда

[latex] T(n) \leq 2 + T(n/2) \leq 2 + 2 log_2 (n/2) = 2 + 2 log_2 n — 2 log_2 2 = 2 + 2 log_2 n — 2 = 2 log_2 n[/latex]

где второе неравенство обосновано индукционным предположением, применённым к числу  [latex] 1 < n/2 < n [/latex].

Таким образом, для всех  [latex] n > 1 [/latex]  справедливо неравенство  [latex] T(n) \leq 2 \log_2 n [/latex]. А это, по определению нотации «О-большое» означает, что  [latex] T(n) = O(\log n) [/latex]. Мы не указываем основание логарифма, поскольку значения логарифма одного и того же числа по двум разным основаниям отличаются на константный множитель, а константы О-большое «съедает».

Вторая оценка времени

У нас есть асимптотическая оценка сложности алгоритма, но мы также можем посчитать точное количество умножений. Сначала, как и выше, проведём неформальный анализ. При каждом рекурсивном вызове показатель степени делится нацело на  [latex] 2 [/latex]. Значит, можем ожидать, что будет примерно  [latex] \log_2 n [/latex]  рекурсивных вызовов. При каждом вызове будет обязательно произведено умножение при вычислении выражения  [latex] a*a [/latex]. Деление нацело на  [latex] 2 [/latex]   равносильно отсечению крайней правой цифры в двоичной записи числа. Если при текущем вызове значение  [latex] n [/latex]  являлось нечётным (т.е. если при текущем отсечении последняя цифры была единицей), то будет произведено дополнительное умножение.

Таким образом, можем ожидать, что  [latex] T(n) \approx \log_2 n + s(n) [/latex], где  [latex] s(n) [/latex]  —  количество единиц в двоичной записи числа  [latex] n [/latex]. Обозначим дополнительно  [latex] l(n) = 1 + \left\lfloor \log_2 n \right\rfloor [/latex]  —  количество цифр в двоичной записи числа  [latex] n [/latex]. Поиграв немного со значениями функции  [latex] T [/latex]  можно прийти к следующему выводу: при любом  [latex] n > 1 [/latex]    [latex] T(n) = l(n) + s(n) — 2 [/latex].

Докажем это равенство индукцией по  [latex] n > 1 [/latex]. При  [latex] n = 1 [/latex]  имеем  [latex] T(1) = 0 = 1 + 1 — 2 = l(n) + s(n) — 2 [/latex]. Пусть  [latex] n \geq 2 [/latex]  и для всех  [latex] 1 \leq k < n [/latex]  доказываемая оценка справедлива. Двоичная запись числа  [latex] n [/latex]  получается из двоичной записи числа  [latex] n/2 [/latex]  приписыванием нуля (если  [latex] n [/latex]  чётное) или единицы (если  [latex] n [/latex]  нечётное). С учётом этого замечания получаем:

• если  [latex] n [/latex]  чётное, то  [latex] l(n) = l(n/2) + 1 [/latex],  [latex] s(n) = s(n/2) [/latex]  и  [latex] T(n) = 1 + T(n/2) = 1 + l(n/2) + s(n/2) — 2 = l(n) + s(n) — 2 [/latex]
• если  [latex] n [/latex]  нечётное, то  [latex] l(n) = l(n/2) + 1 [/latex],  [latex] s(n) = s(n/2) + 1 [/latex]  и  [latex] T(n) = 2 + T(n/2) = 2 + l(n/2) + s(n/2) — 2 = l(n/2) + s(n/2) = l(n) — 1 + s(n) — 1 = l(n) + s(n) — 2 [/latex]

Бинарный алгоритм

Быстрый рекурсивный алгоритм можно превратить в быстрый итеративный, просматривая двоичные цифры числа  [latex] n [/latex]. Просматривать их можно слева направо, а можно справа налево. Соответственно получаем две версии алгоритма. Иногда их называют бинарным алгоритмом возведения в степень слева направо (или справа налево). Кормен и соавторы используют название «метод многократного возведения в квадрат».

В англоязычных источниках описанные методы встречаются под названиями «exponentiation by squaring», «exponentiation by repeated squaring», «square-and-multiply exponentiation», «binary exponentiation». Обычно рассматривается один из методов, а второй предлагается в качестве упражнения. Но мы рассмотрим оба.

Схема «слева направо»

Пусть количество цифр в числе  [latex] n [/latex]  есть  [latex] l(n) = t [/latex]:  [latex] \left( n_{t-1} \dots n_{0} \right)_{2} [/latex].  Для каждого  [latex] k \in \left[ 0 \dots t-1 \right] [/latex]  обозначим  [latex] m_{k} = \left( n_{t-1} \dots n_{k} \right)_{2} [/latex].

Если  [latex] k = 0 [/latex], то  [latex] m_k = n[/latex]  и поэтому  [latex] a ^{m_k} = a^{n} [/latex].

Если  [latex] 0 < k \leq t-1 [/latex], то

[latex] m_{k-1} = \left( n_{t-1} \dots n_{k} n_{k-1} \right)_{2} = [/latex] [latex] 2\left( n_{t-1} \dots n_{k} \right)_{2} + n_{k-1}= 2 m_{k} + n_{k-1}[/latex]

и поэтому

[latex] a^{m_{k-1}} = \left( a^{m_{k}} \right) ^ {2} a^{n_{k-1}} = [/latex][latex]\begin{cases}\left(a^{m_{k}}\right)^{2} & \text{} n_{k-1} = 0 \\a * \left(a^{m_{k}}\right)^{2} & \text{} n_{k-1} = 1 \end{cases} [/latex]

Получаем следующий алгоритм:

Отмечу, что процесс «накопления» показателя степени в данном алгоритме по сути идентичен процессу вычисления значения многочлена  [latex]\sum_{k=0}^{t}n_{k}x^{k} [/latex]  в точке  [latex] x = 2 [/latex]  по схеме Горнера.

Главным преимуществом данного алгоритм перед алгоритмом по схеме «справа налево» является тот факт, что на стадии домножения (если текущая цифра —  единица), умножение всегда происходит на один и тот же элемент. Иногда это позволяет существенно сэкономить время, затрачиваемое на умножение. Например, в тестах простоты встречается ситуация степени с малым основанием и большим показателем. В этом случае умножение на основание является операцией существенно более быстрой, чем умножение на произвольное число, которое может возникнуть в схеме «справа налево».
Оценка времени. Деление на  [latex] 2 [/latex]  можно считать быстрой операцией. Поэтому, игнорируя детали работы со списками, время на препроцессинг можно оценить как линейно зависящее от  [latex] l(n) [/latex]. Оценим количество возведений в квадрат / умножений. В случаях  [latex] n = 0 [/latex]  и  [latex] n = 1 [/latex]   не производится ни возведений в квадрат, ни умножений.  Если же  [latex] n > 1 [/latex], то цикл  for  в стадии вычислений отработает  [latex] l(n)-1 [/latex]  шаг и на каждом из них будет произведено возведение в квадрат. При этом «на каждой единице», кроме старшей, будет произведено умножение на  [latex] a [/latex]. В итоге получаем, что при  [latex] n > 1 [/latex]  алгоритму требуется  [latex] l(n) — 1 [/latex]  возведений в квадрат и  [latex] s(n) — 1 [/latex]  обычных умножений, где  [latex] l(n) [/latex]  —  количество цифр в двоичной записи числа  [latex] n [/latex], [latex] s(n) [/latex]  —  количество единиц в двоичной записи  [latex] n [/latex].

Схема «справа налево»

Пусть количество цифр в числе  [latex] n [/latex]  есть  [latex] l(n) = t [/latex]:  [latex] \left( n_{t-1} \dots n_{0} \right)_{2} [/latex].  Для каждого  [latex] k \in \left[ 0 \dots t-1 \right] [/latex]  обозначим  [latex] m_{k} = \left( n_{k} \dots n_{0} \right)_{2} [/latex].

Если  [latex] k = t-1 [/latex], то  [latex] m_k = n[/latex]  и поэтому  [latex] a ^{m_k} = a^{n} [/latex].

Если  [latex] 0 \leq k < t-1 [/latex], то

[latex] m_{k+1} = \left( n_{k+1} n_{k} \dots n_{0} \right)_{2} = [/latex]  [latex] 2^{k+1} n_{k+1} + \left( n_{k} \dots n_{0} \right)_{2}= 2^{k+1} n_{k+1} + m_{k}[/latex]

и поэтому

[latex] a^{m_{k+1}} = \left( a^{ 2^{k+1} } \right)^{n_{k+1}} a^{m_{k}} =[/latex]  [latex]\begin{cases}a^{m_{k}} & \text{} n_{k+1} = 0 \\ a^{2^{k+1}} a^{m_{k}} & \text{} n_{k+1} = 1 \end{cases} [/latex]

Мы могли бы получить алгоритм, похожий на тот, что приведён в схеме «слева направо», но цифры числа можно просматривать «динамически», без предвычисления.

 На примере этого алгоритма рассмотрим, как проводятся доказательства корректности методом инварианта.

Докажем следующую теорему: для любых  [latex] a \in M [/latex]  и  [latex] n \in \mathbb{N}_{0}[/latex], поданных на вход алгоритма pow, этот алгоритм завершается и возвращает [latex] a ^ {n} [/latex].

Доказательство будем проводить методом инварианта: покажем, что перед циклом и после каждого шага цикла справедлив следующий инвариант: [latex] r a^{n} = A^{N} [/latex], где  [latex] A [/latex]  и  [latex] N [/latex]  —  первоначальные значения переменных  [latex] a [/latex]  и  [latex] n [/latex].

Предусловие. Перед началом цикла  [latex] r a^{n} = e A^{N} = A ^ {N} [/latex].

Сохранение. Предположим, что инвариант справедлив перед некоторым шагом цикла и пусть [latex] \rho, \nu, \alpha [/latex]  —  текущие значения переменных  [latex] r, n, a[/latex]. Если  [latex] \nu [/latex] чётное, то после завершения текущего шага будем иметь: [latex] r = \rho [/latex], [latex] a = \alpha ^ {2} [/latex], [latex] n = \frac{\nu}{2} [/latex], а если нечётное, то  [latex] r = \rho a [/latex], [latex] a = \alpha ^ {2} [/latex], [latex] n = \frac{\nu — 1}{2} [/latex]. Непосредственно проверяется, что в обоих случаях справедливо равенство [latex] r a ^ {n} = \rho \alpha ^ {\nu}[/latex]. По предположению об инварианте, правая часть равна [latex] A ^ {N} [/latex]. Значит, инвариант после завершения текущего шага сохраняется.

Завершение доказательства.

Если  [latex] N = 0 [/latex], то цикл ни разу не запустится и алгоритм завершится, возвратив  [latex] e = A^{0} [/latex], т.е. завершится корректно.

Если  [latex] N = 1 [/latex], то цикл отработает единственный шаг и алгоритм завершится, возвратив  [latex] A = A^{1} [/latex], т.е. завершится корректно.

Если  [latex] N > 0 [/latex], то, поскольку на каждому шаге мы целочисленно делим положительное число [latex] N [/latex]  на  [latex] 2 [/latex], за конечное число шагов цикл завершится. С помощью индукции нетрудно показать, что точное число шагов в этом случае составит  [latex] l(N) [/latex]  —  количество цифр в двоичной записи числа  [latex] N [/latex]. После завершения цикла переменная  [latex] n [/latex]  содержит значение  [latex] 0 [/latex]  и после последнего шага, по выше доказанному, инвариант сохранялся. Значит, [latex] A ^ { N} = r \cdot a ^ {n} = r \cdot a ^ {0} = r \cdot e = r [/latex], т.е. алгоритм к.

Оценка времени. При  [latex] n > 1 [/latex]  цикл  while, как было сказано выше отработает  [latex] l(n) [/latex]  шагов, на каждом из которых, кроме последнего, производится возведение в квадрат. При этом на каждом шаге с единицей мы производим умножение. Однако первое умножение — это фактически умножение на  [latex] e [/latex], а такие умножения принято исключать из подсчёта. В итоге получаем, что при  [latex] n > 1 [/latex]  алгоритму требуется  [latex] l(n) — 1 [/latex]  возведений в квадрат и  [latex] s(n) — 1 [/latex]  обычных умножений, где  [latex] l(n) [/latex]  —  количество цифр в двоичной записи числа  [latex] n [/latex], [latex] s(n) [/latex]  —  количество единиц в двоичной записи  [latex] n [/latex].

Множественное возведение в степень

Напомню, что в записи  [latex] a ^ {b} [/latex]  элемент  [latex] a [/latex]  называют основанием степени, а число  [latex] b [/latex]  —  показателем степени. Кроме случая, когда основание и показатель фиксированы, встречаются задачи, когда оснований/показателей много. Эти задачи встречаются в литературе под общим названием множественного возведения в степень (multiexponentiation). Основные приложения — криптография и тестирование простоты.

Три основные задачи множественного возведения в степень:

• основание фиксировано, показатели меняются:   [latex] a^{b_1}, \dots , a^{b_s} [/latex]
• основания меняются, показатель фиксированный:  [latex] {a_1}^{b}, \dots , {a_s}^{b} [/latex]

• вычисление произведения степеней  [latex] {a_1}^{b_1} \cdot \dots {a_s}^{b_s} [/latex]

В каждом из трёх случаев умения эффективно вычислять одну степень недостаточно.

Некоторые из алгоритмов, разработанных для решения этих задач, являются модификациями  [latex] 2^k [/latex]-арного алгоритма, описанного выше. Для всех трёх задач активно используют алгоритмы с предвычислением некоторых промежуточных степеней и последующим кобминированием их.

Когда показатель степени фиксированный, один из методов — применение аддитивных цепей. Аддитивная цепь длины  [latex] l [/latex]  —  это последовательность  [latex] c_0, \dots, c_{l} [/latex]  целых чисел, начинающаяся с единицы и такая, что каждый следующий член равен сумме двух более ранних:

[latex] c_0 = 1 [/latex]

[latex] \forall k \in \left[ 1 \dots l \right] [/latex]   [latex] \exists i,j \in \left[ 0 \dots k-1 \right] [/latex]   [latex] c_k = c_i + c_j [/latex]

Если задана аддитивная цепь, завершающаяся числом  [latex] n [/latex],  то по ней легко построить процесс возведения в степень  [latex] n [/latex]. Соответственно имеет смысл искать кратчайшую аддитивную цепь, потому что она минимизирует количество умножений. Доказано, что задача нахождения оптимальной аддитивной цепи является NP-полной. Однако иногда поиск такой цепи окупается за счёт последующей экономии на умножениях. Кроме того, существуют эвристические алгоритмы поиска цепей. Обзор аддитивных цепей имеется у Дональда Кнута во втором томе той самой книги.

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

Handbook of Applied Cryptography (Menezes, van Oorschot, Vanstone, глава 14, раздел 6)

A Survey of Fast Exponentiation Methods (Daniel M. Gordon)

Algorithms for multi-exponentiation (Bodo Moller)

Improved Techniques for Fast Exponentiation (Bodo Moller)

Pippenger’s exponentiation algorithm (Daniel J. Bernstein)

Handbook of Elliptic and Hyperelliptic Curve Cryptography (Cohen et al., глава 9)

Модулярное возведение в степень

В языке Python реализована длинная арифметика: мы можем работать с числами произвольной длины, не опасаясь проблем с переполнением. Однако С/С++ и Java этой особенностью похвастаться не могут. Ниже приведена реализация на С++ алгоритма модулярного возведения в степень. Поскольку в основном алгоритме  power  используется умножение по модулю, а модуль предполагается большим, то возникает угроза переполнения.

Чтобы предотвратить переполнение при умножении, заметим, что умножение по модулю — это то же возведение в степень, только в аддитивной терминологии,. Значит, мы можем наравне с алгоритмом  power  возведения в степень реализовать алгоритм  mul  вычисления кратных. Различий будет ровно два: «умножением» будет сложение, «единицей» будет ноль. Теперь проблема переполнения перенеслась на сложение, дав нам больше пространства для манёвра. Поскольку все фигурирующие в алгоритме  mul  числа не превосходят модуль, заключаем: если модуль хотя бы в два раза меньше максимального значения типа (в данном случае — значения типа long long), то переполнение нигде не возникнет.

Приведённый ниже алгоритм также демонстрирует один из возможных способов обработки случая отрицательных показателей степени. Как известно, критерием обратимости элемента в кольцах вычетов является его взаимная простота с модулем. Для проверки этого условия и нахождения обратного элемента я использовал расширенный алгоритм Эвклида (функция ext_gcd), реализованный на основе алгоритма 2.107 на стр.67 в Handbook of Applied Cryptography (Menezes et al.).

Алгоритм (точнее, его неотрицательная часть) протестирован на следующих задачах e-olymp:

273. Возведение в степень

480. Возведение в степень — 2 (требуется обработка переполнения)

Отмечу, что приведённая реализация хотя и является достаточно эффективной для олимпиадных задач, всё же допускает много оптимизаций. Одна из них — ускорение модулярного умножения, например, с помощью алгоритма Монтгомери.

Применение в олимпиадных задачах

Линейные рекуррентные соотношения

Возведение в степень матриц может использоваться для быстрого решения линейных рекуррентных соотношений. Рассмотрим частный случай соотношения глубины [latex] 2 [/latex]:  [latex] x_{n} = A x_{n-1} + B x_{n-2} [/latex]  при  [latex] n \geq 2 [/latex], значения  [latex] x_{0} [/latex]  и  [latex] x_{1} [/latex]  -  произвольные.

Обозначим [latex] X_{n} = \begin{bmatrix}x_{n+1} \\ x_{n}\end{bmatrix} [/latex], [latex] M = [/latex][latex] \begin{bmatrix} A & B \\ 1 & 0 \end{bmatrix} [/latex]. Тогда

[latex] X_{n+1} = \begin{bmatrix}x_{n+2} \\ x_{n+1}\end{bmatrix} = [/latex][latex] \begin{bmatrix}A x_{n+1} + B x_{n} \\ 1 x_{n+1} + 0 x_{n}\end{bmatrix} [/latex][latex] = M X_{n} [/latex]

[latex] X_{n} = M [/latex][latex] X_{n-1} = M^{2} X_{n-2} = \dots = M^{n} X_{0}[/latex]

Предположим, что нам нужно вычислить  [latex] x_n [/latex]. По выше сказанному,

[latex] \begin{bmatrix}x_{n+1} \\ x_{n}\end{bmatrix} [/latex] [latex] = X_{n} = M^{n} X_{0} = [/latex] [latex] \begin{bmatrix} p & q \\ r & s \end{bmatrix} [/latex] [latex] \begin{bmatrix}x_{1} \\ x_{0}\end{bmatrix} [/latex]  [latex] = \begin{bmatrix}p x_{1} + q x_{0} \\ r x_{1} + s x_{0}\end{bmatrix} [/latex]

т.е. нам достаточно возвести матрицу  [latex] M [/latex]  в  [latex] n [/latex]-ую степень и тогда  [latex] x_{n} = r x_{1} + s x_{0} [/latex]

Эта идея легко обобщается на случай рекуррентных соотношений произвольной глубины. Сложность перемножения двух квадратных матриц размера  [latex] k \times k [/latex]  обычным алгоритмом составляет  [latex] k^3 [/latex]  умножений, поэтому получаем оценку  [latex] k^{3} \log N [/latex]  умножений для нахождения [latex] N [/latex]-го члена линейной рекуррентной последовательности глубины  [latex] k [/latex].

Метод проверен на задаче  e-olymp 5197. Числа Фибоначчи

С кодом моего решения на С++ можно ознакомиться здесь.

Разрешимость квадратного уравнения по простому модулю

Пусть  [latex] p > 2 [/latex]  -  простое число и  [latex] a, b, c \in \mathbb{Z}_{p} [/latex]. Можно показать, что разрешимость в  [latex] \mathbb{Z}_{p} [/latex]  квадратного уравнения  [latex] a x^{2} + bx + c = 0 [/latex]  равносильна тому, что в  [latex] \mathbb{Z}_{p} [/latex]  существует квадратный корень из  [latex] D = b^{2} - 4ac [/latex] . По критерию Эйлера,  это имеет место, если и только если  [latex] D^{p-1} \equiv_p 1 [/latex] .

Получаем следующий алгоритм проверки квадратного уравнения на разрешимость: нужно возвести в степень "модуль минус один" дискриминант уравнения и сравнить результат с единицей.

e-olymp 3234. Квадратное уравнение 1

С кодом решения на С++ можно ознакомиться здесь.

Интересно, что при всей простоте проверки на разрешимость, для  нахождения  решения модулярного квадратного уравнения требуются намного более сложные алгоритмы. Более того, эффективного детерминированного алгоритма для этой задачи пока не известно.

Упрощённое возведение в степень по простому модулю

Если модуль  [latex] p [/latex]  -  простое число, то все ненулевые элементы  [latex] \mathbb{Z}_p [/latex]  обратимы. Согласно малой теореме Ферма, для любого ненулевого  [latex] a \in \mathbb{Z}_p [/latex]  справедливо сравнение  [latex] a^{p-1} \equiv_p 1[/latex]. Но тогда  [latex] a ^ {p-2} [/latex]  -  элемент, обратный к  [latex] a [/latex], поскольку  [latex] a \cdot a^{p-2} \equiv_p a^{p-1} \equiv_p 1 [/latex].

Тест простоты Ферма

Малая теорема Ферма утверждает, что для любого простого числа  [latex] p [/latex]  и любого целого  [latex] a [/latex], не делящегося на  [latex] p [/latex], справедливо сравнение  [latex] a^{p-1} \equiv_p 1[/latex]. Верно и обратное утверждение: если целое  [latex] n [/latex]  таково, что для любого  [latex] a [/latex], взаимно простого с  [latex] n [/latex], справедливо сравнение  [latex] a^{n-1} \equiv_n 1 [/latex], то число  [latex] n [/latex]  простое.

На этом наблюдении основан вероятностный тест простоты, называемый тестом Ферма. Для проверки на простоту числа  [latex] n [/latex]  предлагается выбрать случайное  [latex] a \in \left[ 1 \cdots n-1 \right] [/latex], возвести  [latex] a [/latex]  в степень  [latex] n-1 [/latex]  и сравнить результат с единицей: если результат равен единице, возвратить "простое", иначе - возвратить "составное".

Теорема Ферма гарантирует следующие характеристики алгоритма:

• если на вход подано простое число, то алгоритм верно распознает его как простое
• если алгоритм возвратил "составное", то число действительно является составным

Если же алгоритм возвратил "простое", то входное число может как быть, так и не быть простым. Поэтому основной вопрос состоит в следующем: какова вероятность того, что алгоритм, получив на вход составное число, назовёт его простым?

Пусть  [latex] n [/latex]  -  составное число. Рассмотрим множество

[latex] G = \left\{x \in \left[ 1 \cdots n-1 \right] : x^{n-1}\equiv_n 1 \right\} [/latex]

Элементы  [latex] G [/latex]  -  это "лжесвидетели простоты", приводящие к неверному ответу. Алгоритм назовёт составное  [latex] n [/latex]  простым, если и только если случайно выбранное им [latex] a \in \left[ 1 \cdots n-1 \right] [/latex]  окажется элементом множества  [latex] G [/latex]. Вероятность этого события есть  [latex] p = \frac{\left| G \right|}{n-1} [/latex].

Оценим количество элементов во множестве  [latex] G [/latex]. Каждый элемент  [latex] G [/latex]  обратим по модулю  [latex] n [/latex], поскольку  [latex] x x^{n-1} \equiv_n x^{n} \equiv_n 1 [/latex]. Значит, [latex] G [/latex]  является подмножеством  [latex] U(\mathbb{Z}_n) [/latex]  -  группы обратимых по модулю  [latex] n [/latex]  вычетов. При этом если  [latex] x, y \in G [/latex], то

[latex] \left( x \cdot y \right) ^ {n-1} \equiv_n \cdot x^{n-1} y ^{n-1} \equiv 1 \cdot 1 \equiv_n 1[/latex]

т.е.  [latex] xy \in G [/latex]. Как известно, замкнутое подмножество конечной группы является подгруппой, поэтому  [latex] G [/latex]  -  подгруппа  [latex] U(\mathbb{Z}_n) [/latex].

По теореме Лагранжа, порядок подгруппы  [latex] G [/latex]  является делителем порядка группы  [latex] U(\mathbb{Z}_n) [/latex]. Поэтому если  [latex] G \neq U(\mathbb{Z}_n) [/latex], то [latex] \left| G \right| \leq \frac{\left| U(\mathbb{Z}_n) \right|}{2} [/latex]  и в этом случае  [latex] p \leq \frac{1}{2} [/latex]. Если при сделанном относительно множества  [latex] G [/latex]  предположении, применить к числу  [latex] n [/latex]    [latex] k [/latex]  независимых тестов Ферма, то вероятность ошибки составит  [latex] \frac{1}{2^k} [/latex].
Без предположения  [latex] G \neq U(\mathbb{Z}_n) [/latex]  приведённые рассуждения неверны и, к несчастью для теста Ферма, существуют числа (их называют числами Кармайкла), для которых  [latex] G = U(\mathbb{Z}_n) [/latex] и поэтому описанная оценка вероятности для таких чисел не работает. Известно, что чисел Кармайкла бесконечно много, но встречаются они очень редко: статья в Википедии утверждает, что имеется  [latex] 20 138 200 [/latex]  чисел Кармайкла между  [latex] 1 [/latex]  и  [latex] {10}^{21} [/latex].

Кроме того, известно, что каждое число Кармайкла имеет не менее трёх различных простых делителей. Это наблюдение доставляет достаточно эффективный алгоритм проверки на простоту чисел, вмещающихся в стандартные целочисленные типы языков  C/C++  и  Java. Пусть  [latex] n \leq {10}^{21} [/latex]. Если  [latex] n [/latex]  является числом Кармайкла, то у него есть простой делитель, не превосходящий  [latex] {10} ^ {7} [/latex]. Все простые числа до десяти миллионов легко получить, например, решетом Эратосфена. Если число  [latex] n [/latex]  не делится ни на одно из них, то оно не может быть числом Кармайкла. Значит, справедлива полученная выше оценка вероятности и применив к числу  [latex] n [/latex]   три раза тест Ферма, мы

• либо обнаружим, что  [latex] n [/latex] составное
• либо можем предполагать, что оно простое; вероятность ошибки при этом составит не более  [latex] \frac{1}{2^3} = 0.125 [/latex]

Описанный метод проверен на задаче  e-olymp 4316. Простые сложности  (нужно проверить на простоту 5000 чисел, не превосходящих  [latex] {10}^{18} [/latex])

С кодом решения на С++ можно ознакомиться здесь.

На самом деле, для проверки на простоту используют тест Миллера-Рабина, либо ещё более изощрённые алгоритмы. Однако прародителем и идейным вдохновителем многих из этих алгоритмов был именно тест Ферма, поэтому полезно с ним познакомиться и немного поиграться.

Протокол Диффи-Хеллмана

Одним из известнейших приложений алгоритмов возведения в степень является протокол Диффи-Хеллмана. Предположим, что Алиса и Боб хотят иметь общий секрет. Допустим, им нужен общий параметр для некоторого алгоритма шифрования, причём этот параметр должен быть одинаковым для обоих и обменяться им они планируют по открытому каналу связи. Алису и Боба устроит, если параметр мы будем называть «ключом» и он будет элементом некоторой конечной циклической группы. Например, группы обратимых вычетов  [latex] U(\mathbb{Z}_p) [/latex]  или группы точек на некоторой эллиптической кривой, или любой другой.

Итак, пусть  [latex] G [/latex]  —  конечная циклическая группа порядка  [latex] n [/latex]  и  [latex] g \in G [/latex]  —  какой-либо её порождающий элемент. Алиса и Боб всем рассказали про  [latex] g [/latex]. О нём знают они двое и знают злоумышленники, прослушивающие канал связи. Алиса также выбрала случайное число  [latex] a \in \left[1 \dots n \right] [/latex], а Боб выбрал число  [latex] b \in \left[ 1 \dots n \right] [/latex]. Алиса вычислила и отправила Бобу элемент  [latex] g^{a} [/latex], а он ей  —  элемент  [latex] g^{b} [/latex].

Вот тут и начинается магия. В качестве общего секретного параметра Алиса и Боб могут взять элемент  [latex] g^{ab} [/latex]. Хотя Алиса не знает число  [latex] b [/latex], но ей известен элемент  [latex] g^b [/latex]  и она знает, что  [latex] \left(g^{b} \right) ^ {a} = g^{ab} [/latex]. Аналогично Боб может вычислить  [latex] \left(g^{a} \right) ^ {b} = g^{ab} [/latex].

Однако возникает вопрос: что помешает узнать секретный(!) параметр  [latex] g^{ab} [/latex]  злоумышленникам, которые прослушивали канал связи? Оказывается, что в некоторых группах задача нахождения элемента  [latex] g^{ab} [/latex]  по известным  [latex] g^{a} [/latex]  и  [latex] g^{b} [/latex]  при неизвестных  [latex] a [/latex]  и  [latex] b [/latex]  является вычислительно очень сложной и на сегодняшний день не может быть решена за адекватное время. В частности, в качестве  [latex] G [/latex]  можно использовать уже упомянутую группу  [latex] U(\mathbb{Z}_p) [/latex]