Теория чисел и кодирование с открытым ключом

Содержание

• Кольцо вычетов по модулю m
• Малая теорема Ферма, китайская теорема об остатках, теорема Эйлера
• Основные алгоритмы теории чисел:
  • алгоритм Евклида,
  • расширенный алгоритм Евклида,
  • быстрое возведение в степень,
  • алгоритм в китайской теореме об остатках,
  • извлечение квадратного корня в поле вычетов по модулю p
• Тест простоты Ферма и кармайкловы числа
• Вероятностный тест простоты Миллера-Рабина
• Кодирование с открытым ключом, схема RSA
• Алгоритмы факторизации (разложения числа на множители):
  • ρ-алгоритм Полларда,
  • p-1 алгоритм Полларда.

В криптографии две идеи привели к настоящей революции в этой области, случившейся вскоре после окончания второй мировой войны. Первая идея — это блочное кодирование. Любое сообщение, любая информация представима в виде последовательности нулей и единиц. Эта последовательность разбивается на блоки фиксированной длины (например, 128 бит), каждый блок шифруется и передается по отдельности. И вторая идея — это сведение шифрования к работе с целыми числами (вернее, элементами колец вычетов). Каждый блок представляет собой двоичную запись некоторого целого числа в ограниченном диапазоне. Шифрование — это инъективное отображение, которое исходное целое число отображает в другое целое число. Поскольку целые числа ограничены, их удобно представлять как элементы колец вычетов $\mathbb{Z}_m$.

Напомним, что кольцо $\mathbb{Z}_m$ — это фактор-кольцо кольца целых чисел по идеалу, порожденному числом $m$. Его элементами являются классы эквивалентности целых чисел по отношению эквивалентности $$ x \sim y \iff m | (x - y) $$ В каждом классе эквивалентности можно выбрать по одному представителю. Такое множество представителей называется системой остатков. Чаще всего используют либо неотрицательную систему остатков, либо симметричную, последняя включает отрицательные и положительные остатки, по абсолютной величине не превосходящие $m/2$. Например, кольцо $\mathbb{Z}_5$ представимо двумя системами остатков: $$ \mathbb{Z}_5 = \{ 0, 1, 2, 3, 4 \} \\ \mathbb{Z}_5 = \{ -2, -1, 0, 1, 2 \} $$ Чтобы выполнить действие с остатками (например, сложение или умножение), надо выполнить это действие с ними как с целыми числами и определить, в какой класс эквивалентности попадает результат. Остаток, представляющий этот класс, и будет результатом операции в $\mathbb{Z}_m$.

(Заметим, что систем остатков бесконечно много, и не слишком грамотно отождествлять кольцо $\mathbb{Z}_m$ с какой-либо конкретной системой остатков. Но так часто делают в книгах, адресованных тем, кто не являются математиками и не слишком уверенно владеют понятием фактор-множества.)

Криптографические алгоритмы работают, как правило, с элементами колец вычетов, а не с целыми числами. Целые числа очень быстро растут при умножении и особенно при возведении в степень, тогда как элементы кольца вычетов всегда представимы остатками в ограниченном диапазоне. Отсюда вытекает их огромное практическое значение.

Самая замечательная, хотя и совсем простая теорема в теории чисел — это Малая теорема Ферма.

Теорема 1 (Малая теорема Ферма). Пусть $p\in \mathbb{Z}$ — простое число. Пусть $b\in \mathbb{Z}$ — произвольное целое число, которое не делится на $p$. Тогда $$ b^{p-1} \equiv 1 \pmod{p}. $$

Можно отказаться от требования, чтобы число $b$ было бы взаимно просто с $p$. Справедлива следующая формулировка Малой теоремы Ферма.

Теорема 1'. Пусть $p\in \mathbb{Z}$ — простое число. Пусть $b\in \mathbb{Z}$ — произвольное целое число. Тогда для любого $k\in \mathbb{Z}$, $k\ge 0$ справедливо сравнение $$ b^{(p-1)k + 1} \equiv b \pmod{p}. $$

В частности, $b^p \equiv b \pmod{p}$.

Теорема 2 (Китайская теорема об остатках). Пусть $m_1, m_2,\dots, m_k\in \mathbb{Z}$ — попарно взаимно простые целые числа. Пусть $m = m_1m_2\dots m_k$. Тогда имеет место изоморфизм колец вычетов: $$ \mathbb{Z}_m \cong \mathbb{Z}_{m_1} \oplus \mathbb{Z}_{m_2} \oplus \dots \oplus \mathbb{Z}_{m_k}. $$

Значение Китайской теоремы об остатках состоит в том, что она сводит изучение колца $\mathbb{Z}_m$ для произвольного $m$ к изучению колец вычетов по модулю степеней простых чисел, т.е. колец вида $\mathbb{Z}_{p^s}$, где $p$ — простое число. Поскольку любое целое число $m$ раскладывается в произведение степеней простых чисел $$ m = p_1^{e_1} p_2^{e_2} \dots p_k^{e_k}, $$ то кольцо $\mathbb{Z}_m$ изоморфно прямой сумме примарных колец: $$ \mathbb{Z}_m \cong \mathbb{Z}_{p_1^{e_1}} \oplus \mathbb{Z}_{p_2^{e_2}} \oplus \dots \oplus \mathbb{Z}_{p_k^{e_k}}. $$

Элементами прямой суммы являются строки длины $k$, где на $i$-м месте стоит элемент из $\mathbb{Z}_{m_i}$. Китайская теорема об остатках утверждает, что все действия с элементами кольза $\mathbb{Z}_m$ реализуются как действия со строками. Элементу $x\in \mathbb{Z}_m$ соответствует строка $$ x \to (x_1, x_2, \dots, x_k),\quad x_i \equiv x \pmod{m_i}. $$ Операции со строками выполняются покомпонентно.

Важно, что изоморфизм в Китайской теореме об остатках конструктивный. В сторону слева направо компоненты $x_i$ вычисляются просто как остатки от деления $x$ на $m_i$. Но и обратное отображение также вычисляется с помощью алгоритма, который будет описан в следующем разделе. Пусть дан набор остатков $r_1, r_2,\dots, r_k$. Тогда этот алгоритм вычисляет число $x\in\mathbb{Z}$ такое, что $x \equiv x_i \pmod{m_i}$ для $i=1, 2,\dots, k$.

Как пример применения Китайской теоремы об остатках, рассмотрим следующую задачу. Сколько квадратных корней из единицы в кольце вычетов $\mathbb{Z}_{105}$? (Как вариант — найти все корни из единицы.) Ответ: 8, потому что $105 = 3\cdot 5\cdot 7$ и по Китайской теореме об остатках $$ \mathbb{Z}_{105} \cong \mathbb{Z}_{3} \oplus \mathbb{Z}_{5} \oplus \mathbb{Z}_{7}. $$ Кольцо $\mathbb{Z}_p$ является полем тогда и только тогда, когда $p$ простое. Но в поле по теореме Безу может быть не больше двух квадратных корней из единицы (т.к. они являются корнями уравнения $x^2 - 1 = 0$, а число корней многочлена по теореме Безу не превосходит его степени). Поэтому все корни из единицы $r_j$ в кольце $\mathbb{Z}_{3} \oplus\mathbb{Z}_{5} \oplus \mathbb{Z}_{7}$ имеют вид $$ r_j = (\pm 1, \pm 1, \pm 1), $$ и всего их $2^3=8$. Алгоритм в Китайской теореме об остатках позволяет вычислить их все: $$ \begin{array}{rcl} 1 & \cong & (1, 1, 1) \\ 76 & \cong & (1, 1, -1) \\ 64 & \cong & (1, -1, 1) \\ 34 & \cong & (1, -1, -1) \\ 71 & \cong & (-1, 1, 1) \\ 41 & \cong & (-1, 1, -1) \\ 29 & \cong & (-1, -1, 1) \\ 104 & \cong & (-1, -1, -1) \\ \end{array} $$

Отметим, что элемент $x\in \mathbb{Z}_m$ обратим тода и только тогда, когда $x$ взаимно прост с $m$, т.е. НОД$(x, m)=1$. Отсюда вытекает эквивалентное определение функции Эйлера.

Эквивалентное определение. Функция Эйлера $\phi(m)$ равна количеству чисел $x\in \mathbb{Z}$ в диапазоне $1 \le x \le m-1$, взаимно простых с $m$.

Обратимые элементы образуют группу в кольце $\mathbb{Z}_m$. Поэтому из теоремы Лагранжа вытекает теорема Эйлера (порядок элемента делит порядок группы).

Китайская теорема об остатках позволяет вычислить функцию Эйлера $\phi(m)$. Пусть $$ m = p_1^{e_1} p_2^{e_2} \dots p_k^{e_k}. $$ Изоморфизм колец в Китайской теореме об остатках индуцирует изоморфизм групп обратимых элементов. Обозначим через $U_m$ группу обратимых элементов кольца $\mathbb{Z}_m$. Имеем: $$ U_m \cong U_{p_1^{e_1}}\times U_{p_2^{e_2}}\times \dots \times U_{p_k^{e_k}}, $$ поэтому $$ \phi(m) = \phi(p_1^{e_1}) \phi(p_2^{e_2})\dots \phi(p_k^{e_k}). $$

