Теория чисел и математика для алгоритмов и криптографии

Теория чисел долго была витриной «бесполезной» математики. Годфри Харди в 1940 году в «Апологии математика» с гордостью писал, что теория чисел не имеет и не будет иметь никаких военных или практических применений. Через тридцать шесть лет вышла статья Диффи и Хеллмана, ещё через год — RSA, и оказалось, что вся электронная коммерция мира держится ровно на том, что Харди считал чистым искусством: на простых числах, на арифметике остатков и на асимметрии между «перемножить два простых» и «разложить произведение обратно».

Сегодня этот раздел нужен инженеру по трём совершенно разным причинам, и их полезно не путать.

1. Как источник примитивов. НОД, быстрое возведение в степень, решето, обратный элемент по модулю — это такие же базовые кирпичи, как бинарный поиск. Они всплывают в задачах, которые внешне не про числа: сокращение дробей, циклические сдвиги, периодичность, раскладка данных по шардам, генерация детерминированных перестановок.
2. Как арифметика конечных полей. Полиномиальное хеширование строк, NTT-умножение многочленов, коды Рида — Соломона, zk-SNARK, гомоморфное шифрование — всё это работает в кольце Z/m или в поле GF(p). Не понимая, почему модуль должен быть простым и что такое первообразный корень, вы будете копировать эти алгоритмы, не имея возможности их отладить.
3. Как основа криптографии. RSA, Диффи — Хеллман, DSA, эллиптические кривые. Здесь важна не только корректность, но и то, чего в обычных алгоритмах нет вовсе: устойчивость к атакам по побочным каналам, качество генератора случайных чисел, поведение при сбое оборудования.

Мы предполагаем знакомство с языком асимптотики и инвариантов (Анализ алгоритмов) и с идеей разделяй-и-властвуй (Рекурсия и разделяй-и-властвуй) — на неё мы обопрёмся, когда дойдём до быстрого умножения.

1. Делимость и алгоритм Евклида

1.1 Основная теорема арифметики

Всякое целое n > 1 раскладывается в произведение простых единственным способом с точностью до порядка множителей. Это утверждение выглядит очевидным ровно до тех пор, пока вы не увидите кольцо, где оно ложно: в Z[√−5] число 6 раскладывается двумя разными способами — 2 · 3 и (1 + √−5)(1 − √−5). Единственность разложения в Z — не тавтология, а теорема, и доказывается она через лемму Евклида, которая, в свою очередь, — через соотношение Безу. Поэтому начинать надо с НОД.

1.2 Евклид

Алгоритм Евклида — вероятно, самый старый нетривиальный алгоритм в истории (Начала, книга VII, около 300 года до н.э.). Он опирается на одно наблюдение: если d делит a и b, то d делит и a − b, а значит и a mod b. Множество общих делителей пары (a, b) совпадает с множеством общих делителей пары (b, a mod b) — а вторая пара «меньше».

def gcd(a, b):
    """НОД. Инвариант: множество общих делителей (a, b) не меняется."""
    while b:
        a, b = b, a % b
    return a          # b == 0 → НОД равен a

Сложность. Ключевая лемма: a mod b < a / 2, когда b ≤ a. Действительно, если b ≤ a/2, то a mod b < b ≤ a/2; если b > a/2, то a mod b = a − b < a/2. Значит за две итерации первый аргумент уменьшается минимум вдвое, и итераций не больше 2·log₂ a. Точная граница — теорема Ламе (1844): худший случай достигается на соседних числах Фибоначчи, число делений не превосходит log_φ(√5 · a) ≈ 4.785 · log₁₀ a. Для чисел до 10¹⁸ это ~90 итераций.

Для машинных слов операции — O(1), итого O(log min(a, b)). Для длинных чисел n бит наивный Евклид даёт O(n²); на практике GMP использует полу-рекурсивный алгоритм Лемера — Кнута со сложностью O(M(n) log n), где M(n) — стоимость умножения.

Практическая деталь: бинарный НОД (алгоритм Стайна) заменяет деление на сдвиги и вычитания. На процессорах, где деление дорого, он быстрее, но на современных x86-64 с быстрым div и на маленьких числах разница мала. В GMP используется гибрид.

1.3 Расширенный Евклид и соотношение Безу

Теорема Безу. Для любых целых a, b существуют x, y такие, что ax + by = gcd(a, b). Причём gcd(a, b) — наименьшее положительное число вида ax + by.

Расширенный алгоритм ведёт коэффициенты параллельно с остатками:

def ext_gcd(a, b):
    """Возвращает (g, x, y) такие, что a*x + b*y == g == gcd(a, b).

    Инвариант на каждой итерации:
        a*old_s + b*old_t == old_r
        a*s     + b*t     == r
    Обе строки — линейные комбинации исходных a и b, поэтому
    вычитание одной из другой сохраняет свойство.
    """
    old_r, r = a, b
    old_s, s = 1, 0
    old_t, t = 0, 1
    while r:
        q = old_r // r
        old_r, r = r, old_r - q * r
        old_s, s = s, old_s - q * s
        old_t, t = t, old_t - q * t
    return old_r, old_s, old_t     # old_r == gcd, коэффициенты — old_s, old_t

