Разбор Proggy-Buggy Contest

5 декабря компанией DataArt проводилась международная  юмористическая олимпиада по программированию Proggy-Buggy Contest.

Задачи на ней не были сложными, но решать их нужно было очень быстро: на 13 задач отводилось 42 минуты времени.

В Одесском офисе DataArt было рекордное количество команд-участинков среди всех городов, принимавших мероприятие. Среди участников была команда ONU_Rabid_Pandas (Марченко, Илларионова, Вустянюк), которая подготовила разбор задач соревнования:

Условия

Разбор

Пожалуйста, сообщайте в комментариях обо всех найденных очепятках и непонятках. Мы исправим и ответим на вопросы.

Возведение в целую степень

В заметке описаны четыре простейших алгоритма возведения в целую степень: наивный алгоритм, быстрый рекурсивный алгоритм и две схемы бинарного итеративного алгоритма. На этих примерах проиллюстрированы основные методы доказательства корректности алгоритма: по индукции и методом инварианта.

Приведена реализация модулярного возведения в степень на  С++  с предотвращением переполнения.

Дан краткий обзор материалов по множественному возведению в степень, описаны некоторые приложения темы в олимпиадных задачах.

Описан тест простоты Ферма и основанный на нём эффективный тест простоты для чисел, не превосходящих  [latex] {10}^{21} [/latex].

Прелюдия

Прелюдия

  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]. Для обработки этого случая достаточно в начале каждого из приведённых ниже алгоритмов реализовать следующий псевдокод:
    где функция 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].

Рекурсивно PythonC++Итеративно PythonC++

Доказательство корректности обеих версий алгоритма легко проводится индукцией с учётом отмеченного выше рекуррентного соотношения. В итеративной реализации очевидно, что для завершения алгоритма при  [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].

PythonC++

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

Докажем следующую теорему: для любых  [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. Возведение в степень

273

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

480

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

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

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

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

Возведение в степень матриц может использоваться для быстрого решения линейных рекуррентных соотношений. Рассмотрим частный случай соотношения глубины [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. Числа Фибоначчи

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

3234

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

С++

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

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

Если модуль  [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])

4316

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

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

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

Одним из известнейших приложений алгоритмов возведения в степень является протокол Диффи-Хеллмана. Предположим, что Алиса и Боб хотят иметь общий секрет. Допустим, им нужен общий параметр для некоторого алгоритма шифрования, причём этот параметр должен быть одинаковым для обоих и обменяться им они планируют по открытому каналу связи. Алису и Боба устроит, если параметр мы будем называть «ключом» и он будет элементом некоторой конечной циклической группы. Например, группы обратимых вычетов  [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]

А717б

Условие

Две правые треугольные матрицы [latex] A [/latex] и [latex] B [/latex] порядка [latex] n [/latex] заданы в виде последовательности [latex] \frac{\left( n + 1 \right) n}{2}[/latex] чисел: сначала идёт [latex] n [/latex] элементов первой строки, затем [latex] n-1 [/latex] элемент второй строки, и т.д. (из последней, [latex] n [/latex] -ой строки берётся только [latex] n [/latex] -ый элемент). Нужно получить в аналогичном виде матрицу

б) [latex] A \left( I + B^{2} \right) [/latex], где [latex] I [/latex] — единичная матрица порядка [latex] n [/latex]

Тест

[latex] \begin{bmatrix} 2 & 2 & 2 \\ 0 & 0 & 3 \\ 0 & 0 & 1 \end{bmatrix} \cdot \begin{bmatrix} 1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{bmatrix} = \begin{bmatrix} 4 & 54 & 236 \\
0 & 0 & 111 \\
0 & 0 & 37 \end{bmatrix} [/latex]

Решение

Создадим класс для работы с треугольными матрицами указанного вида и снабдим его основными методами, необходимыми для решения задачи.

У нас будут конструктор матрицы из стандартного потока ввода, конструктор скалярной матрицы, а также будут перегружены операторы сложения, умножения и присваивания.

Код на С++

Ideone (C++)

 

Код на Java

Ideone (Java)

e-olimp 3837. Выражения

Задача

Арифметические выражения, как правило, записываются с операторами между двумя операндами (такая запись называется инфиксной нотацией). Например, (x + y) * (zw) — арифметическое выражение в инфиксной нотации. Однако легче написать программу, вычисляющую значение выражения, когда оно находится в постфиксной нотации (известная как обратная польская нотация). В постфиксной нотации оператор записывается за своими двумя операндами, которые и сами в могут быть выражениями. Например, x y + z w — * является постфиксной нотацией приведенного выше выражения. Заметим, что для такой записи скобки не нужны.

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

  1. push: число кладется на верх стека.
  2. pop: число снимается с вершины стека.

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

a := pop();

b := pop();

push(b O a);

Результатом выражения будет единственное число, оставшееся в стеке.

Предположим, что мы используем вместо стека очередь. Очередь также имеет операции push иpop, но их смысл немного другой:

  1.    push: число вставляется в конец очереди.
  2.    pop: из начала очереди извлекается число.

Можете ли Вы переписать заданное выражение так, чтобы алгоритм, использующий очередь при его вычислении, давал тот же результат, что и алгоритм вычисления со стеком?

 

Пояснения к решению

Если задана правильная строка в постфиксной нотации, то алгоритм её вычисления на основе стека выглядит следующим образом:

  1. Просматриваем строку слева направо
  2. Если встречаем переменную, помещаем её в конец стека/очереди
  3. Если встречаем символ оператора  * , то выполняем следующие действия: снимаем со стека a, снимаем со стека b  и  кладём на стек  b * a.

