Простота и разложение больших чисел

Перебор делителей до корня проверяет простоту за $O(\sqrt{n})$. Для $n$ до $10^{12}$ это миллион операций — прекрасно. Для $n$ до $10^{18}$ это миллиард на каждый запрос — уже нет.

Существуют алгоритмы, которые справляются с 64-битными числами практически мгновенно. Оба вероятностные по природе, но для нашего диапазона их удаётся сделать полностью детерминированными.

Тест Ферма и почему он не работает

Начнём с наивной идеи. Малая теорема Ферма говорит: если $p$ простое, то $a^{p-1} \equiv 1 \pmod p$ для любого $a$, не кратного $p$.

Значит, можно взять несколько случайных $a$ и проверить. Не выполнилось — число точно составное. Выполнилось — вроде бы простое.

Ловушка в том, что существуют числа Кармайкла — составные, для которых сравнение выполняется при всех взаимно простых $a$. Наименьшее из них $561 = 3 \cdot 11 \cdot 17$, дальше $1105$, $1729$. Их бесконечно много, и тест Ферма на них ошибается всегда, сколько баз ни бери.

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

Тест Миллера—Рабина усиливает проверку одним наблюдением: у единицы по простому модулю ровно два квадратных корня, $1$ и $-1$. Составной модуль обычно даёт лишние корни, и на этом попадается.

Представим $n - 1 = d \cdot 2^s$ с нечётным $d$. Для честного простого последовательность

$$a^d, \; a^{2d}, \; a^{4d}, \; \dots, \; a^{2^{s-1} d}$$

обязана либо начинаться с $1$, либо где-то содержать $-1$. Если ни того, ни другого — $n$ составное, и это доказано, а не предположено.

long long mulmod(long long a, long long b, long long m) {
    return (__int128)a * b % m;
}

long long powmod(long long a, long long e, long long m) {
    long long r = 1; a %= m;
    while (e) { if (e & 1) r = mulmod(r, a, m); a = mulmod(a, a, m); e >>= 1; }
    return r;
}

bool isPrime(long long n) {
    if (n < 2) return false;
    for (long long p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37})
        if (n % p == 0) return n == p;

    long long d = n - 1;
    int s = 0;
    while (d % 2 == 0) { d /= 2; s++; }

    for (long long a : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
        long long x = powmod(a, d, n);
        if (x == 1 || x == n - 1) continue;
        bool composite = true;
        for (int i = 1; i < s; i++) {
            x = mulmod(x, x, n);
            if (x == n - 1) { composite = false; break; }
        }
        if (composite) return false;
    }
    return true;
}

Про набор баз стоит сказать отдельно. Первые двенадцать простых — $2, 3, \dots, 37$ — доказанно дают верный ответ для всех $n < 3.3 \cdot 10^{24}$. То есть для всего диапазона long long этот тест не вероятностный, а точный. Случайные базы брать не нужно.

Меньшие наборы тоже известны: до $3.2 \cdot 10^{18}$ хватает баз $2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37$; до $3.2 \cdot 10^{9}$ — всего $2, 3, 5, 7$. А вот наивный набор $2, 3, 5, 7$ на бо́льших числах ошибается: наименьший контрпример — $3215031751$.

mulmod через __int128 обязателен. Произведение двух чисел около $10^{18}$ — это $10^{36}$, никакой unsigned long long этого не удержит.

Ро-алгоритм Полларда

Простоту проверили. Если число составное, его ещё нужно разложить — и перебор до корня опять слишком медленный.

Идея Полларда красива. Возьмём псевдослучайную последовательность $x_{k+1} = (x_k^2 + c) \bmod n$. Она обязана зациклиться, и если нарисовать её как путь, получится греческая буква ρ: хвост и петля.

Пусть $p$ — какой-то делитель $n$. Та же последовательность по модулю $p$ зацикливается раньше — примерно через $\sqrt{p}$ шагов по парадоксу дней рождения. Значит, найдётся пара $x_i \ne x_j$ по модулю $n$, но $x_i \equiv x_j$ по модулю $p$. Тогда $\gcd(|x_i - x_j|, n)$ — нетривиальный делитель.

Пару ищут алгоритмом «черепаха и заяц»: один указатель делает шаг, другой два.

long long pollard(long long n) {
    if (n % 2 == 0) return 2;
    for (long long c = 1;; c++) {
        auto f = [&](long long x) { return (mulmod(x, x, n) + c) % n; };
        long long x = 2, y = 2, d = 1;
        while (d == 1) {
            x = f(x);
            y = f(f(y));
            d = __gcd(llabs(x - y), n);
        }
        if (d != n) return d;
    }
}

void factor(long long n, map<long long, int> &out) {
    if (n == 1) return;
    if (isPrime(n)) { out[n]++; return; }
    long long d = pollard(n);
    factor(d, out);
    factor(n / d, out);
}

Внешний цикл по $c$ нужен, потому что при неудачном выборе константы алгоритм упирается в $d = n$ и ничего не находит. Тогда берут следующую $c$ и пробуют снова.

Ожидаемое время — $O(n^{1/4})$. Для $n$ около $10^{18}$ это порядка тридцати тысяч операций: разложение произведения двух девятизначных простых занимает миллисекунды.

В factor важен порядок: сначала проверка на простоту, потом Поллард. Запускать ро-алгоритм на простом числе бессмысленно — он будет крутиться до бесконечности, перебирая $c$.

Когда что применять

ситуация	инструмент
много чисел до $10^7$	линейное решето, $\mathrm{lp}$
одно число до $10^{12}$	перебор до корня
одно число до $10^{18}$, нужна простота	Миллер—Рабин
одно число до $10^{18}$, нужно разложение	Миллер—Рабин + Поллард
простые в отрезке $[L, R]$, $R$ до $10^{12}$	сегментное решето

Отдельно стоит сказать: на школьных олимпиадах Поллард встречается редко. Но знать про его существование полезно — иначе задача с ограничением $10^{18}$ выглядит нерешаемой, хотя решается в двадцать строк.