Итак, достаточно вычислить функцию Эйлера для примарного кольца, т.е. кольца вида $\mathbb{Z}_{p^s}$, где $p$ — простое число. В таком кольце обратимы элементы, которые не делятся на $p$. Напротив, необратимы элементы, которые делятся на $p$, т.е. каждый $p$-й элемент, их число равно $$ |\mathbb{Z}_{p^s}|/p = p^s/p = p^{s-1}. $$ Отсюда число обратимых элементов равно числу всех элементов минус число необратимых: $$ \phi(p^s) = p^s - p^{s-1} = (p-1)p^{s-1}. $$

Окончательно получаем формулу Эйлера для числа $m$. Пусть $$ m = p_1^{e_1} p_2^{e_2} \dots p_k^{e_k}. $$ Тогда $$ \phi(m) = (p_1 - 1)p_1^{e_1 - 1} (p_2 - 1)p_2^{e_2 - 1} \dots (p_k - 1)p_k^{e_k - 1}. $$ В частности, когда число $m$ свободно от квадратов, т.е. $$ m = p_1 p_2 \dots p_k, $$ функция Эйлера записывается особенно просто: $$ \phi(m) = (p_1 - 1)(p_2 - 1)\dots (p_k - 1). $$

Для чисел $m$, свободных от крадратов, справедлив вариант теоремы Эйлера, аналогичный сформулированному выше варианту Малой теоремы Ферма.

Теорема 3' (частный случай теоремы Эйлера). Пусть положительное целое число $m$ свободно от квадратов, т.е. $$ m = p_1 p_2 \dots p_k, $$ где $p_i$ — различные простые числа, $i=1, 2, \dots, k$. Пусть $s\ge 0$ — произвольное неотрицательное целое число. Тогда для всякого $b\in \mathbb{Z}$ справедливо сравнение $$ b^{\phi(m)s + 1} \equiv b \pmod m. $$

Основные алгоритмы теории чисел

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

Дано: пара целых чисел $m$, $n$, не равных нулю одновременно.
Надо: вычислить наибольший общий делитель чисел $m$, $n$.

Запишем алгоритм на языке Python3. Отметим, что наибольший общий делитель определен с точностью до знака, для однозначности мы всегда будем возвращать положительное значение НОД.

def gcd(m, n):
    (a, b) = (m, n)
    while b != 0:
        (a, b) = (b, a%b)
    if a >= 0:
        return a
    else:
        return (-a)

Алгоритм основан на том, что при замене пары $(a, b)$ на пару $(b, r)$, где $r$ — остаток от деления $a$ на $b$, наибольший общий делитель пары сохраняется. Это следует из того, что множества всех общих делителей пар $(a, b)$ и $(b, r)$ совпадают. Последнее вытекает из равенств: $$ a = qb + r,\\ r = a - qb. $$ А наибольший общий делитель пары $(a, 0)$ равен $a$.

Расширенный алгоритм Евклида

Дано: пара целых чисел $m$, $n$, не равных нулю одновременно.
Надо: вычислить наибольший общий делитель $d$ чисел $m$, $n$ и найти его выражение в виде линейной комбинации этих чисел с целыми коэффициентами, т.е. найти целые числа $u$, $v$ такие, что $$ d = um + un. $$

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

Идея расширенного алгоритма Евклида заключается в том, что мы применяем обычный алгоритм Евклида, заменяя пару $(a, b)$ на $(b, r)$, где $r$ — остаток от деления $a$ на $b$. При этом в любой момент времени мы храним коэффициенты выражения чисел $a$, $b$ через исходные числа $m$, $n$: $$ a = u_1 m + v_1 n,\\ b = u_2 m + v_2 n. $$

При замене пары $(a, b)$ на пару $(b, r)$ мы перевычисляем коэффициенты $u_1, v_1, u_2, v2$. Пусть $q$, $r$ — целая часть частного и остаток от деления $a$ на $b$: $$ a = qb + r, \quad |r| < |b|. $$ Тогда $$ r = a - qb. $$ Подставляя выражения $a$ и $b$ в виде линейных комбинаций чисел $m$, $n$, получим: $$ r = a - qb = u_1 m + v_1 n - q(u_2 m + v_2 n) = \\ = (u_1 - q u_2)m + (v_1 - q v_2)n. $$ Таким образом, при замене $(a, b) \to (b, r)$ новые значения $u'_1, v'_1, u'_2, v'_2$ коэффициентов $u_1, v_1, u_2, v_2$ вычисляются по формулам $$ u'_1 = u_2, \quad v'_1 = v_2,\\ u'_2 = u_1 - q u_2, \quad v'_2 = v_1 - q v_2. $$

Запишем алгоритм на языке Python3. Алгоритм для входных параметров $m, n$ возвращает тройку чисел $(d, u, v)$ таких, что $d$ — наибольший общий делитель пары $(m, n)$ и $d = um + vn$.

def extEucl(m, n):
    (a, b) = (m, n)
    u1 = 1; v1 = 0
    u2 = 0; v2 = 1

    while b != 0:
        assert (a == u1*m + v1*n and b == u2*m + v2*n)

        q = a // b; r = a % b
        assert (a == q*b + r)

        (a, b) = (b, r)
        (u1, u2) = (u2, u1 - q*u2)
        (v1, v2) = (v2, v1 - q*v2)

    if a >= 0:
        return (a, u1, v1)
    else:
        return (-a, -u1, -v1)

Вычисление обратного элемента в кольце вычетов

Значение расширенного алгоритма Евклида состоит в том, что он позволяет вычислить обратный элемент в кольце вычетов по модулю $m$. Хорошо известно, что в этом кольце обратимы элементы, взаимно простые с $m$. Напомним, что элемент $x$ обратим в кольце, если существует $y$ такой, что $xy = 1$. Пусть $x$ обратим в кольце $\mathbb{Z}_m$, значит, он взаимно прост с $m$, т.е. $d =$ НОД$(x, y) = 1$. Тогда, применив расширенный алгоритм Евклида, можно вычислить числа $u, v$ такие, что $$ 1 = um + vx. $$ Но $um \equiv 0 \pmod{m}$, поэтому $$ 1 \equiv vx \pmod{m}. $$ Таким образом, $v$ является обратным элементом к $x$ в кольце $\mathbb{Z}_m$.

Запишем алгоритм на языке Python3. Алгоритм для входных параметров $x, m$ возвращает обратный к $x$ элемент, если $x$ обратим (т.е. целые числа $x, m$ взаимно простые); в противном случае возбуждается исключение типа ValueError.

def invmod(x, m):
    assert(m != 0)
    (d, u, v) = extEucl(m, x)
    if d == 1:
        return v%m
    else:
        raise ValueError("Element is not invertible modulo m")

Алгоритм быстрого возведения в степень — это еще одна жемчужина среди множества алгоритмов наряду с расширенным алгоритмом Евклида. Скорость работы этих алгоритмов (число арифметических операций) равна $O(\log n)$, где $n$ — максимум из аргументов; $\log n$ — это примерно число разрядов в двоичной записи числа $n$, т.е. длина записи входа. Понятно, что линейная зависимость времени работы алгоритма от длины входа — это очень хороший показатель (поскольку почти всегда требуется по крайней мере полностью прочесть входные данные). Алгоритм быстрого возведения в степень применим к любым объектам, на множестве которых определена ассоциативная операция умножения (числа, матрицы, кольца вычетов и т.д.). Особенно часто алгоритм быстрого возведения в степень применяется в кольцах вычетов $\mathbb{Z}_m$, поскольку, в отличие от целых чисел, результат не растет катастрофически быстро, можно, например, возводить число в 1000-значную степень практически мгновенно даже на самом медленном компьютере. Отметим также, что возведение в степень применяется очень часто, буквально на каждом шагу в криптографии.

Итак, пусть надо возвести элемент $a$ (принадлежащий в общем случае некоторой полугруппе, т.е. множеству с ассоциативной операцией умножения) в неотрицательную целую степень $n\in\mathbb{Z}$, $n\ge 0$. Идея алгоритма состоит в том, что мы строим цикл так, чтобы сохранялась величина $$ p b^k = a^n $$ (инвариант цикла). В цикле меняются 3 переменные $p, b, k$, причем указанный инвариант сохраняется при каждой итерации. Начальные значения переменных выбираются так, чтобы инвариант выполнялся: $$ p = 1,\ b = a,\ k = n. $$ Цель состоит в том, чтобы число $k$ уменьшалось. При $k=0$ из выполнения инварианта следует, что ответ содержится в переменной $p$: $$ a^n = p b^k, \ k = 0 \ \Rightarrow a^n = p. $$ Если $k$ четное, то $k$ делится пополам, а $b$ возводится в квадрат; при нечетном $k$ от $k$ вычитается единица, а $p$ домножается на $b$. Вот реализация алгоритма в случае кольца вычетов по модулю $m$ на языке Python3:

def powermod(a, n, m):
    p = 1; b = a; k = n
    while (k > 0):
        #Invariant: p*b^k == a^n
        if k%2 == 0:
            k //= 2
            b = b*b%m
        else:
            k -= 1
            p = p*b%m
    return p

Дано: набор взаимно простых целых чисел $m1, m2, \dots, m_k\in \mathbb{Z}$, НОД$(m_i, m_j) = 1$ при $i\ne j$, и набор произвольных целых чисел $r1, r2, \dots, r_k\in \mathbb{Z}$.
Надо: вычислить число $x\in \mathbb{Z}$ такое, что $$ x \equiv r_i \pmod{m_i},\ i=1, 2,\dots, k. $$

Рассмотрим сначала простейший случай $k=2$. Надо найти $x$ такой, что $$ x \equiv r_1 \pmod m_1\\ x \equiv r_2 \pmod m_2 $$ Первое сравнение справедливо тогда и только тогда, когда $$ x = r_1 + s m_1. $$ Подберем число $s\in\mathbb{Z}$ так, чтобы выполнялось второе сравнение. $$ r1 + s m_1 \equiv r_2 \pmod{m_2}. $$ Отсюда получаем $$ s m_1 \equiv r_2 - r_1 \pmod{m_2}. $$ Поскольку $m_1$ взаимно просто с $m_2$, то элемент $m_1$ обратим в кольце $\mathbb{Z}_{m_2}$. Пусть $v$ — обратный элемент к $m_1$ (он вычисляется с помощью рассмотренного выше алгоритма): $$ m_1 v \equiv 1 \pmod{m_2}. $$ Домножив сравнение выше справа на $v$, получим: $$ s m_1 v \equiv (r_2 - r_1)v \pmod{m_2}\\ s \equiv (r_2 - r_1)v \pmod{m_2}. $$ Таким образом, мы вычислили искомый элемент $s$.

В общем случае, для произвольного $k$ мы действуем по индукции. Для $k = 1$ ответом служит $x = r_1$. Пусть мы решили задачу для $k$, т.е. вычислили такое $x$, что $x\equiv r_i \pmod{m_i}$ при $i=1, 2,\dots, k$. Нужно "подправить" $x$ так, чтобы выполнялось и следующее сравнение $x\equiv r_{k+1} \pmod{m_{k+1}}$. Ищем новое значение $x'$ в виде $$ x' = x + s\cdot m_1 m_2\dots m_k, \quad s\in\mathbb{Z}. $$ Поскольку $m_1 m_2\dots m_k \equiv 0 \pmod{m_i}$ для $i = 1,\dots, k$, то $x'$, как и $x$, удовлетворяет всем предыдущим сравнениям. Вычислим множитель $s$. Поскольку числа $m_1 m_2\dots m_k$ и $m_{k+1}$ взаимно просты, то элемент $m_1 m_2\dots m_k$ обратим в кольце $\mathbb{Z}_{m_{k+1}}$. Вычислим обратный элемент $v$: $$ m_1 m_2\dots m_k \cdot v \equiv 1 \pmod{m_{k+1}}. $$ Из уравнения $$ x' \equiv r_{k+1} \pmod{m_{k+1}} $$ вычисляем искомое $s$ и $x'$: $$ x + s\cdot m_1 m_2\dots m_k \equiv r_{k+1} \pmod{m_{k+1}}\\ s\cdot m_1 m_2\dots m_k \equiv r_{k+1} - x \pmod{m_{k+1}}\\ s\cdot m_1 m_2\dots m_k\cdot v \equiv (r_{k+1} - x)v \pmod{m_{k+1}}\\ s \equiv (r_{k+1} - x)v \pmod{m_{k+1}}\\ x' = x + s\cdot m_1 m_2\dots m_k. $$

Программу на языке Python, реализующую данный алгоритм, предлагается написать самостоятельно в качестве упражнения.

Все указанные в тексте алгоритмы содержатся в файле numth.py

Один из интересных алгоритмов в теории чисел — это алгоритм извлечения квадратного корня в поле вычетов по модулю $p$, где $p$ — нечетное простое число. Постановка задачи:
пусть $p \ne 2$ — простое число, $a\in \mathbb{Z}_p$ — элемент поля $\mathbb{Z}_p$, являющийся квадратом. Требуется найти элемент $x\in \mathbb{Z}_p$ такой, что $$ x^2 \equiv a \pmod{p}. $$

Замечание. То, что элемент $a\in \mathbb{Z}_p$ является квадратом по модулю $p$, легко проверяется. Справедливо следующее утверждение.

Предложение. Пусть $p\in \mathbb{Z}$ — нечетное простое число. Элемент $a\in \mathbb{Z}$ является квадратом в поле $\mathbb{Z}_p$ тогда и только тогда, когда $$ a^{(p-1)/2} \equiv 1 \pmod{p}. $$

Предложение легко выводится из хорошо известного утверждения о том, что для любого конечного поля группа его ненулевых элементов по умножению является циклической. Пусть $b$ — порождающий элемент группы для поля $\mathbb{Z}_p$. Тогда любой ненулевой элемент $a$ поля $\mathbb{Z}_p$ представляется в виде $a = b^k$. Если $a$ — квадрат, то $k$ четное, $k=2s$, и $a^{(p-1)/2} = (b^{2s})^{(p - 1)/2} = (b^{p-1})^s \equiv 1^s = 1$ по Малой теореме Ферма. Обратно, пусть $a^{(p-1)/2} \equiv 1$. Имеем $a = b^k$, покажем, что $k$ четное. Допустим противное. Имеем $(b^k)^{(p-1)/2} = b^{k(p-1)/2} \equiv 1$, значит, поскольку порядок элемента $b$ равен $p - 1$, то $k(p-1)/2$ делится на $(p - 1)$. Но поскольку мы допустили, что $k$ нечетное, получаем противоречие, т.к. число $2$ входит в разложение числа $(p - 1)$ на простые множители в степени на единицу большей, чем в разложении $k(p-1)/2$.

Идея алгоритма. Рассмотрим два случая, в зависимости от остатка от деления числа $p$ на $4$. Поскольку $p$ нечетное, этот остаток может принимать два значения: $1$ и $3$ (или $-1$, поскольку $3\equiv -1 \pmod{4}$).

Случай 1: $p \equiv 3 \pmod{4}$. Это простой случай. Поскольку число $p + 1$ в этом случае делится на 4, то можем рассмотреть элемент $$ x \equiv a^{(p + 1)/4} \pmod{p}, $$ он и дает искомое решение нашей задачи, поскольку $$ x^2 = a^{(p + 1)/2} = a^{(p - 1)/2} \cdot a \equiv a \pmod{p}. $$

Случай 2: $p \equiv 1 \pmod{4}$. Это сложный случай. Решение, предложенное в предыдущем случае, уже не годится, т.к. число $p + 1$ не делится на 4. Так как $a$ является квадратом по модулю $p$, то $$ a^{(p - 1)/2} \equiv 1 \pmod{p}. $$ Представим число $(p - 1)/2$ в виде $$ (p - 1)/2 = q\cdot 2^s, $$ где $q$ — нечетное число, $s > 0$. Рассмотрим следующую последовательность элементов поля $\mathbb{Z}_p$: $$ a^q, a^{2q}, a^{4q}, \dots, a^{q\cdot 2^s} = a^{(p-1)/2}. $$ В этой последовательности каждый следующий элемент является квадратом предыдущего по модулю $p$. Так как $a$ является квадратом, то эта последовательность заканчивается единицей. Заметим, что в любом поле есть только два квадратных корня из минус единицы: $1$ и $-1$. Поэтому последовательность имеет один из двух видов, в зависимости от того, начинается ли она с единицы или нет: $$ 1, 1, 1, \dots, 1 $$ или $$ *, *, \dots, *, -1, 1, \dots, 1. $$

В первом случае имеем $a^q \equiv 1 \pmod{p}$ и $q$ — нечетное число. Тогда $$ r \equiv a^{(q + 1)/2} \pmod{p} $$ является искомым квадратным корнем из $a$, поскольку $$ r^2 = a^{q + 1} = a^q\cdot a \equiv a \pmod{p}. $$

Осталось справиться со вторым случаем. В нем последовательность не начинается с единицы и число $$ r \equiv a^{(q + 1)/2} \pmod{p} $$ сразу не является решением нашего уравнения $x^2 \equiv a \pmod{p}$. Обозначим $$ t \equiv a^q \pmod{p}. $$ Тогда выполняется равенство $$ r^2 \equiv t\cdot a \pmod{p}, $$ причем последовательность квадратов, порождаемая элементом $t$, $$ t, t^2, t^4, ..., t^{2^s} $$ имеет вид $$ *, \dots, *, -1, 1, ..., 1 $$ Идея состоит в том, чтобы так подправить числа $r$ и $t$, чтобы при сохранении этих двух условий число единиц в конце последовательности квадратов, задаваемых элементом $t$, увеличилось бы. Таких шагов может быть несколько, после каждого шага число единиц в конце последовательности будет увеличиваться по крайней мере на единицу. Тогда в конце концов $t$ станет равным единице, и из выполнения условия $r^2 \equiv t\cdot a \pmod{p}$ будет следовать, что $r$ является искомым корнем из $a$.

Для этого найдем произвольный элемент $z\in \mathbb{Z}_p$, не являющийся квадратом. Таких элементов очень много (каждый второй ненулевой элемент поля), поэтому можно просто перебирать подряд элементы $\mathbb{Z}_p$ и проверять, являются ли они квадратом. Пусть $z\in \mathbb{Z}_p$ не является квадратом, т.е. $z^{(p-1)/2} \equiv -1 \pmod{p}$. Тогда последовательность $$ z^q, z^{2q}, z^{4q}, \dots, z^{q\cdot 2^s} = z^{(p-1)/2}. $$ имеет вид $$ *, *, *, \dots, *, -1, $$ где звездочкой обозначены элементы, отличные от $\pm 1$. В последовательности квадратов, порожденной элементом $t$, минус единица стоит строго левее минус единицы в последовательноси, порожденной $z^q$. Например: $$ \begin{array}{lllllllll} t: &\underline{*}, &*, &*, &\dots, &*, &-1, &1, &1 \\ z^q: &*, &*, &\underline{*}, &\dots, &*, &*, &*, &-1 \\ \end{array} $$ В данном случае в первой последовательности минус единица стоит левее, чем во второй, на 2 позиции.

Теперь домножим $t$ на третий элемент последовательности, порожденной $z^q$, в формуле выше это элемент отмечен подчеркиванием. Элемент $r$ домножим на предыдущий элемент из последовательности для $z^q$ — это квадратный корень из подчеркнутого элемента. Обозначим подчеркнутый элемент через $b^2$, а предыдущий элемент через $b$. Тогда новые значения переменных $r$, $t$ будут $$ r' \equiv r\cdot b,\qquad t' \equiv t\cdot b^2 \pmod{p}. $$ Выполнение равенста $r'^2 \equiv t'\cdot a \pmod{p}$ сохранится, поскольку $$ r'^2 \equiv (rb)^2 \equiv r^2 b^2 \equiv ta b^2 \equiv t'a \pmod{p}. $$ Но последовательность квадратов для элемента $t'$ уже будет содержать больше единиц в конце, поскольку она получается перемножением последовательности для $t$ и последовательности для $z^q$, сдвинутой на 2 позиции влево. В результате минус единица первой последовательности умножится на минус единицу во второй последовательности, что даст искомую дополнительную единицу.

Для иллюстрации этого алгоритма рассмотрим следующий пример. Пусть $p = 1009$ (это простое число) и $a = 71$. Нужно найти число $x$ такое, что $x^2 \equiv 71 \pmod{1009}$. Имеем $$ 1009 \equiv 1 \pmod{4}, $$ значит, нужно рассмотреть второй, более сложный случай. Выделим из числа $p - 1$ максимальную степень двойки: $$ p - 1 = 1008 = 63\cdot 2^4,\\ q = 63,\ s = 4. $$ Вычислим $t \equiv a^q \pmod{p}$ и последовательность квадратов, порожденную элементом $t$: $$ t\equiv 247,\ t^2\equiv 469,\ t^4\equiv -1,\ t^8\equiv 1 \pmod{1009}, \\ r \equiv a^{(q + 1)/2} \equiv 383 \pmod{1009}, \\ r^2 \equiv t\cdot a \pmod{1009}. $$ Найдем элемент $z$ — не квадрат по модулю $p$, и вычислим последовательность квадратов, начинающуюся с $z^q$: $$ z \equiv 11 \pmod{1009},\\ z^q \equiv 179,\ z^{2q} \equiv 762, \ z^{4q} \equiv 469, z^{8q} \equiv -1 \pmod{1009}. $$ Подправим теперь элементы $t$, $r$: поскольку в $z$-последовательности элемент $-1$ стоит на 4-м месте, а в $t$-последовательности на 3-м, что означает сдвиг на одну позицию, то $t$ надо домножить на второй элемент $z$-последовательности 762, а $r$ — на первый элемент 179 (корень из второго). Получим $$ t' \equiv t\cdot 762 \equiv 540 \pmod{1009}, \\ r' \equiv r\cdot 179 \equiv 954 \pmod{1009}. $$ Новая последовательность квадратов для нового значения $t = t'$ будет $$ t\equiv 540,\ t^2\equiv -1,\ t^4\equiv 1,\ t^8\equiv 1 \pmod{1009}, $$ новое значение $r = r' = 954$. Теперь элемент $-1$ в $t$-последовательности стоит на втором месте, а в $z$-последовательности на четвертом, что означает сдвиг на две позиции; поэтому $t$ надо домножить на третий элемент $z$-последовательности 469, а $r$ — на второй элемент 762: $$ t' \equiv t\cdot 469 \equiv 1 \pmod{1009}, \\ r' \equiv r\cdot 762 \equiv 468 \pmod{1009}. $$ Поскольку новое значение $t' = 1$, то алгоритм на этом заканчивается, и последнее значение $r' = 468$ будет являться искомым квадратным корнем из $a = 71$.

Малая теорема Ферма является необходимым условием простоты числа $p$. Долгое время существовало заблуждение, что и обратная теорема верна: если для всех $b\in \mathbb{Z}$, взаимно простых с $p$, выполняется сравнение $$ b^{p-1} \equiv 1 \pmod{p}, $$ то $p$ — простое число. Причем это предполагалось верным даже для конкретного основания $b = 2$: $$ 2^{p-1} \equiv 1 \pmod{p} \Rightarrow p\ \mbox{простое}. $$ Первым показал, что это не так, французский математик П.Ф.Саррус в 1819 г. (спустя почти 200 лет после Ферма!): число $341 = 31\cdot 11$ не является простым, но удовлетворяет заключению Малой теоремы Ферма для основания 2. Действительно, $$ 2^{340} = (2^{10})^{34} = 1024^{34} = (3\cdot 341 + 1)^{34} \equiv\\ \equiv 1^{34} \equiv 1 \pmod{341}. $$ Подобные числа называют псевдопростыми по основанию $b$ (в данном случае $b=2$): они не являются простыми, но ведут себя как простые в Малой теореме Ферма. Минимальное подобное число — это $91=7\cdot 13$, оно является псевдопростым по основанию 3. Действительно, $$ 3^{90} = (3^6)^{15} = 729^{15} = (8*91 + 1)^{15} \equiv\\ \equiv 1^{15} \equiv 1 \pmod{91}. $$

Отметим, что, выбрав другие основания, можно показать, что числа 341 и 91 не являются простыми, применив Малую теорему Ферма: $$ 3^{340} \equiv 56 \pmod{341},\quad 2^{90} \equiv 64 \pmod{91}. $$ Однако существуют такие парадоксальные числа $m$, которые не являются простыми, но ведут себя как простые в Малой теореме Ферма для всех оснований $b$, взаимно простых с $m$. Такие числа называются кармайкловыми в честь американского математика Роберта Кармайкла, открывшего первое такое число 561 в 1910 г. Кармайклово число $m$ обладает следующими свойствами:

1) $m$ — составное число;
2) для всякого основания $b\in \mathbb{Z}$, взаимно простого с $m$, выполняется сравнение $b^{m-1} \equiv 1 \pmod{m}$.

Справедлива следующая теорема, описывающая устройство кармайкловых чисел.

Теорема 4 (о кармайкловых числах). Число $m\in\mathbb{Z}$ кармайклово тогда и только тогда, когда
1) $m = p_1 p_2 \dots p_k$, где $k\ge 3$ и $p_i$ — простые числа (т.е. число $m$ свободно от квадратов);
2) $(p_i - 1) | (m - 1)$ для любого $i=1, 2, \dots, k$.

Теорема легко выводится из Китайской теоремы об остатках и Малой теоремы Ферма.

Минимальные кармайкловы числа — это $561, 1105, 1729, \dots$ Проверим, например, что число 561 удовлетворяет условию теоремы 4. Действительно, $561 = 3\cdot 11\cdot 17$ и $(3 - 1)|560$, $(11 - 1)|560$, $(17 - 1)|560$. Совсем недавно доказано, что множество кармайкловых чисел бесконечно (Алфильд, Гранвиль, Померанс, 1994 г.).

Хотя такие числа не простые, это невозможно доказать, используя в качестве теста Малую теорему Ферма. Разве что нам повезет и мы выберем в качестве основания $b$ число, не взаимно простое с $m$, но для больших чисел вероятность подобного события близка к нулю. Так что Малая теорема Ферма не годится в качестве теста простоты числа. Тем не менее тест Ферма можно немного подправить, чтобы получить вероятностный тест простоты числа.

В тесте простоты Миллера-Рабина проверяется простота числа $m$. Для этого выбирается случайное основание $b\in \mathbb{Z}$, $2\le b < m-1$. Тест для фиксированного основания $b$ выдает один из двух ответов:
(1) число $m$ составное;
(2) не знаю.

Если тест выдает ответ (1), то он приводит и доказательство того, что число $m$ составное. Вероятность ответа (2) для составного числа не превышает $1/4$, так что в этом случае число $m$ скорее всего простое. Если тест 20 раз подряд выдал ответ "не знаю", то вероятность того, что число $m$ составное, не превышает $(1/4)^{20} = 2^{-40} \approx 10^{-12}$, т.е. совершенно ничтожна (и делая большее число тестов, ее можно сделать сколь угодно малой). Тем не менее тест не выдает доказательства, что число $m$ простое. Однако для практических целей криптографии такие числа вполне приемлемы.

В тесте Миллера-Рабина проверяется нечетное число $m$ (поскольку все простые числа, кроме 2, нечетные). Пусть $$ m - 1 = t \cdot 2^s,\ t \text{ — нечетное число}. $$ Выберем случайное основание $b\in \mathbb{Z}_m$, $2 \le b < m-1$. Применив к $b$ алгоритм быстрого возведения в степень, вычислим последовательность $$ x_0 \equiv b^t,\ x_1 \equiv x_0^2,\ x_2 \equiv x_1^2, \dots, x_s \equiv x_{s-1}^2 \equiv b^{m-1} \pmod{m}. $$ В последовательности $$ x_0, x_1, x_2,\dots, x_s $$ каждый следующий элемент является квадратом предыдущего.

Тест выдает ответ "$m$ — составное число" в двух случаях:
1) последний элемент последовательности не равен 1 — в этом случае не выполняется Малая теорема Ферма;
2) в последовательности есть фрагмент $x_i$, $x_{i+1} \equiv x_i^2$, где $$ x_i \not\equiv \pm 1,\quad x_{i+1} \equiv 1 \pmod{m}. $$ В этом случае также $m$ является составным числом, поскольку, если бы $m$ было простым, то $\mathbb{Z}_m$ было бы полем, а в поле нет других квадратных корней из единицы, кроме $\pm 1$.

В противном случае тест выдает ответ "не знаю", в этом случае $m$ скорее всего является простым числом.

Итак, ответ "$m$ составное" выдается в случаях, когда последовательность имеет один из трех видов: $$ *, *, *, \dots, *, * \\ *, *, *, \dots, *, -1 \\ *, *, \dots, *, 1, 1, \dots, 1 $$ где звездочками обозначены элементы, отличные от $\pm 1$.

Ответ "не знаю" (когда число $m$, скорее всего, является простым) выдается в случаях, когда последовательность имеет один из двух видов: $$ 1, 1, \dots, 1 \\ *, *, \dots, *, -1, 1, \dots, 1. $$

Нетривиальная часть результата Миллера-Рабина состоит в том, что вероятность ответа "не знаю" для составного числа $m$ не превышает $1/4$ (впрочем, доказательство этого факта несложное, оно основано на тривиальном соображении, что число элементов любой собственной подгруппы не превышает половины числа элементов группы). Таким образом, выполняя несколько тестов для разных случайных оснований $b$, можно сделать вероятность ошибки сколь угодно малой. Например, если тест для числа $m$ для разных случайных оснований 20 раз выдал ответ "не знаю", то число $m$ является, скорее всего, простым, вероятность того, что оно составное, не превышает $(1/4)^{20} \approx 10^{-12}$.

Запишем алгоритм теста Миллера-Рабина на псевдокоде. Реализация алгоритма на языке Python оставлена читателю в качестве упражнения. Отметим, что в Python очень удобно организована работа со случайными числами. Для получения случайного целого числа $x$ в заданном диапазоне $a \le x < b$ используется функция random.randrange(a, b) модуля random:

import random
. . .
    x = randow.randrange(a, b)

Входными параметрами алгоритма являются целое число $m$ и количество случайных тестов $k$, по умолчанию $k=20$. Алгоритм выдает ответ True, если число $m$, скорее всего, простое, и False, если $m$ составное (в этом случае алгоритм получает доказательство этого факта).

алгоритм primeTest(m, k=20)
| если m принадлежит списку [2, 3, 5, 7, 11, 13, 17, 19]
| | то вернуть True
| иначе если m четное или m < 20
| | то вернуть False
| конец если
|
| t = m - 1; s = 0
| цикл пока t четное
| | t = t/2; s = s + 1
| конец цикла
| утверждение: m - 1 == t*2^s
|
| цикл k раз выполнить
| | b = случайное целое число в диапазоне 2 <= b < m-1
| | x = b^t (mod m)
| | если x == 1 (mod m) или x == -1 (mod m)
| | | continue (к следующей итерации цикла)
| | конец если
| | i = 0
| | цикл пока i < s
| | | x = x*x (mod m)
| | | i = i + 1
| | | если x == 1 (mod m)
| | | | # нашли нетривиальный квадратный корень из 1
| | | | вернуть False
| | | иначе если x == -1 (mod m)
| | | | break (выйти из цикла)
| | | конец если
| | конец цикла
| | если i == s и x != 1 (mod m)
| | | # не выполняется Малая теорема Ферма
| | | вернуть False
| | конец если
| конец цикла
|
| вернуть True
конец алгоритма

Схема кодирования с открытым ключом совершила революцию в криптографии. Все классические криптографические алгоритмы были симметричными, т.е. один и тот же ключ использовался как для шифрования, так и для расшифровки. Кроме того, и сами алгоритмы шифрования имели статус секретной или служебной информации. В отличие от них, сам по себе алгоритм шифрования с открытым ключом не представляет никакого секрета, он очень простой и основан на сведениях из теории чисел, восходящих еще к Древней Греции и к работам Ферма. При кодировании с открытым ключом ключ шифрования, так же как и алгоритм шифрования, известен всем, секрет представляет лишь ключ дешифровки. Причем этот ключ невозможно вычислить по ключу шифрования существующими алгоритмами (пока не реализован квантовый компьютер). Вычисление секретного ключа сводится к хорошо известной математической проблеме — разложению на множители большого целого числа. Эту проблему, конечно же, не способен решить никакой хакер, подобный результат может быть достигнуть только усилиями всего человечества, и его невозможно скрыть в случае успеха (который, скорее всего, никогда не будет достигнут).

Отметим, что схема RSA была изобретена в MIT (Massachusetts Institute of Technology) тремя авторами — R.Rivest, A.Shamir, L.Adleman, она названа по первым буквам их фамилий; хотя идея асимметричного кодирования высказывалась и раньше другими учеными. Это типично научное достижение, опубликованное открыто, в отличие от многих работ в области криптографии, выполняемых с той или иной степенью секретности.

В схеме RSA используется целое число $m$, которое равно произведению двух больших простых чисел $p$ и $q$ (длины чисел $p$ и $q$ не меньше 150 десятичных цифр): $$ m = pq,\quad p, q\in \mathbb{Z} \text{ — простые числа}. $$ При этом $m$ известно всем, но разложение $m$ на множители (числа $p$ и $q$) секретное.

Для генерации пары ключей вычисляется функция Эйлера $$ \phi(m) = (p - 1)(q - 1). $$ В кольце вычетов $\mathbb{Z}_{\phi(m)}$ выбирается случайный элемент $e$, взаимно простой с $\phi(m)$. С помощью расширенного алгоритма Евклида вычисляется обратный к $e$ элемент $d$ в кольце $\mathbb{Z}_{\phi(m)}$: $$ e\cdot d \equiv 1 \pmod{\phi(m)}. $$

Пара $(m, e)$ является открытым ключом, пара $(m, d)$ — секретным. (Секретным является число $d$.)

Напомним, что в современной криптографии используется идея блочного шифрования: любая передаваемая информация сначала представляется в виде последовательности нулей и единиц (битов), эта последовательность делится на блоки фиксированной длины, блоки шифруются и передаются по отдельности. Каждый блок можно рассматривать как двоичную запись некоторого числа, причем эти числа ограничены из-за ограниченности длины блока. Значит, эти числа можно считать элементами кольца вычетов по фиксированному числу $m$ (его величина зависит от размера блока).

Шифрование — это отображение $$ E: \mathbb{Z}_m \to \mathbb{Z}_m, $$ которое по исходному сообщению $x\in \mathbb{Z}_m$ вычисляет зашифрованное сообщение $y = E(x)$. (Мы используем букву $E$ как первую букву слова Encrytion.) Отображение $E$ состоит в возведении элемента $x$ в степень $e$ в кольце вычетов по модулю $m$: $$ E: x \mapsto x^e \pmod m. $$ Расшифровка — это отображение $D$ (от слова Decryption), которое элемент $y$ возводит в степень $d$: $$ D: y \mapsto y^d \pmod m. $$ Фактически $D$ извлекает корень степени $e$ в кольце $\mathbb{Z}_m$ — оказывается, что для извлечения корня в этом кольце нужно возвести элемент в степень, обратную к $e$ в кольце вычетов по модулю функции Эйлера $\mathbb{Z}_{\phi(m)}$.

Докажем, что отбражения $E$ и $D$ взаимно обратные. Для любого $x\in \mathbb{Z}_m$ имеем $$ y = E(x) \equiv x^e \pmod{m} \\ D(y) \equiv (x^e)^d \equiv x^{ed} \pmod{m}. $$ Поскольку $ed \equiv 1 \pmod{\phi(m)}$, то $$ ed = 1 + s\cdot \phi(m), \quad s\in \mathbb{Z},\ s \ge 0. $$ По Теореме 3' (частный случай теоремы Эйлера), поскольку $m=pq$ свободно от квадратов, то для любого $x\in \mathbb{Z}_m$ $$ x^{ed} \equiv x^{1 + s\cdot \phi(m)} \equiv x \pmod{m}. $$ Мы показали, что отображение $E$ имеет левый обратный элемент (композиция $E$ и $D$ есть тождественное отображение в $\mathbb{Z}_m$): $$ % D\circ E = \mathbb{1}_{\mathbb{Z}_m}. D(E(x)) = x\quad \text{для любого } x\in \mathbb{Z}_m. $$ Значит, отображение $E$ инъективно и, поскольку множество $\mathbb{Z}_m$ конечно, то и сюръективно; поэтому $E$ и $D$ являются взаимно обратными биективными отображениями.

Генерация и взлом ключей в схеме RSA

Для генерации пары ключей — открытого и секретного — надо выполнить следующие шаги.

1. Сгенерировать два больших простых числа $p$ и $q$. Величина их должна быть такой, чтобы произведение $m = pq$ заведомо невозможно было бы разложить на множители существующими алгоритмами, т.е. в настоящее время каждое из этих чисел должно иметь длину не меньше 150 десятичных цифр.
2. Вычисляется функция Эйлера $\phi(m) = (p - 1)(q - 1)$.
3. Выбирается случайное число $e\in \mathbb{Z}_{\phi(m)}$, взаимно простое с $\phi(m)$, и с помощью расширенного алгоритма Евклида вычисляется обратный к $e$ элемент $d\in \mathbb{Z}_{\phi(m)}$: $$ ed \equiv 1 \pmod{\phi(m)}. $$
4. Открытый ключ — это пара $(m, e)$, секретный ключ — пара $(m, d)$. Шифрование — это отображение $$ E: x \mapsto x^e \pmod m, $$ дешифровка — отображение $$ D: y \mapsto y^d \pmod m. $$

Итак, для генерации пары ключей надо уметь генерировать большие простые числа. Это можно сделать, используя вероятностный тест простоты Миллера-Рабина. Дела в том, что простых чисел очень много. Количество простых чисел, не превосходящих $x$, в математике обозначается через $\pi(x)$, известно, что эта функция эквивалентна $$ \pi(x) \sim x/\ln x \quad \text{при } x\to\infty. $$ Грубо говоря, это означает, что для больших чисел $x$ — например, 200-значных десятичных — в отрезке длиной $\ln x$ содержится в среднем одно простое число. (Для 200-значных чисел $\ln x \approx 460.5$, это значит, что в среднем каждое 461-е число простое.) Можно просто генерировать абсолютно случайные нечетные числа и проверять их на простоту, довольно быстро таким образом найдется простое число. Конечно, с помощью теста Миллера-Рабина мы не получаем доказательства простоты, но вероятность ошибки ничтожна (может быть сделана сколь угодно малой), так что такие числа вполне можно использовать на практике. (Отметим, что существует и способ генерации больших простых чисел с сертификатами, т.е. легко проверяемыми доказательствами простоты, но на практике чаще используются вероятностные алгоритмы.)

Реализация алгоритма генерации RSA-ключей на Python'е

Приведем реализацию алгоритма генерации RSA-ключей на Python'е. Функция generateRSAKeys(keyLength) генерирует пару ключей (открытый и секретный) заданной длины в битах. Результатом работы этой функции является тройка $$ (m, e, d),\quad \text{где }m = pq,\quad ed\equiv 1 \pmod{\phi(m)}, $$ и число $m$ имеет заданную длину keyLength его двоичной записи. Для реализации функции мы воспользовались генератором случайных чисел

import random
. . .
    p = random.randrange(k0, k1)

    primeTest(p)

    (_, d) = invmod(e, phi)

Вот полный текст алгоритма:

import random

def generateRSAKeys(keyLength):
    '''Generate RSA public and secret keys
    Input: key length in bits.
    Return: a triple (m, e, d), where
        m is a product of two prime numbers of length keyLength/2,
        (e, d) is a pair of mutually inverse numbers modulo phi(m).
    The pair (m, e) is a public key,
    the pair (m, d) is a secret key.'''

    # Generate 2 primes of length keyLength/2
    k0 = 2**(keyLength//2 - 1); k1 = 2*k0
    k0 += 1
    primes = []
    steps = 0
    while len(primes) < 2 and steps < 10000:
        p = random.randrange(k0, k1)
        if primeTest(p):
            primes.append(p)
        steps += 1
    if len(primes) < 2:
        return None     # Falure...

    m = primes[0]*primes[1]
    phi = (primes[0] - 1)*(primes[1] - 1)

    # Generate numbers e and d such that e d = 1 (mod phi(m))
    generated = False; steps = 0
    while not generated and steps < 10000:
        e = random.randrange(2, phi-1)
        if gcd(e, phi) != 1:
            steps += 1; continue
        else:
            generated = True
    if not generated:
        return None     # Falure...

    (_, d) = invmod(e, phi)
    return (m, e, d)

Проверим работу алгоритма. Сгенерируем пару ключей длины 20 десятичных цифр (66 двоичных).

>>> keys = generateRSAKeys(66)
>>> keys
(31878197109862343987, 18964933914993953353, 23894192283802500217)
>>> (m, e, d) = keys
>>> x = 123456789
>>> y = pow(x, e, m)
>>> y
8514718178434504648
>>> pow(y, d, m)
123456789

Конечно, в данном примере длина ключа в 20 десятичных цифр совершенно не достаточна для надежного шифрования — при современном состоянии компьютерной науки длина ключа должна быть не меньше 300 десятичных цифр (а лучше еще в 4 раза больше, т.е. 4096 двоичных цифр).

Взлом ключа в схеме RSA

Идея любой схемы асимметричного кодирования состоит в том, что по открытому ключу невозможно вычислить секретный ключ — вернее, существующие алгоритмы потребуют чудовищного времени для этого, миллиарды миллиардов лет — ничто по сравнению с ним. Разве что будет действительно реализован квантовый компьютер (если это произойдет, то это будет означать крах существующих банковских систем, электронных подписей, электронных валют, секретных каналов передачи, систем "свой-чужой" и т.п. — серьезные проблемы для всего человечества).

Для того, чтобы в схеме RSA вычислить секретный ключ по открытому, нужно вычислить функцию Эйлера $$ \phi(m) = (p - 1)(q - 1). $$ Зная функцию Эйлера, мы по отрытому ключу $(m, e)$ легко вычислим секретный ключ $(m, d)$, поскольку $d$ является обратным элементом к $e$ в кольце $\mathbb{Z}_{\phi(m)}$, а обратный элемент легко вычисляется с помощью расширенного алгоритма Евклида.

Для вычисления функции Эйлера нажно знать разложение числа $m$ на множители, которое держится в секрете. Верно и обратное: зная функцию Эйлера $\phi(m)$, можно вычислить разложение числа $m$ на множители. Действительно, равенства $$ \begin{cases} p\cdot q = m \\ (p - 1)(q - 1) = \phi(m) \end{cases} $$ можно рассматривать как систему из двух уравнений с двумя неизвестными $(p, q)$, которая легко решается: $$ \begin{array}{l} p = \frac{m + \phi(m) + \sqrt{m^2 - 2m(\phi(m) + 1) + % (\phi(m) - 1)^2} - 1}{% m - \phi(m) + \sqrt{m^2 - 2m(\phi(m) + 1) + (\phi(m) - 1)^2} - 1} \\ q = m/p \end{array} $$

Итак, разложение числа $m$ на множители и вычисление функции Эйлера $\phi(m)$ — это две эквивалентные по трудности задачи, каждая из которых сводится к другой. Поэтому задача взлома секретного ключа эквивалентна задаче разложения числа $m$ на множители.

Схема кодирования с отрытым ключом RSA дает возможность удостоверить отправителя сообщения. Для этого отправитель шифрует сообщение с помощью своего секретного ключа. Воспользовавшись открытым ключом отправителя, каждый может расшифровать его сообщение, поскольку отображения $E$ и $D$ взаимно обратные и не только $D\circ E$, но и $E\circ D$ — тождественное отображение: $$ E(D(x)) = x\quad \text{для любого } x\in \mathbb{Z}_m. $$ Сформировать сообщение может только обладатель секретного ключа, таким образом подтверждается его авторство.

Подобное "обратное" применение схемы кодирования с открытым ключом даже более важно, чем прямое. В частности, оно дает возможность сформировать электронную подпись документа, которая исключает возможность его подмены. Сам документ передается открыто, вместе с ним передается электронная подпись, которую все могут проверить, воспользовавшись открытым ключом отправителя документа. Электронная подпись формируется слудеющим образом. По документу вычисляются следующие его характеристики.

1. Размер документа в байтах.
2. Хэш-функция от содержимого документа.
3. Дата и время создания документа.
4. К этим данным добавляется номер транзакции, можно также добавить и другие дополнительные данные.

Вся эта информация шифруется с помощью секретного ключа отправителя. Зашифрованная информация представляет собой электронную подпись документа.

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

Для подмены документа нужно сгенерировать фальшивый документ той же длины и с тем же значением хэш-функции, а эта задача практически неразрешимая (хотя анализ ее трудности непрост, она сильно зависит от свойств используемой хэш-функции).

Термин "факторизация" используется для алгоритмов разложения на множители (от англ. factor — множитель). Для взлома ключа RSA требуется разложить на множители число $m = p\cdot q$, где $p$, $q$ — большие простые числа. Изучение подобных алгоритмов нужно не только преступникам: создатели криптосистем должны быть уверены в их стойкости, для этого надо хорошо представлять пределы применимости современных алгоритмов факторизации. И, конечно, задача факторизации интересна и сама по себе, с чисто научной точки зрения.

Рекордсменом среди алгоритмов факторизации чисел до недавнего времени являлся алгоритм "квадратичное решето". Он основан на совсем простом соображении, восходящем еще к Эйлеру. Если мы найдем два числа $a, b\in \mathbb{Z}_m$, квадраты которых совпадают по модулю $m$, то $m$ раскладывается на множители. $$ a^2 \equiv b^2 \pmod{m}, \quad a\not\equiv \pm b \pmod{m}. $$ Отсюда следует, что $$ 0 \equiv a^2 - b^2 \equiv (a - b)(a + b) \pmod{m}. $$ Это означает, что $$ m | (a - b)(a + b). $$ Но, поскольку $a\not\equiv \pm b \pmod{m}$, то $$ m \nmid (a - b),\quad m \nmid (a + b). $$ Это означает, что $$ \text{НОД}(m, a - b)\ \text{ и }\ \text{НОД}(m, a + b) $$ являются нетривиальными делителями числа $m$.

В частном случае, если мы найдем нетривиальный квадратный корень из единицы в кольце $\mathbb{Z}_m$, то $m$ раскладывается на множители: $$ a^2 \equiv 1,\ a \not \equiv \pm 1 \pmod m \Rightarrow \\ \text{НОД}(m, a - 1)\ \text{— нетривиальный делитель } m. $$

В алгоритме "квадратичное решето" осуществляется поиск двух квадратов, совпадающих по модулю $m$. Для этого генерируются квадраты по модулю $m$ и эти числа раскладываются на множители по мультипликативной базе $\cal B$, состоящей из набора небольших простых чисел. Пусть $|\cal B| = k$. Если мы найдем $k+1$ квадратов, разложимых по базе $\cal B$, то степени простых множителей из $\cal B$ в разложении этих квадратов образуют матрицу размера $(k+1)\times k$. Строки ее линейно зависимы над полем $Z_2$. Это означает, что сумма некоторого подмножества строк этой матрицы есть строка с четными числами. Значит, из произведения соответствующих квадратов можно извлечь квадратный корень двумя способами — во-первых, взяв произведение исходных чисел, возводимых в квадрат, во-вторых, взяв произведение простых чисел из базы $\cal B$ в степенях, равных половинам элементов указанной суммы строк матрицы. Таким образом получается нетривиальное равенство $$ a^2 \equiv b^2 \pmod{m} $$ и $m$ раскладывается на множители. Идея этого алгоритма принадлежит Дж.Диксону (John D. Dixon, Carleton University), работа которого была опубликована в 1981 г.

В настоящее время наиболее сильным из известных алгоритмов является общий метод решета над конечным расширением поля рациональных чисел (General Number Field Sieve), обобщающий классический метод квадратичного решета. Мы, однако, здесь не будем рассматривать эти достаточно сложные алгоритмы факторизации чисел и ограничимся совсем простыми, но при этом очень красивыми методами факторизации Полларда.

Алгоритмы факторизации Полларда часто называют методами типа Монте-Карло, поскольку положительный результат зависит от везения и достигается не всегда. Достоинством этих алгоритмов является то, что их реализация очень простая и они хорошо извлекают небольшие простые множители даже из очень больших целых чисел. Имеется два алгоритма целочисленной факторизации Полларда: $\rho$-алгоритм и $p-1$-алгоритм.

$\rho$-алгоритм факторизации Полларда основан на поиске цикла в рекуррентной последовательности. Пусть $m$ — составное число, которое мы хотим разложить на множители. Пусть $p|m$ — минимальный простой делитель числа $m$ ($m$ нам известно, $p$ нет). Расcмотрим следующее отображение $$ f: Z_m \to Z_m, \quad f(x) \equiv x^2 + 1 \pmod{m} $$ (в качестве $f(x)$ можно взять любой многочлен степени больше 1). Выберем случайное начальное число $x_0\in \mathbb{Z}_m$ и рассмотрим бесконечную последовательность $$ x_0,\ x_1=f(x_0),\ x_2=f(x_1),\ x_3=f(x_2),\dots $$ В последовательности $$ x_0, x_1, x_2, x_3, \dots $$ каждое следующее число является образом предыдущего под действием отображения $f$. Поскольку все элементы $x_i$, $i=0, 1, 2,\dots$ принадлежат конечному множеству $\mathbb{Z}_m$, эта последовательность с некоторого момента становится циклической. Рассмотрим ту же самую последовательность, но уже по модулю $p$, а не $m$: $$ \overline{x}_0, \overline{x}_1, \overline{x}_2, \overline{x}_3, \dots \quad \overline{x}_i \equiv x_i \pmod p. $$ Эта последовательность также является циклической. Но, поскольку $p < m$, то с большой вероятностью длина цикла в последовательности $\{\overline{x}_i,\ i=0,1,2,\dots\}$ меньше, чем длина цикла в последовательности $\{x_i,\ i=0,1,2,\dots\}$. Значит, найдутся индексы $i < j$ такие, что $$ \overline{x}_i = \overline{x}_j, \quad x_i \ne x_j $$ (в последовательности $\{\overline{x}_i\}$ мы уже вошли в цикл, в последовательности $\{x_i\}$ еще нет). Это означает, что $$ x_i \equiv x_j \pmod{p}, \quad x_i \not\equiv x_j \pmod{m}. $$ Отсюда вытекает, что $$ p | (x_i - x_j), \quad m \nmid (x_i - x_j). $$ Значит, $$ p | \text{НОД}(x_i - x_j, m),\quad m \nmid \text{НОД}(x_i - x_j, m), $$ то есть $d = \text{НОД}(x_i - x_j, m)$ — нетривиальный делитель числа $m$.

Существует несколько схем поиска цикла в рекуррентной последовательности. Самая простая из них состоит в том, что мы последовательно сравниваем элементы $$ \begin{array}{lll} x_0 \leftrightarrow x_1,\quad & x_1 \leftrightarrow x_3,\quad & x_2 \leftrightarrow x_5,\\ x_3 \leftrightarrow x_7,\quad & x_4 \leftrightarrow x_9,\quad & x_5 \leftrightarrow x_{11},\ \dots \end{array} $$ То есть в потенциально бесконечном цикле сравниваются элементы $x$ и $у$. На каждой итерации к элементу $x$ применяется отображение $f$, к элементу $y$ отображение $f$ применяется дважды: $$ x = f(x), \quad y = f(f(y)). $$ Любая рекуррентная последовательность, состоящая из элементов конечного множества, имеет в своем начале отрезок апериодичности, после чего входит в цикл. За счет того, что элемент $x$ на каждой итерации цикла сдвигается на одину позицию вперед по последовательности, он рано или поздно войдет в цикл. Элемент $y$ на каждой итерации сдвигается на 2 позиции, значит, войдя в цикл, он рано или поздно догонит $x$ и мы придем к ситуации, когда $x = y$, т.е. $x_i = x_j$ при $i < j$.

В $\rho$-алгоритме Полларда вместо сравнения элементов $x$ и $y$ мы вычисляем $d = \text{НОД}(x - y, m)$. Если мы вошли в цикл в последовательности $\{\overline{x}_i\}$, но еще не вошли в цикл в последовательности $\{x_i\}$, то $$ x \equiv y \pmod p, \quad x \not\equiv y \pmod m. $$ Это означает, что $p | d = \text{НОД}(x - y, m)$, $m \nmid d$, т.е. $d$ — нетривиальный делитель числа $m$.

Запишем $\rho$-алгоритм факторизации на псевдокоде. Алгоритм rhoFactorization(m, maxSteps) возвращает собственный делитель числа $m$ или 1, если разложить число на множители не удается. Чтобы алгоритм не работал неопределенно долго, количество итераций ограничено числом maxSteps, по умолчанию maxSteps=1000000.

алгоритм rhoFactorization(m, maxSteps = 1000000)
| x = случайное целое число в диапазоне 2 <= x < m-1
| y = x*x + 1 (mod m)
| d = НОД(x - y, m)
| steps = 0
|
| цикл пока d == 1 и d != m и steps < maxSteps:
| | x = x*x + 1 (mod m)
| | y = y*y + 1 (mod m)
| | y = y*y + 1 (mod m)
| | d = НОД(x - y, m)
| | steps = steps + 1
| конец цикла
|
| если d == m   # Неудача...
| | то d = 1
| конец если
|
| вернуть d
конец алгоритма

Реализация алгоритма на языке Python3 предлагается читателю в качестве упражнения.

Приведем пример использования этого алгоритма. Сначала мы генерируем два случайных 10-значных простых числа $p$ и $q$, используя вероятностный тест простоты, и вычисляем $m = p\cdot q$:

>>> p = random.randrange(10**9, 10**10)
>>> while not primeTest(p):
...     p -= 1
... 
>>> p
6658968373
>>> q = random.randrange(10**9, 10**10)
>>> while not primeTest(q):
...     q -= 1
... 
>>> q
4312161011
>>> m = p*q
>>> m
28714543791532705103

>>> rhoFactorization(m)
1

>>> rhoFactorization(m, maxSteps=10000000)
4312161011

Оценка времени работы ρ-алгоритма Полларда

Количество итераций в $\rho$-алгоритме Полларда пропорционально длине цикла в последовательности $\{\overline{x}_i, i=0, 1, 2,\dots\}$. Эта последовательность является орбитой элемента $\overline{x}_0\in \mathbb{Z}_p$ при действии отображения $$ f: \mathbb{Z}_p \to \mathbb{Z}_p,\quad f(x) \equiv x^2 + 1 \pmod{p}. $$ Известно, что при полиномиальном отображении $f: \mathbb{Z}_n\to \mathbb{Z}_n$, где $f$ — полином степени выше 1, средняя длина орбиты произвольного элемента пропорциональна $\sqrt{n}$. Таким образом, время работы $\rho$-алгоритма Полларда при разложении числа $m$ на множители в среднем пропорционально $\sqrt{p}$, где $p$ — минимальный простой делитель числа $m$. В худшем случае, когда $m$ есть произведение двух простых чисел, почти одинаковых по величине, время работы оценивается как $O(\sqrt[4]{m})$.

P-1-алгоритм факторизации основан на Малой теореме Ферма. Пусть нужно разложить на множители число $m$ и пусть $p | m$ — простой делитель числа $m$. Предположим, что число $p-1$ гладкое (smooth) — число называется $N$-гладким. если оно раскладывается в произведение не очень больших степей простых чисел: $$ p - 1 = q_1^{\alpha_1} q_2^{\alpha_2} \dots q_k^{\alpha_k}, \text{ причем } q_i^{\alpha_i} \le N \text{ для всех }i=1,\dots, k. $$

Пусть $S$ — множество всех максимальных степеней простых чисел, не превосходящих $N$: $$ \begin{array}{rcl} S & = & \{ t_i^{\beta_i},\ % i=1, 2,\dots, k,\text{ где }t_i \le N \text{ —} \\ & & \text{всевозможные простые числа и } t_i^{\beta_i} \le N,\ t_i^{\beta_i + 1} \gt N\}. \end{array} $$ Обозначим элементы множества $S$ через $s_i$: $$ S = \{s1, s2,\dots, s_k\}. $$ Например, при $N = 20$ $$ S = \{ 16, 9, 5, 7, 11, 13, 17, 19 \}. $$ Верхняя граница $N$ является параметром алгоритма: чем больше $N$, тем большие числа может раскладывать алгоритм, но и тем дольше он работает.

Выберем случайное основание $b_0\in\mathbb{Z}$, $2 \le b \lt m-1$, и рассмотрим последовательность $$ b_0,\ b1 \equiv b_0^{s_1},\ b_2 \equiv b_1^{s_2},\ \dots,\ % b_k \equiv b_{s_{k-1}}^{s_k} \pmod{m}. $$ Каждый следующий элемент в этой последовательности получается из предыдущего возведением в степень $s_i$ в кольце $\mathbb{Z}_m$: $$ b_i \equiv b_{i-1}^{s_i} \pmod m. $$ Последний элемент в этой последовательности равен $$ b_k \equiv b_0^{s_1 s_2\dots s_k} \pmod{m}. $$

Из предположения гладкости числа $p - 1$ вытекает, что $$ (p - 1) | s_1 s_2\dots s_k, $$ то есть $$ s_1 s_2\dots s_k = (p - 1)r\ \text{для некоторого } r\in\mathbb{Z},\ r\ge 1. $$ Из Малой теоремы Ферма следует, что $$ b_k \equiv b_0^{s_1 s_2\dots s_k} \equiv b_0^{(p - 1)r} \equiv (b_0^{p - 1})^r \equiv 1 \pmod{p}. $$ Значит, $$ 0 \equiv b_k - 1 \pmod p,\ \text{т.е. } p \mid (b_k - 1). $$ Поскольку также $p\mid m$, то $$ p \mid \text{НОД}(b_k - 1, m) $$ и $d = \text{НОД}(b_k - 1, m)$ является нетривиальным делителем числа $m$.

Поскольку мы заранее не знаем степени гладкости простого делителя числа $m$, то $\text{НОД}(b_i - 1, m)$ вычисляется на каждом шаге $i = 0, 1, 2,\dots, k$. Алгоритм заканчивается, как только вычисленный наибольший общий делитель станет отличен от 1 (если при этом он равен $m$, то это означает неудачу).

Проиллюстрируем алгоритм на простом примере. Разложим на множители число $m = 102691 = 103\cdot 997$. Возьмем $N = 20$ и выпишем все степени простых чисел, не превосходящие 20: $$ S = \{ 16, 9, 5, 7, 11, 13, 17, 19 \}. $$ Выберем случайное число $b_0 = 2$. Последовательно вычисляем $$ \begin{array}{ll} & \text{НОД}(2 - 1, m) = 1, \\ 2^{16} \equiv 65536 \pmod{m}, & \text{НОД}(65536 - 1, m) = 1, \\ 65536^9 \equiv 62040 \pmod{m}, & \text{НОД}(62040 - 1, m) = 1, \\ 62040^5 \equiv 86069 \pmod{m}, & \text{НОД}(86069 - 1, m) = 1, \\ 86069^7 \equiv 72031 \pmod{m}, & \text{НОД}(72031 - 1, m) = 1, \\ 72031^{11} \equiv 46049 \pmod{m}, & \text{НОД}(46049 - 1, m) = 1, \\ 46049^{13} \equiv 101691\pmod{m}, & \text{НОД}(101691 - 1, m) = 1,\\ 101691^{17} \equiv 45115\pmod{m}, & \text{НОД}(45115 - 1, m) = 103. \end{array} $$ Алгоритм сработал на 7-м шаге, поскольку для делителя $p = 103$ числа $m$ число $p - 1$ гладкое (для $N = 20$): $$ p - 1 = 102 = 2 \cdot 3 \cdot 17, \quad 2, 3, 17 \le N. $$

Запишем алгоритм на псевдокоде. Параметром алгоритма является верхняя граница $N$, задающая меру гладкости числа $p - 1$, где $p$ — простой делитель числа $m$, которое раскладывается на множители. Назовем этот параметр алгоритма upperBound, по умолчанию upperBound=100000.

алгоритм p1Factorization(m, upperBound = 100000)
| b = случайное число в интервале 2 <= b < m - 1
| d = НОД(b - 1, m)
| p = 2         # простое число
|
| цикл пока d == 1 и p <= upperBound
| | # вычислим макс. степень p <= upperBound
| | s = p; s1 = s*p
| | цикл пока s1 <= upperBound
| | | s = s1; s1 = s1*p
| | конец цикла
| | b = b^s (mod m)
| | d = НОД(b - 1, m)
| |
| | # найдем следующее простое число
| | если p == 2
| | | p = 3
| | иначе
| | | p = p + 2
| | | цикл пока не primeTest(p)
| | | | p = p + 2
| | | конец цикла
| | конец если
| |
| конец цикла
|
| если d == m   # Неудача...
| | d = 1
| конец если
| вернуть d
конец алгоритма

Здесь мы использовали ранее реализованный алгоритм primeTest(p), определяющий простоту числа $p$ методом Миллера-Рабина. Реализация алгоритма на языке Python3 оставлена в качестве упражнения.

Приведем пример использования алгоритма. Сгенерируем сначала два 12-значных простых числа $p$ и $q$ и вычислим их произведение $m = pq$:

>>> p = random.randrange(10**11, 10**12)
>>> while not primeTest(p):
...     p -= 1
... 
>>> p
943187169361
>>> q = random.randrange(10**11, 10**12)
>>> while not primeTest(q):
...     q -= 1
... 
>>> q
585260672951
>>> m = p*q
>>> m
552010357458967668654311

>>> p1Factorization(552010357458967668654311)
585260672951

Алгоритм сработал, потому что $$ q - 1 = 585260672951 - 1 = 2\cdot 5^2\cdot 227\cdot 1129\cdot 45673, $$ т.е. число $q - 1$ гладкое для $N = 100000$ (значение параметра upperBound по умолчанию).

При реализации p-1 алгоритма Полларда p1Factorization мы воспользовались алгоритмом Миллера-Рабина primeTest(p), проверяющим простоту числа $p$. Можно, однако, обойтись и без этого теста. В упрощенном алгоритме p1Factorization2(m, upperBound) вместо произведения всех простых степеней, не превосходящих $N$, мы используем просто число $N! = 2\cdot 3\cdot 4\cdots N$. Количество итераций цикла при этом увеличивается незначительно, потому что простых чисел очень много $\pi(N) \sim N/\ln N$, так что число итераций возрастает не более чем в $\ln N$ раз. Зато не тратится время на проверку простоты числа $p$. Вот реализация упрощенного варианта p-1 алгоритма Полларда на псевдокоде. Реализация на Python'е оставлена как упражнение читателю.

алгоритм p1Factorization2(m, upperBound = 100000)
| b = случайное число в интервале 2 <= b < m - 1
| d = НОД(b - 1, m)
| s = 2
|
| цикл пока d == 1 и s <= upperBound
| | b = b^s (mod m)
| | d = НОД(b - 1, m)
| | s = s + 1
| конец цикла
|
| если d == m   # Неудача...
| | d = 1
| конец если
| вернуть d
конец алгоритма

Пример использования упрощенного алгоритма: генерируем два 12-значных простых числа $p$, $q$ и вычисляем $m = pq$:

>>> p = random.randrange(10**11, 10**12)
>>> while not primeTest(p):
...     p -= 1
... 
>>> p
286850613971
>>> q = random.randrange(10**11, 10**12)
>>> while not primeTest(q):
...     q -= 1
... 
>>> q
991236419323
>>> m = p*q
>>> m
284336775473218158161633

>>> p1Factorization2(284336775473218158161633)
1

>>> p1Factorization2(284336775473218158161633, upperBound=1000000)
991236419323