Сложность та же, O(log min(a, b)) операций. Коэффициенты не взрываются: |x| ≤ b / (2g), |y| ≤ a / (2g), так что переполнения при 64-битных входах не будет, если считать в 64 битах со знаком (но осторожно — промежуточные q*s требуют внимания при a, b под 2⁶³).

Из Безу немедленно следуют два факта, на которых стоит вся дальнейшая статья: лемма Евклида (если простое p делит ab, то p делит a или b) и критерий обратимости: a обратим по модулю m тогда и только тогда, когда gcd(a, m) = 1.

2. Модульная арифметика

2.1 Кольцо вычетов

Z/m — множество остатков {0, 1, …, m−1} со сложением и умножением по модулю m. Сложение, вычитание и умножение работают «как обычно»: (a + b) mod m, (a · b) mod m. Деление — не работает. Точнее, оно определено не для всех элементов.

Обратимые элементы Z/m (те, для которых существует a⁻¹ с a · a⁻¹ ≡ 1) — это ровно числа, взаимно простые с m. Их количество обозначают φ(m) — функция Эйлера. Они образуют мультипликативную группу (Z/m)* порядка φ(m).

Когда m = p простое, φ(p) = p − 1, обратим каждый ненулевой элемент, и Z/p — поле. Именно поэтому в соревновательном программировании модули берут простыми (10⁹+7, 998244353), а в хешировании строк простой модуль спасает от систематических коллизий.

2.2 Обратный элемент — три способа

def inv_ext(a, m):
    """Обратный через расширенный Евклид. Работает для ЛЮБОГО m при gcd(a, m) == 1.
    O(log m)."""
    g, x, _ = ext_gcd(a % m, m)
    if g != 1:
        raise ValueError(f"обратного нет: gcd({a}, {m}) = {g}")
    return x % m                      # приводим в [0, m)

def inv_fermat(a, p):
    """Только для ПРОСТОГО p. Малая теорема Ферма: a^(p-1) ≡ 1, значит a^(p-2) ≡ a^(-1).
    O(log p) умножений — чуть медленнее Евклида, но короче в записи."""
    return pow(a, p - 2, p)

def inv_batch(n, p):
    """Обратные ко всем 1..n за O(n) — трюк с рекуррентой.
    Из p = (p // i) * i + p % i следует i^(-1) ≡ -(p // i) * (p % i)^(-1)."""
    r = [0] * (n + 1)
    r[1] = 1
    for i in range(2, n + 1):
        r[i] = (p - (p // i) * r[p % i] % p) % p
    return r

В Python 3.8+ есть встроенный вариант: pow(a, -1, m) — он реализован на C через расширенный Евклид и работает для любого модуля. В Go — new(big.Int).ModInverse(a, m), в Java — BigInteger.modInverse.

2.3 Быстрое возведение в степень

def binpow(a, e, m):
    """a^e mod m за O(log e) умножений. Инвариант: result * a^e == исходное a^e."""
    result = 1
    a %= m
    while e:
        if e & 1:
            result = result * a % m
        a = a * a % m
        e >>= 1
    return result

Это тот же приём «разделяй-и-властвуй», что и в быстром умножении матриц: a^e = (a^(e/2))² при чётном e. Его применяют не только к числам — к матрицам (линейные рекурренты, числа Фибоначчи за O(log n)), к перестановкам, к любой ассоциативной операции (Рекурсия и разделяй-и-властвуй).

Осторожно в криптографии. Наивная реализация ветвится по битам секретной экспоненты. Пол Кохер в 1996 году показал, что по времени работы и по энергопотреблению эти биты восстанавливаются (Timing Attacks on Implementations of Diffie-Hellman, RSA, DSS). Криптобиблиотеки используют «лестницу Монтгомери» — вариант, выполняющий одинаковый набор операций независимо от бита, — плюс ослепление экспоненты.

2.4 Умножение по модулю и переполнение

В Python целые длинные, поэтому a * b % m всегда корректно. В C++/Java/Go при m около 10¹⁸ произведение не помещается в 64 бита. Варианты:

// 1) __int128 — просто и быстро на x86-64 и ARM64 (gcc/clang)
uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
    return (__uint128_t)a * b % m;
}

// 2) Барретт: заменяем деление на умножение на предвычисленную "обратную" константу.
//    Компилятор делает это сам, если модуль — константа времени компиляции.
// 3) Монтгомери: держим числа в специальном представлении aR mod m,
//    тогда редукция — только умножения и сдвиги, без деления вовсе.

Редукция Монтгомери (Montgomery, 1985) — основной приём во всех современных криптобиблиотеках и в zk-стеках: элементы поля хранятся в форме a·R mod m (где R = 2⁶⁴ᵏ), и умножение сводится к трём умножениям без деления. Переход в форму и обратно стоит по одной редукции, поэтому выгода появляется, когда умножений много подряд — что и происходит в возведении в степень.

Ещё одна ловушка: знак остатка. В C, C++, Java, Go, Rust и C# результат % для отрицательного левого операнда отрицателен (-7 % 3 == -1), а в Python — нет (-7 % 3 == 2). Каноническая нормализация: ((a % m) + m) % m.

3. Простые числа: решёта

3.1 Классическое решето

def sieve(n):
    """Все простые до n. Время O(n log log n), память O(n) байт."""
    is_composite = bytearray(n + 1)
    primes = []
    for i in range(2, n + 1):
        if not is_composite[i]:
            primes.append(i)
            # старт с i*i: всё, что меньше, уже вычеркнуто меньшими простыми
            for j in range(i * i, n + 1, i):
                is_composite[j] = 1
    return primes

Откуда log log n. Внутренний цикл для простого p делает n/p шагов, суммарно n · Σ_{p ≤ n} 1/p. По второй теореме Мертенса Σ_{p ≤ n} 1/p = ln ln n + M + o(1), где M ≈ 0.2615 — константа Мейсселя — Мертенса. Отсюда Θ(n log log n). Практически log log n для n = 10⁹ равно примерно 3 — то есть решето почти линейно.

Плотность простых даёт теорема о распределении простых чисел: π(n) ~ n / ln n. Для n = 10⁹ простых около 50.8 миллионов. Отсюда же полезная оценка «k-е простое ≈ k ln k» — она нужна, когда надо заранее выделить массив.

3.2 Линейное решето и таблица наименьших делителей

Классическое решето вычёркивает число столько раз, сколько у него разных простых делителей. Линейное решето вычёркивает каждое ровно один раз — его наименьшим простым делителем, что даёт бесплатный побочный продукт: таблицу spf (smallest prime factor).

def linear_sieve(n):
    """Простые + таблица наименьших простых делителей за O(n) времени и памяти."""
    spf = [0] * (n + 1)
    primes = []
    for i in range(2, n + 1):
        if spf[i] == 0:                 # i не вычеркнуто — значит простое
            spf[i] = i
            primes.append(i)
        for p in primes:
            # ключевое условие: p не больше наименьшего делителя i,
            # тогда p — наименьший делитель числа i*p, и мы попадём в него ровно раз
            if p > spf[i] or i * p > n:
                break
            spf[i * p] = p
    return primes, spf

def factorize(x, spf):
    """Факторизация за O(log x) при готовой таблице spf."""
    result = {}
    while x > 1:
        p = spf[x]
        while x % p == 0:
            result[p] = result.get(p, 0) + 1
            x //= p
    return result

Тот же каркас позволяет за O(n) посчитать любую мультипликативную функцию: φ, μ, число делителей d(n), сумму делителей σ(n). Это стандартный приём в задачах, где нужны все значения функции до n.

Trade-off. Линейное решето асимптотически лучше, но на практике для n ≤ 10⁸ классическое часто быстрее: оно работает с bytearray последовательно и дружелюбно к кэшу, а линейное делает случайные записи в spf и хранит int-массив (в 4–8 раз больше памяти). Если нужны только простые — берите классическое с оптимизацией «только нечётные» и битовой упаковкой. Если нужна факторизация многих чисел — линейное. Это ровно тот случай, о котором говорит статья Практическая оптимизация: асимптотика не единственный критерий.

3.3 Сегментированное решето

Когда нужны простые в диапазоне [L, R] при R ~ 10¹², массив на R не влезет в память. Решение — просеивать окно: базовых простых до √R хватает, чтобы вычеркнуть всё составное.

import math

def segmented_sieve(lo, hi):
    """Простые в [lo, hi]. Память O(sqrt(hi) + (hi - lo)), время O((hi-lo) log log hi)."""
    base = sieve(math.isqrt(hi))
    mark = bytearray(hi - lo + 1)
    for p in base:
        # первое кратное p, которое >= lo и >= p*p
        start = max(p * p, (lo + p - 1) // p * p)
        for j in range(start, hi + 1, p):
            mark[j - lo] = 1
    return [lo + i for i, m in enumerate(mark) if not m and lo + i >= 2]

Сегментированное решето — это ещё и способ уложиться в кэш: если брать окна по 32 КБ, все вычёркивания идут внутри L1, и решето ускоряется в разы. На этой идее построен primesieve — самая быстрая открытая реализация, считающая π(10¹⁵) за минуты.

4. Функция Эйлера, порядки и первообразные корни

Функция Эйлера φ(n) — количество чисел из [1, n], взаимно простых с n. Она мультипликативна (φ(ab) = φ(a)φ(b) при gcd(a,b)=1) и вычисляется по разложению: φ(n) = n · Π_{p | n} (1 − 1/p).

def phi(n):
    """Функция Эйлера через пробное деление. O(sqrt n)."""
    result = n
    p = 2
    while p * p <= n:
        if n % p == 0:
            while n % p == 0:
                n //= p
            result -= result // p     # умножаем на (1 - 1/p) без дробей
        p += 1
    if n > 1:                          # остался большой простой делитель
        result -= result // n
    return result

def phi_sieve(n):
    """Все значения φ(1..n) за O(n log log n) — решето по простым."""
    f = list(range(n + 1))
    for p in range(2, n + 1):
        if f[p] == p:                  # p простое
            for j in range(p, n + 1, p):
                f[j] -= f[j] // p
    return f

Теорема Эйлера. Если gcd(a, n) = 1, то a^φ(n) ≡ 1 (mod n). Частный случай при простом n = p — малая теорема Ферма: a^(p−1) ≡ 1 (mod p).

Следствие, которое постоянно нужно на практике: показатели можно приводить по модулю φ(n). a^e mod n = a^(e mod φ(n)) mod n — но только при gcd(a, n) = 1. Это одна из самых частых ошибок: при gcd(a, n) ≠ 1 формула ломается, и нужен обобщённый вариант (теорема о «поднятии показателя»): при e ≥ log₂ n верно a^e ≡ a^(φ(n) + e mod φ(n)) (mod n).

Порядок элемента ord_n(a) — минимальное k > 0 с a^k ≡ 1 (mod n). По теореме Лагранжа ord_n(a) делит φ(n). Элемент, у которого порядок равен ровно φ(n), называется первообразным корнем — он порождает всю группу.

def primitive_root(p, prime_factors_of_p_minus_1):
    """Первообразный корень по простому модулю p.
    Проверяем: g^((p-1)/q) != 1 для каждого простого делителя q числа p-1.
    Этого достаточно — порядок делит p-1 и не делит ни один максимальный делитель."""
    for g in range(2, p):
        if all(pow(g, (p - 1) // q, p) != 1 for q in prime_factors_of_p_minus_1):
            return g

Первообразные корни «плотны»: наименьший обычно очень мал (для 998244353 это 3), поэтому перебор сходится за единицы итераций, а дорогая часть — факторизация p − 1.

5. Китайская теорема об остатках

Формулировка. Если m₁, …, m_k попарно взаимно просты и M = Π m_i, то отображение x ↦ (x mod m₁, …, x mod m_k) — биекция между Z/M и Z/m₁ × … × Z/m_k, причём согласованная со сложением и умножением (изоморфизм колец).

Практический смысл: арифметика по большому модулю распадается на независимые операции по маленьким. Это и ускорение, и способ обойти переполнение.

def crt2(r1, m1, r2, m2):
    """Решает x ≡ r1 (mod m1), x ≡ r2 (mod m2) БЕЗ требования взаимной простоты.
    Возвращает (остаток, модуль) или None, если система несовместна. O(log min(m1, m2))."""
    g, p, _ = ext_gcd(m1, m2)
    if (r2 - r1) % g:
        return None                      # несовместно: r1 и r2 расходятся по модулю g
    lcm = m1 // g * m2                   # делим до умножения — бережём разрядность
    # ищем k из r1 + m1*k ≡ r2 (mod m2)
    k = (r2 - r1) // g * p % (m2 // g)
    return (r1 + m1 * k) % lcm, lcm

Классическая формула для попарно взаимно простых модулей: x = Σ r_i · M_i · (M_i⁻¹ mod m_i), где M_i = M / m_i. На SVG выше видно, откуда берутся «индикаторы» M_i · (M_i⁻¹ mod m_i): это числа, равные 1 по своему модулю и 0 по всем остальным — базис координат.

Где применяется:

• RSA-CRT. Расшифрование c^d mod n при n = pq считают отдельно mod p и mod q с экспонентами d mod (p−1), d mod (q−1), а потом собирают. Экспоненты и модули вдвое короче, а стоимость возведения в степень примерно кубическая по длине — итого ускорение почти в 4 раза. Так делают все библиотеки; формат приватного ключа PKCS#1 прямо хранит p, q, dP, dQ, qInv.
• Multi-modular арифметика. Свёртка больших чисел считается NTT по нескольким простым модулям вида k·2ⁿ + 1, результат собирается КТО. Так работают flint, NTL, и так же устроены реализации умножения полиномов в задачах на большие коэффициенты.
• Хеширование. Двойной полиномиальный хеш (два разных простых модуля) — это КТО: пара хешей эквивалентна одному хешу по модулю-произведению, и вероятность коллизии падает мультипликативно (см. Строковые алгоритмы).

Опасность RSA-CRT: атака Боне — ДеМилло — Липтона (Bellcore, 1997). Если во время вычисления одной из двух половин произойдёт аппаратный сбой, то gcd(m'^e − m, n) мгновенно выдаёт p. Одна ошибка в одном бите — и приватный ключ раскрыт. Поэтому реализации верифицируют результат (пересчитывают m^e mod n и сравнивают) — оригинальная работа.

6. Тесты простоты

Проверка «простое ли число» и «разложи число» — задачи разной сложности. Первая решается за полиномиальное время, вторая — нет (насколько мы знаем). На этом разрыве стоит вся криптография с открытым ключом.

6.1 Почему тест Ферма недостаточен

Малая теорема Ферма даёт необходимое условие: если n простое, то a^(n−1) ≡ 1 (mod n) для любого a, не делящегося на n. Обратное неверно. Составные числа, проходящие тест по какой-то базе, называются псевдопростыми Ферма; хуже того, существуют числа Кармайкла, проходящие тест по всем базам, взаимно простым с n. Наименьшее — 561 = 3 · 11 · 17. Их бесконечно много (Alford, Granville, Pomerance, 1994). Поэтому тест Ферма непригоден как самостоятельный критерий.

6.2 Миллер — Рабин

Тест Миллера — Рабина усиливает Ферма ещё одним наблюдением: в поле Z/p уравнение x² ≡ 1 имеет ровно два решения, x = ±1. Значит, если по дороге возведения a^(n−1) в квадрат мы встретим нетривиальный квадратный корень из единицы, число составное.

Записываем n − 1 = d · 2^s с нечётным d и смотрим на последовательность a^d, a^(2d), a^(4d), …, a^(2^s·d) = a^(n−1). Для простого n она обязана либо начинаться с 1, либо содержать n − 1.

def is_prime(n):
    """Детерминированный тест простоты для всех n < 3.3 * 10^24 (покрывает весь uint64).
    Время: O(k log^3 n) битовых операций, k = 12 баз."""
    if n < 2:
        return False
    small = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)
    for p in small:
        if n % p == 0:
            return n == p

    d, s = n - 1, 0
    while d % 2 == 0:
        d //= 2
        s += 1

    for a in small:                      # эти 12 баз доказанно достаточны для n < 3.3e24
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue                     # свидетель не найден, база бесполезна
        for _ in range(s - 1):
            x = x * x % n
            if x == n - 1:
                break
        else:
            return False                 # ни разу не встретили -1 → составное, гарантированно
    return True

Свойства. Для случайной базы вероятность, что составное число пройдёт тест, не больше 1/4 (Рабин), на практике сильно меньше. При k независимых базах ошибка ≤ 4^(−k); k = 40 даёт 2^(−80) — надёжнее, чем вероятность аппаратного сбоя. При этом тест односторонний: «составное» — всегда правда, «простое» — вероятностное утверждение (если только базы не фиксированы для проверенного диапазона, как выше).

Проверенный набор из 12 первых простых покрывает n < 3 317 044 064 679 887 385 961 981 (результат Sorenson & Webster, 2015) — то есть весь 64-битный диапазон и заметно больше. Для 32-битных чисел хватает трёх баз: 2, 7, 61.

BPSW. Промышленный стандарт — тест Бейли — Померанца — Селфриджа — Вагстаффа: одна база Миллера — Рабина (2) плюс сильный тест Люка. Контрпримеров не известно ни одного, при том что искали целенаправленно и перебрали все числа до 2⁶⁴. Именно BPSW реализован в Go math/big.Int.ProbablyPrime (после n раундов Миллера — Рабина всегда добавляется Люка) — документация, и в GMP mpz_probab_prime_p.

AKS. В 2002 году Агравал, Каял и Саксена показали, что простота проверяется детерминированно за полиномиальное время без всяких гипотез (PRIMES is in P). Это огромный теоретический результат и абсолютно бесполезный практически: даже улучшенные варианты Õ(log⁶ n) проигрывают BPSW на порядки.

6.3 Генерация простых для криптографии

import secrets

def gen_prime(bits):
    """Простое заданной битности. Ожидаемое число попыток ~ ln(2^bits)/2 ≈ 0.35 * bits."""
    while True:
        n = secrets.randbits(bits) | (1 << (bits - 1)) | 1   # старший бит и нечётность
        if is_prime(n):
            return n

Здесь secrets, а не random: random использует Mersenne Twister, чьё состояние восстанавливается по 624 выходам. Реальная катастрофа на этой почве — Debian OpenSSL (2008), где патч сузил энтропию до 15 бит, и все сгенерированные за два года ключи оказались перебираемыми. Другая — исследование Mining Your Ps and Qs (USENIX 2012): собрав RSA-ключи со всего интернета, авторы посчитали попарные НОД и взломали 0.5% хостов — устройства с плохой энтропией порождали общие простые. gcd(n₁, n₂) ≠ 1 — и оба ключа мгновенно раскладываются.

7. Факторизация

7.1 Пробное деление и метод Ферма

Пробное деление до √n — O(√n), годится до 10¹²–10¹⁴. Метод Ферма ищет представление n = a² − b² = (a−b)(a+b), перебирая a от ⌈√n⌉: он быстр, когда делители близки друг к другу, и катастрофически медленный иначе. Отсюда практическое правило генерации RSA-ключей: |p − q| должно быть велико — иначе метод Ферма разложит модуль за секунды. FIPS 186-5 прямо требует |p − q| > 2^(nlen/2 − 100).

7.2 ρ-метод Полларда

Идея красивая и вероятностная. Возьмём псевдослучайную функцию f(x) = x² + c mod n и итерируем её. Последовательность рано или поздно зациклится. Если p — делитель n, то та же последовательность по модулю p зациклится гораздо раньше: по парадоксу дней рождения — через O(√p) = O(n^(1/4)) шагов. Момент, когда xᵢ ≡ xⱼ (mod p), но xᵢ ≠ xⱼ (mod n), детектируется как gcd(|xᵢ − xⱼ|, n) ∈ (1, n).

import math, random

def pollard_rho(n):
    """Находит НЕтривиальный делитель составного n. Ожидаемое время O(n^(1/4))."""
    if n % 2 == 0:
        return 2
    while True:
        x = random.randrange(2, n)
        c = random.randrange(1, n)
        y, d = x, 1
        while d == 1:
            x = (x * x + c) % n              # черепаха — один шаг
            y = (y * y + c) % n
            y = (y * y + c) % n              # заяц — два шага (цикл Флойда)
            d = math.gcd(abs(x - y), n)
        if d != n:
            return d                          # неудача d == n → меняем c и пробуем снова

def factor(n):
    """Полная факторизация. Сначала проверяем простоту — иначе рекурсия не остановится."""
    if n == 1:
        return []
    if is_prime(n):
        return [n]
    d = pollard_rho(n)
    return sorted(factor(d) + factor(n // d))

Улучшение Брента (1980): вместо цикла Флойда — поиск цикла удвоением, и главное — батчинг НОД: перемножаем 128 разностей по модулю n и берём gcd один раз. gcd дороже умножения на порядок, так что это ускоряет в разы. С Брентом и __int128 факторизация 64-битного числа занимает микросекунды.

Для чисел больше 10²⁵ нужны субэкспоненциальные методы: квадратичное решето (msieve) и общее решето числового поля GNFS (CADO-NFS) со сложностью exp(((64/9)^(1/3) + o(1)) · (ln n)^(1/3) (ln ln n)^(2/3)). Рекорд для RSA-модуля — RSA-250 (829 бит), разложен в 2020 году примерно за 2700 core-years. RSA-2048 остаётся недосягаемым для классических компьютеров; алгоритм Шора на достаточно большом квантовом компьютере сломает его за полиномиальное время — именно поэтому NIST в 2024 году стандартизировал постквантовые схемы ML-KEM и ML-DSA.

8. Дискретный логарифм

Задача: найти x из g^x ≡ b (mod p). В Z/p она субэкспоненциальна (index calculus), в группе точек эллиптической кривой — только O(√n), и именно поэтому ECC даёт ту же стойкость при вчетверо более коротких ключах: 256-битная кривая ≈ 3072-битный RSA.

Baby-step giant-step (Шенкс) — классический пример размена памяти на время в духе метода встречи посередине. Пишем x = i·n − j, где n = ⌈√p⌉, 0 ≤ j < n. Тогда g^(i·n) = b · g^j: слева n «больших шагов», справа n «малых».

def bsgs(g, b, p):
    """Решает g^x ≡ b (mod p) для простого p и gcd(g, p) = 1.
    Время O(sqrt(p)) операций, память O(sqrt(p)) — типичный meet-in-the-middle."""
    g %= p
    b %= p
    n = math.isqrt(p) + 1

    table = {}                       # малые шаги: b * g^j
    cur = b
    for j in range(n + 1):
        table[cur] = j               # при коллизии остаётся больший j → минимальный x
        cur = cur * g % p

    gn = pow(g, n, p)                # большие шаги: (g^n)^i
    cur = 1
    for i in range(1, n + 1):
        cur = cur * gn % p
        if cur in table:
            x = i * n - table[cur]
            if x > 0:
                return x
    return None

ρ-метод Полларда для дискретного логарифма даёт то же O(√p) времени, но O(1) памяти, и распараллеливается методом «различимых точек» (van Oorschot — Wiener). Именно он ломал рекорды ECDLP — например, Certicom ECC2K-130.

9. Криптография на практике

9.1 Диффи — Хеллман

Ключевой момент: DH защищает только от пассивного наблюдателя. Активный человек-посередине подменяет A и B своими и договаривается с обоими по отдельности — поэтому в TLS обмен ключами всегда подписан долгоживущим ключом сервера, чей сертификат проверяется по цепочке доверия.

Практические уроки:

• Logjam (2015) — многие серверы использовали одну и ту же 512-битную группу. Предвычисление index calculus для одной группы стоит месяцы, но потом каждый отдельный лог считается за минуты (статья). Отсюда правило: не меньше 2048 бит и стандартные группы из RFC 3526 / RFC 7919.
• Проверяйте чужой публичный ключ. Если принять A = 0, 1, p−1, общий секрет становится предсказуемым (small subgroup confinement).
• Современный выбор — не Z/p, а X25519 на кривой Curve25519 (Bernstein, 2006): нет «плохих» точек, реализация естественно постоянна по времени.

9.2 RSA

def rsa_keygen(bits=2048):
    """Учебная генерация ключа. В проде — используйте библиотеку, не это."""
    p = gen_prime(bits // 2)
    q = gen_prime(bits // 2)
    n = p * q
    lam = (p - 1) * (q - 1) // math.gcd(p - 1, q - 1)   # функция Кармайкла λ(n)
    e = 65537                                            # мало единиц в двоичной записи → быстро
    assert math.gcd(e, lam) == 1
    d = pow(e, -1, lam)
    return (n, e), (n, d, p, q)

def rsa_crt_decrypt(c, n, d, p, q):
    """Расшифрование через КТО — примерно в 4 раза быстрее прямого pow(c, d, n)."""
    dp, dq = d % (p - 1), d % (q - 1)
    qinv = pow(q, -1, p)
    m1, m2 = pow(c, dp, p), pow(c, dq, q)
    h = qinv * (m1 - m2) % p
    return (m2 + h * q) % n

Корректность: (m^e)^d = m^(ed) = m^(1 + kλ(n)) = m по теореме Эйлера — Кармайкла.

«Учебный» RSA небезопасен, и это не педантизм. Он детерминирован (одинаковые сообщения дают одинаковые шифртексты), мультипликативен (E(m₁)·E(m₂) = E(m₁m₂) — отсюда атаки с выбранным шифртекстом), уязвим для малых сообщений при малом e (если m^e < n, корень берётся просто так). Реальные схемы — OAEP для шифрования и PSS для подписи, и даже с ними история знает атаку Блайхенбахера на PKCS#1 v1.5 (padding oracle, 1998), которая воскресала в 2017 году как ROBOT.

ROCA (2017). Библиотека Infineon генерировала простые вида p = k·M + (65537^a mod M) ради скорости. Такая структура позволила разложить ключи методом Копперсмита за часы — затронуты миллионы паспортов и TPM-чипов (статья). Мораль: любое сужение пространства ключей ради производительности — потенциальная катастрофа.

10. Теория чисел вне криптографии

10.1 Комбинаторика по простому модулю

MOD = 10**9 + 7
MAXN = 200_001

fact = [1] * MAXN
for i in range(1, MAXN):
    fact[i] = fact[i - 1] * i % MOD

# обратные факториалы за ОДНО возведение в степень: ifact[i-1] = ifact[i] * i
ifact = [1] * MAXN
ifact[MAXN - 1] = pow(fact[MAXN - 1], MOD - 2, MOD)
for i in range(MAXN - 1, 0, -1):
    ifact[i - 1] = ifact[i] * i % MOD

def C(n, k):
    """Биномиальный коэффициент по простому модулю за O(1) после O(n) предподсчёта."""
    if k < 0 or k > n:
        return 0
    return fact[n] * ifact[k] % MOD * ifact[n - k] % MOD

Если n больше модуля, факториал обнуляется, и нужна теорема Люка: C(n, k) mod p равен произведению C(nᵢ, kᵢ) mod p по цифрам n и k в системе счисления с основанием p.

def lucas(n, k, p, f, finv):
    """C(n, k) mod p для простого p и произвольно больших n. O(log_p n)."""
    result = 1
    while n or k:
        ni, ki = n % p, k % p
        if ki > ni:
            return 0                 # цифра k больше цифры n → коэффициент делится на p
        result = result * f[ni] % p * finv[ki] % p * finv[ni - ki] % p
        n //= p
        k //= p
    return result

Когда модуль составной, деления делать нельзя вовсе: раскладывайте модуль, считайте по каждой степени простого (обобщение Эндрю Гранвилла) и собирайте КТО.

10.2 NTT — свёртка без чисел с плавающей точкой

FFT над комплексными числами быстра, но накапливает погрешность: при больших коэффициентах результат перестаёт округляться к правильному целому. NTT — это тот же алгоритм Кули — Тьюки, но в поле Z/p, где роль комплексного корня из единицы играет элемент порядка 2^k. Отсюда требование к модулю: p = c · 2^k + 1. Классический выбор — 998244353 = 119 · 2²³ + 1 с первообразным корнем 3: он поддерживает свёртки длиной до 2²³.

NMOD, ROOT = 998244353, 3

def ntt(a, invert=False):
    """Числовое преобразование Фурье в Z/998244353. O(n log n), n — степень двойки."""
    n = len(a)
    j = 0
    for i in range(1, n):                 # перестановка бит-реверса
        bit = n >> 1
        while j & bit:
            j ^= bit
            bit >>= 1
        j |= bit
        if i < j:
            a[i], a[j] = a[j], a[i]

    length = 2
    while length <= n:
        w = pow(ROOT, (NMOD - 1) // length, NMOD)   # корень степени length из единицы
        if invert:
            w = pow(w, NMOD - 2, NMOD)
        for i in range(0, n, length):
            wn, half = 1, length // 2
            for k in range(i, i + half):
                u = a[k]
                v = a[k + half] * wn % NMOD
                a[k] = (u + v) % NMOD
                a[k + half] = (u - v) % NMOD
                wn = wn * w % NMOD
        length <<= 1

    if invert:
        ninv = pow(n, NMOD - 2, NMOD)
        for i in range(n):
            a[i] = a[i] * ninv % NMOD
    return a

NTT — не академическая экзотика. Она в ядре постквантовых стандартов: ML-KEM (Kyber) работает в кольце Z_3329[x]/(x²⁵⁶+1) и умножает многочлены именно через NTT. Она же — основной потребитель времени в схемах гомоморфного шифрования (Microsoft SEAL, OpenFHE) и в zk-SNARK, где нужны свёртки над полем BLS12-381.

10.3 Мелочи, которые часто выручают

• Периодичность. Задачи вида «что будет на 10¹⁸-м шаге» почти всегда сводятся к нахождению порядка элемента или к возведению матрицы перехода в степень.
• Циклические сдвиги и шардирование. Шардирование по hash(key) mod N при смене N перемешивает всё; консистентное хеширование решает это, но арифметика остатков объясняет, почему проблема возникает.
• Сокращение дробей, рациональная арифметика. gcd — единственный способ держать числители и знаменатели ограниченными.
• Хеш-таблицы. Простой размер таблицы разрушает регулярность входных ключей — подробнее в курсе Структуры данных.

11. Типичные ошибки

1. Отрицательный остаток. (a - b) % m в C++/Java/Go может быть отрицательным. Всегда нормализуйте: ((a - b) % m + m) % m. В Python эта ловушка отсутствует, но при портировании кода она возвращается.
2. Переполнение при умножении. a * b % m при m > 2³¹ в 64-битных типах молча даёт мусор. Нужен __int128, Монтгомери или переход к другому языку.
3. Обратный элемент по составному модулю через Ферма. pow(a, m-2, m) верно только для простого m. Для составного нужен расширенный Евклид — и он же честно скажет, что обратного нет.
4. Деление в комбинаторике по составному модулю. C(n, k) mod 10⁹ (не простое) через обратные факториалы даст неправильный ответ — факториалы необратимы.
5. Приведение показателя по φ(n) без проверки gcd(a, n) = 1. Классический источник тонких багов в задачах на «башни степеней».
6. Тест Ферма вместо Миллера — Рабина. Число 561 радостно пройдёт.
7. Забытый i*i в решете. Работает, но в разы медленнее — и это заметно на n = 10⁸.
8. random вместо secrets в криптографии. Или time() как seed. Это не «немного хуже», это полное отсутствие безопасности.
9. Ветвление по секретным данным. Любой if от бита ключа — потенциальный timing-канал.
10. Слишком узкий модуль в хешировании строк. При m ~ 10⁹ парадокс дней рождения даёт коллизию уже на ~10⁵ строках; против злонамеренного входа нужен случайный модуль и случайное основание.

12. Что использовать в проде

Правило номер один: не пишите собственную криптографию. Всё из раздела 9 приведено, чтобы вы понимали, что происходит внутри, а не чтобы вы это использовали.

Отдельно про Python: pow(a, e, m) реализован на C и для больших модулей использует скользящее окно; pow(a, -1, m) даёт обратный элемент; math.isqrt — точный целочисленный корень без ошибок округления, которые даёт int(math.sqrt(n)) уже около 10¹⁶. Три эти функции покрывают 90% практических нужд, и писать их руками стоит только в учебных целях либо когда нужен свой Монтгомери в горячем цикле.

13. Мини-итог

• Всё начинается с Евклида: НОД, соотношение Безу и, как следствие, критерий обратимости по модулю. O(log n) — и это фундамент.
• Z/m — кольцо; Z/p — поле. Разница определяет, можно ли делить, и потому диктует выбор модуля почти во всех прикладных задачах.
• Решето даёт простые за Θ(n log log n), линейное решето — ещё и таблицу наименьших делителей за Θ(n); сегментированное позволяет работать с диапазонами до 10¹⁵.
• Проверка простоты дёшева (Миллер — Рабин, детерминированно для 64 бит), факторизация — нет. Этот разрыв и есть криптография с открытым ключом.
• КТО превращает большой модуль в набор маленьких: RSA-CRT, multi-modular NTT, двойное хеширование — это всё она.
• Дискретный логарифм и факторизация — две «трудные» задачи, на которых стоят DH и RSA; обе ломаются алгоритмом Шора, отсюда переход на решёточную криптографию.
• В криптографии корректности мало. Нужны качественная энтропия, постоянное время выполнения и защита от сбоев — почти все известные взломы били именно туда, а не в математику.

14. Источники

• Кормен, Лейзерсон, Ривест, Штайн. Алгоритмы: построение и анализ, глава 31 «Теоретико-числовые алгоритмы» — каноническое изложение с доказательствами.
• Victor Shoup. A Computational Introduction to Number Theory and Algebra — бесплатная и, вероятно, лучшая книга по теме: https://shoup.net/ntb/
• Menezes, van Oorschot, Vanstone. Handbook of Applied Cryptography — свободно доступен целиком: https://cacr.uwaterloo.ca/hac/
• Crandall, Pomerance. Prime Numbers: A Computational Perspective — про решёта, тесты простоты и факторизацию максимально подробно.
• Agrawal, Kayal, Saxena. PRIMES is in P, 2002 — https://www.cse.iitk.ac.in/users/manindra/algebra/primality_v6.pdf
• Pollard. A Monte Carlo method for factorization, BIT 1975 — https://link.springer.com/article/10.1007/BF01933667
• Brent. An improved Monte Carlo factorization algorithm, 1980 — https://maths-people.anu.edu.au/~brent/pd/rpb051i.pdf
• Montgomery. Modular multiplication without trial division, 1985 — https://www.ams.org/journals/mcom/1985-44-170/S0025-5718-1985-0777282-X/
• Heninger et al. Mining Your Ps and Qs, USENIX Security 2012 — https://factorable.net/weakkeys12.extended.pdf
• Nemec et al. The Return of Coppersmith’s Attack (ROCA), CCS 2017 — https://crocs.fi.muni.cz/public/papers/rsa_ccs17
• Adrian et al. Imperfect Forward Secrecy (Logjam), CCS 2015 — https://weakdh.org/imperfect-forward-secrecy-ccs15.pdf
• NIST FIPS 186-5 (цифровые подписи, требования к генерации простых) — https://csrc.nist.gov/pubs/fips/186-5/final
• NIST FIPS 203 (ML-KEM, постквантовый обмен ключами) — https://csrc.nist.gov/pubs/fips/203/final
• cp-algorithms — компактные эталонные реализации всего из этой статьи — https://cp-algorithms.com/algebra/module-inverse.html

Что дальше

Обратите внимание, что в этой статье впервые появился алгоритм, который может ошибиться: Миллер — Рабин с случайными базами не гарантирует правильный ответ, он гарантирует лишь, что вероятность ошибки меньше 4^(-k). И ρ-метод Полларда не имеет детерминированной оценки времени — только математическое ожидание. Это не недостаток реализации, а осознанный инженерный размен: отказавшись от абсолютных гарантий, мы получили алгоритмы, которые на порядки быстрее любых детерминированных аналогов.

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

Рандомизированные алгоритмы и вероятностный анализ