По условию задачи, мы должны обрабатывать строки, используя вместо стека очередь, поэтому описанные выше действия примут следующий вид:  снимаем с очереди a, снимаем с очереди b  и  кладём в конец очереди  b * a.

Постфиксную запись, рассчитанную на обработку посредством стека, будем называть s-постфиксной, а рассчитанную на обработку очередью — q-постфиксной.

Две записи «вычислятся одинаково» в том и только в том случае, когда соответствующие им деревья синтаксического разбора (англ.: parse tree, далее — ДСР) одинаковы. Идея приведённого ниже алгоритма состоит примерно в следующем: на основе s-постфиксной записи строим ДСР, обходя которое специальным образом, получаем q-постфиксную запись. Нужно обходить дерево по уровням справа налево (под уровнем дерева подразумевается множество всех вершин дерева, равноудалённых от корня; уровни дерева естественным образом нумеруются: корень расположен на уровне 0 и так далее). В приведённом ниже коде функция  get_levels  генерирует уровни ДСР на основе s-постфиксной нотации в виде вектора строк, k-ая компонента которого соответствует k-му уровню ДСР заданной строки, вершины которого (если таковые имеются) перечислены справа налево. Поскольку у ДСР может быть не больше уровней, чем символов в строке, то создаём вектор как раз с таким числом компонент и инициализируем его компоненты пустыми строками.

Корректность функции get_levels можно доказать индукцией по длине строки, отметим, что если строка содержит более одного элемента, то она имеет вид ФGО, где Ф и G — правильные строки, О — символ оператора. Запустить рекурсивно функцию с предпоследней позиции — всё равно, что применить её к строке G. ДСР этой строки является левым поддеревом для ДСР исходной строки и k-ый уровень этого дерева является (k+1)-ым уровнем ДСР исходной строки (именно поэтому вызываем функцию с параметром depth на 1 больше). По индукционному предположению, дойдя до начала строки G, функция корректно «расставит» уровни ДСР этой строки (с учётом + 1) и завершит свою работу.  Аналогично будут правильно расставлены уровни ДСР строки Ф, которое является правым поддеревом исходной строки. Таким образом, все вершины ДСР исходной строки будут правильно расставлены, а поскольку рекурсия вызывается сначала для левого, а потом для правого её поддерева, то и перечислены вершины на каждом уровне будут справа налево.

 

Код на С++

Ссылка на программу в ideone.

 

Код на Java

Ссылка на программу в ideone.

 

e-olimp 122. Горные маршруты

Задача. Горный туристический комплекс состоит из n турбаз, соединенных между собой k горными переходами (другие маршруты в горах опасны). Любой переход между двумя базами занимает 1 день. Туристическая группа находится на базе a и собирается попасть на базу b не более чем за d дней. Сколько существует разных таких маршрутов (без циклов) между a и b?
Ссылка на задачу

Пояснение к решению

Поиск путей представляет из себя модифицированный обход в глубину. Вместо массива «глобально посещённых» вершин, в которые уже нельзя возвращаться при обходе, я использую массив «локально посещённых» вершин. Это позволяет использовать одну и ту же вершину в разных, в том числе очень длинных и извращённых путях (в стандартном варианте ищется кратчайший путь, а нам нужны все: в том числе длинные и извращённые).

Кроме того, добавлен параметр depth — длина рассматриваемого в данный момент пути (или глубина текущего рекуррентного вызова, если угодно). Этот параметр изначально имеет значение 0 и при каждом вызове сравнивается с максимально допустимым значением. Перед рекуррентным вызовом мы увеличиваем значение depth, ведь длина нового пути будет на 1 больше. После вызова, т.е. возвращаясь на предыдущий уровень рекурсии, мы уменьшаем длину пути (ведь мы «отбросили» последнюю вершину, т.е. длина пути уменьшилась на 1).

Алгоритм работает примерно так: сначала посещаем начальную вершину и помечаем её как локально посещённую. Если она — искомая, увеличиваем количество найденных путей на 1 и алгоритм завершает свою работу (путей без циклов, содержащих более одной вершины, в этом случае быть не может). Если же начальная вершина не совпадает с целевой, то для всех вершин, в которые можно попасть из неё, применяем рекуррентно тот же алгоритм. Если таких вершин нет, алгоритм завершается, возвращая 0 (нет ни одного пути из начальной вершины в конечную). Если же такие вершины имеются, то увеличиваем параметр «глубина» и вызываем для них рекуррентно наш алгоритм. После выхода из рекурсии, возвращаем параметр «глубина» к последнему его значению (тем, что было перед вызовом) и, чтобы мы могли использовать данную вершину в других путях помечаем её как «локально свободную».

Стоит отметить, что в общем случае задача поиска всех простых путей (пусть даже удовлетворяющих определённым условиям, как в данной задаче) является NP-полной, т.е. решается исключительно переборными алгоритами, работающими за неполиномиальное время (более эффективных алгоритмов для решения таких задач на сегодняшний день неизвестно; кто их найдёт, получит миллион долларов; я не нашёл). Также отметим, что простая модификация приведённого алгоритма позволяет не только вычислять количество простых путей указанной длины, но и перечислять их.

Решение на С++

Решение на Java

 

АА18

Задача

В строке найти первое слово, которое встречается два раза подряд. Слова разделяются одним или несколькими пробелами. Напечатать его.

Пояснение к решению

Идея простая:

  1. последовательно разбиваем строку на отдельные слова с помощью функций стандартной библиотеки
  2. смотрим, не встретилось ли некоторое слово два раза подряд

Решение

Основной вариант:

Перевод на Java:

 

Старый неаккуратный и медленный вариант: