Мотивация и контекст

Проверка простоты числа — одна из самых фундаментальных задач в программировании. Наивный алгоритм пробных делений работает за $O(\sqrt{N})$

​) и непрактичен при $N \approx 10^{18}$. Решето Эратосфена ограничено малыми $N \le 10^8$. Для соревновательных задач требуется эффективная проверка простоты одиночных чисел до $10^{18}$.

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

Числа, "обманывающие" тест Ферма для всех оснований $\gcd(a, n) = 1$, называются числами Кармайкла. Их существование разрушает наивную идею "тест Ферма = тест простоты". В этой главе мы изучим природу этих чисел и построим надёжный детерминированный тест простоты.

Числа Кармайкла на практике:

• $561 = 3 \cdot 11 \cdot 17$ — наименьшее.
• $1105, 1729, 2465, \ldots$ — их бесконечно много (теорема Алфорда–Гранвилля–Померанса, 1994).
• До $10^{21}$ их $\approx 246683$ — редки, но существуют в любом диапазоне.

Что мы построим:

1. Критерий Корселта — характеристика чисел Кармайкла.
2. Тест Миллера–Рабина — надёжный вероятностный тест.
3. Детерминированный тест для $n < 3.3 \cdot 10^{24}$ через фиксированные основания.
4. Тест Бейли–Псевдо–Простых (BPSW) — практически надёжный.

---

Теория

Малая теорема Ферма и псевдопростые

Теорема (Малая теорема Ферма). Если $p$ — простое, $\gcd(a, p) = 1$, то $a^{p-1} \equiv 1 \pmod{p}$.

Доказательство. Числа $a, 2a, 3a, \ldots, (p-1)a$ образуют перестановку $\{1, 2, \ldots, p-1\}$ по модулю $p$ (т.к. $a$ обратимо). Перемножая: $(p-1)! \cdot a^{p-1} \equiv (p-1)! \pmod{p}$, откуда $a^{p-1} \equiv 1 \pmod{p}$. $\square$

Определение. Составное $n$ называется псевдопростым по основанию $a$ (base-$a$ Fermat pseudoprime), если $a^{n-1} \equiv 1 \pmod{n}$.

Пример. $n = 341 = 11 \cdot 31$. Проверим: $2^{340} \equiv 1 \pmod{341}$. Действительно, $2^{10} = 1024 = 3 \cdot 341 + 1 \equiv 1 \pmod{341}$, и $340 = 10 \cdot 34$, поэтому $2^{340} = (2^{10})^{34} \equiv 1^{34} = 1$. Число $341$ — псевдопростое по основанию 2.

Числа Кармайкла

Определение. Составное $n$ называется числом Кармайкла, если $a^{n-1} \equiv 1 \pmod{n}$ для всех $a$ с $\gcd(a, n) = 1$.

Иначе говоря, числа Кармайкла — "абсолютные" псевдопростые Ферма.

Теорема Корселта (1899). $n > 1$ является числом Кармайкла тогда и только тогда, когда:

1. $n$ squarefree (не делится на $p^2$ ни для какого простого $p$).
2. Для каждого простого $p | n$: $(p-1) | (n-1)$.

Доказательство. ($\Rightarrow$) Пусть $n$ — число Кармайкла.

Squarefree: Предположим $p^2 | n$. Возьмём $a \equiv p \pmod{p^2}$ с $\gcd(a, n) = 1$ (можно найти $a = p + p^2 k$ с $\gcd(a, n/p^2) = 1$). Тогда $a^{n-1} \equiv p^{n-1} \pmod{p^2}$. Если $n - 1 \ge 1$ и $p^2 | n$, то $p^{n-1} \not\equiv 1 \pmod{p^2}$ при $p | a$... точнее, возьмём $a$ такое, что $p | a$ и $\gcd(a, n) = 1$. Но $p | a$ и $p | n$ означает $\gcd(a, n) \ge p > 1$ — противоречие. Значит, для $a = p$ выполнено $\gcd(p, n) = p > 1$, что не допускается. Поэтому рассмотрим другой $a$...

Строгое доказательство: при $p^2 | n$, если $\gcd(a, n) = 1$, то $p \nmid a$. По CRT можно взять $a \equiv 1 + p \pmod{p^2}$ (если $\gcd(1+p, n/p) = 1$). Тогда $a^{n-1} \equiv (1+p)^{n-1} \pmod{p^2}$. Биномом: $(1+p)^{n-1} \equiv 1 + (n-1)p \pmod{p^2}$. Для этого быть $\equiv 1$ нужно $p | (n-1)$. Но $p^2 | n$ значит $p | n$, значит $p | n - 1$ iff $p | -1$ — невозможно. Значит $a^{n-1} \not\equiv 1 \pmod{n}$ — противоречие с определением. Значит $n$ squarefree.

(p-1) | (n-1): Пусть $p | n$. Так как $n$ squarefree, $p \| n$. Возьмём $a$ — примитивный корень по модулю $p$. По CRT $\exists a$ с $a \equiv g \pmod{p}$ и $\gcd(a, n) = 1$. Тогда $a^{n-1} \equiv 1 \pmod{p}$, т.е. $g^{n-1} \equiv 1 \pmod{p}$, т.е. $\text{ord}_p(g) = p-1 | n-1$. $\square$

($\Leftarrow$) Пусть $n$ squarefree и $(p-1) | (n-1)$ для всех $p | n$. Для любого $a$ с $\gcd(a, n) = 1$: по каждому $p | n$ имеем $a^{n-1} = a^{(p-1) \cdot k} = (a^{p-1})^k \equiv 1^k = 1 \pmod{p}$ (т.к. $(p-1) | (n-1)$, $k = (n-1)/(p-1)$). По КТО: $a^{n-1} \equiv 1 \pmod{n}$. $\square$

Следствие. Число Кармайкла имеет не менее трёх различных простых делителей.

Доказательство. Пусть $n = p \cdot q$ с $p < q$ простыми. Тогда $(q-1) | (n-1)$. Но $n-1 = pq - 1$. Если $(q-1) | (pq-1)$, то $(q-1) | (pq - 1 - p(q-1)) = pq - 1 - pq + p = p - 1$. Значит $q - 1 \le p - 1$, т.е. $q \le p$ — противоречие с $p < q$. $\square$

Наименьшее число Кармайкла

$561 = 3 \cdot 11 \cdot 17$. Проверим критерий Корселта:

• Squarefree: $3, 11, 17$ — различные простые. ✓
• $(3-1) | 560$: $2 | 560$. ✓
• $(11-1) | 560$: $10 | 560$. ✓ ($560 = 56 \cdot 10$)
• $(17-1) | 560$: $16 | 560$. ✓ ($560 = 35 \cdot 16$)

Сильно псевдопростые числа (SPRP)

Тест Ферма недостаточен. Усилим его, используя свойство квадратных корней из 1.

Лемма. Если $p$ — нечётное простое и $x^2 \equiv 1 \pmod{p}$, то $x \equiv \pm 1 \pmod{p}$.

Доказательство. $(x-1)(x+1) \equiv 0 \pmod{p}$, и $p$ — простое, значит $p | (x-1)$ или $p | (x+1)$. $\square$

Алгоритм Миллера–Рабина. Пусть $n - 1 = 2^s \cdot d$, где $d$ нечётно. Число $n$ проходит тест Миллера–Рабина по основанию $a$ (является strong pseudoprime to base $a$), если выполнено одно из:

• $a^d \equiv 1 \pmod{n}$, или
• $a^{2^r d} \equiv -1 \pmod{n}$ для некоторого $r \in \{0, 1, \ldots, s-1\}$.

Теорема. Если $p$ — нечётное простое и $\gcd(a, p) = 1$, то $p$ всегда проходит тест Миллера–Рабина по основанию $a$.

Доказательство. $a^{p-1} \equiv 1 \pmod{p}$ (Ферма). Разложим: $a^{p-1} = (a^d)^{2^s}$. Рассмотрим последовательность $a^d, a^{2d}, a^{4d}, \ldots, a^{2^s d} = a^{p-1} \equiv 1$. Первый член $\equiv 1$: условие выполнено. Иначе первый член $\equiv -1$: условие выполнено. Иначе $\exists r$: $a^{2^r d} \not\equiv \pm 1$ но $a^{2^{r+1} d} \equiv 1$. Тогда $(a^{2^r d})^2 \equiv 1 \pmod{p}$, но $a^{2^r d} \not\equiv \pm 1$ — противоречие с леммой. $\square$

Теорема (Рабин, 1980). Если $n$ составное, то доля оснований $a \in [1, n-1]$, по которым $n$ проходит тест МР, не превосходит $1/4$.

Следствие: тест МР с $k$ случайными основаниями ошибается с вероятностью $\le 4^{-k}$.

Детерминированный тест Миллера–Рабина

Для детерминированного теста нужно знать фиксированный набор оснований, гарантирующий правильный ответ для всех $n < B$.

Теорема (Payam Pourali, 2020; сводная таблица):

Для задач соревновательного программирования (числа до $10^{18}$) достаточно оснований $\{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37\}$ или просто $\{2, 3, 7, 61, 24251\}$ (достаточно для $n < 3.2 \cdot 10^{18}$).

BPSW тест

Тест Бейли–Псевдо–Простых–Числа–Вайсса (BPSW) — комбинация:

1. Тест МР по основанию 2.
2. Тест Лукаса с сильным условием (strong Lucas pseudoprime test).

Ни одного BPSW-псевдопростого числа до $10^{15}$ не известно. Предположительно, их нет до $10^{30}$ и выше. На практике BPSW считается надёжным для любых $n$ в пределах вычислительных возможностей.

Числа Кармайкла и тест МР

Числа Кармайкла проходят тест Ферма по любому основанию, взаимно простому с $n$. Но они не проходят тест Миллера–Рабина по большинству оснований.

Пример. $n = 561 = 3 \cdot 11 \cdot 17$, $n - 1 = 560 = 2^4 \cdot 35$.

Тест МР по основанию $a = 2$: $2^{35} \pmod{561}$: $2^{35} = 2^{32} \cdot 2^3 = 4294967296 \cdot 8$. $2^{35} \bmod 561 = 263$. $263 \ne 1$ и $263 \ne 560$ ($= -1 \bmod 561$). $2^{70} \bmod 561 = 263^2 \bmod 561 = 69169 \bmod 561 = 166$. $166 \ne 560$. $2^{140} \bmod 561 = 166^2 \bmod 561 = 27556 \bmod 561 = 67$. $67 \ne 560$. $2^{280} \bmod 561 = 67^2 \bmod 561 = 4489 \bmod 561 = 1$. $1 = 1$, но ни на каком предыдущем шаге не было $-1$.

По условию теста МР: $a^d = 2^{35} \equiv 263 \not\equiv 1$ и ни одно $a^{2^r d}$ не равно $-1$ (хотя последнее равно $1$). Значит, $n = 561$ проваливает тест МР по основанию 2 — правильно выявлен как составное!

---

Ключевые формулы (сводная таблица)

---

Реализация на C++20 (с анализом сложности)

#include <bits/stdc++.h>
using namespace std;
using ll  = long long;
using ull = unsigned long long;
using u128 = unsigned __int128;

// ============================================================
// Модульное умножение без переполнения (n до 2^63)
// ============================================================
ll mulmod(ll a, ll b, ll m){
    return (__int128)a * b % m;
}

// Быстрое возведение в степень по модулю
ll powmod(ll a, ll b, ll m){
    ll res = 1;
    a %= m;
    if(a < 0) a += m;
    for(; b > 0; b >>= 1){
        if(b & 1) res = mulmod(res, a, m);
        a = mulmod(a, a, m);
    }
    return res;
}

// ============================================================
// Тест Миллера–Рабина: проверяет n по основанию a
// Возвращает true если n проходит тест (возможно простое)
// ============================================================
bool miller_rabin_test(ll n, ll a){
    if(n % a == 0) return n == a;

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

    ll x = powmod(a, d, n);
    if(x == 1 || x == n - 1) return true;

    for(int r = 1; r < s; r++){
        x = mulmod(x, x, n);
        if(x == n - 1) return true;
    }
    return false;
}

// ============================================================
// Детерминированный тест простоты
// Корректен для n < 3,317,044,064,679,887,385,961,981
// (покрывает все числа, встречающиеся в задачах)
// Время: O(log^2 n) — 12 итераций MR
// ============================================================
bool is_prime(ll n){
    if(n < 2) return false;
    if(n == 2 || n == 3 || n == 5 || n == 7) return true;
    if(n % 2 == 0 || n % 3 == 0 || n % 5 == 0) return false;

    // Для n < 3,215,031,751 достаточно {2, 3, 5, 7}
    // Для n < 3,317,044,064,679,887,385,961,981 — следующий набор:
    for(ll a : {2LL, 3LL, 5LL, 7LL, 11LL, 13LL, 17LL, 19LL, 23LL, 29LL, 31LL, 37LL}){
        if(n == a) return true;
        if(!miller_rabin_test(n, a)) return false;
    }
    return true;
}

// Компактная версия для n < 3.2 * 10^18 (достаточно для long long)
bool is_prime_fast(ll n){
    if(n < 2) return false;
    for(ll a : {2LL, 3LL, 5LL, 7LL, 11LL, 13LL, 17LL, 19LL, 23LL, 29LL, 31LL, 37LL}){
        if(n == a) return true;
        if(n % a == 0) return false;
    }
    for(ll a : {2LL, 3LL, 5LL, 7LL, 11LL, 13LL, 17LL, 19LL, 23LL, 29LL, 31LL, 37LL}){
        if(!miller_rabin_test(n, a)) return false;
    }
    return true;
}

// ============================================================
// Алгоритм Полларда–Ро для факторизации
// Время: O(n^{1/4}) в среднем
// ============================================================
ll pollard_rho(ll n){
    if(n % 2 == 0) return 2;
    if(n % 3 == 0) return 3;

    // Случайные значения для нескольких попыток
    static mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());

    while(true){
        ll x = rng() % (n - 2) + 2;
        ll y = x;
        ll c = rng() % (n - 1) + 1;
        ll d = 1;

        while(d == 1){
            x = (mulmod(x, x, n) + c) % n;
            y = (mulmod(y, y, n) + c) % n;
            y = (mulmod(y, y, n) + c) % n;
            d = __gcd(abs(x - y), n);
        }
        if(d != n) return d;
    }
}

// Полная факторизация n
void factorize(ll n, map<ll, int>& factors){
    if(n == 1) return;
    if(is_prime(n)){
        factors[n]++;
        return;
    }
    ll d = pollard_rho(n);
    factorize(d, factors);
    factorize(n / d, factors);
}

// ============================================================
// Проверка: является ли n числом Кармайкла?
// ============================================================
bool is_carmichael(ll n){
    if(n < 2 || is_prime(n)) return false;
    map<ll, int> f;
    factorize(n, f);
    // Проверка squarefree
    for(auto [p, e] : f)
        if(e > 1) return false;
    // Проверка (p-1) | (n-1)
    for(auto [p, e] : f)
        if((n - 1) % (p - 1) != 0) return false;
    return true;
}

// ============================================================
// Нахождение k-го числа Кармайкла
// ============================================================
ll kth_carmichael(int k){
    ll n = 3;
    int cnt = 0;
    while(true){
        if(is_carmichael(n)) { cnt++; if(cnt == k) return n; }
        n += 2; // все числа Кармайкла нечётны
    }
}

// ============================================================
// Оптимизированный is_prime с предварительной проверкой
// малых делителей (даёт ускорение в 2-5 раз на практике)
// ============================================================
bool is_prime_optimized(ll n){
    if(n < 2) return false;
    if(n < 4) return true;
    if(n % 2 == 0 || n % 3 == 0) return false;

    // Пробные деления до min(sqrt(n), 1000)
    for(ll i = 5; i <= min((ll)1000, (ll)sqrtl((long double)n)); i += (i % 6 == 5) ? 2 : 4){
        if(n % i == 0) return false;
    }
    if(n <= 1000000LL) {
        // Уже проверили все делители до 1000 >= sqrt(n)
        return true;
    }

    // Тест МР
    for(ll a : {2LL, 3LL, 5LL, 7LL, 11LL, 13LL, 17LL, 19LL, 23LL, 29LL, 31LL, 37LL}){
        if(n == a) return true;
        if(n % a == 0) return false;
        if(!miller_rabin_test(n, a)) return false;
    }
    return true;
}

// ============================================================
// Сильный тест Лукаса (для BPSW)
// ============================================================

// Символ Якоби (a/n)
ll jacobi(ll a, ll n){
    a %= n;
    ll res = 1;
    while(a != 0){
        while(a % 2 == 0){
            a /= 2;
            if(n % 8 == 3 || n % 8 == 5) res = -res;
        }
        swap(a, n);
        if(a % 4 == 3 && n % 4 == 3) res = -res;
        a %= n;
    }
    return (n == 1) ? res : 0;
}

// Выбор параметра D для теста Лукаса (метод Сельфриджа)
ll find_D(ll n){
    ll D = 5;
    while(true){
        ll j = jacobi(D, n);
        if(j == 0 && abs(D) < n) return 0; // n составное (общий делитель)
        if(j == -1) return D;
        D = (D > 0) ? -(D + 2) : -(D - 2);
    }
}

// Последовательности Лукаса U_k(P, Q), V_k(P, Q) по модулю n
// P, Q определяются из D: P = 1, Q = (1 - D) / 4
struct Lucas {
    ll U, V, Qk;
};

// Вычисляет U_{k} и V_{k} по модулю n через удвоение
Lucas lucas_sequence(ll P, ll Q, ll k, ll n){
    // Используем алгоритм удвоения:
    // U_{2m} = U_m * V_m (mod n)
    // V_{2m} = V_m^2 - 2*Q^m (mod n)
    ll U = 1, V = P, Qk = Q;
    ll U2 = 0, V2 = 2; // U_0 = 0, V_0 = 2
    // Бит за битом k:
    ll D = mulmod(P, P, n) - 4 * Q % n;
    D = ((D % n) + n) % n;

    ll Um = 0, Vm = 2, Qm = 1;
    ll Uk = 1, Vk = P, Qkk = Q % n;

    int bits = 63 - __builtin_clzll(k);
    for(int i = bits; i >= 0; i--){
        // Удвоение: (Um, Vm, Qm) -> (U_{2m}, V_{2m}, Q^{2m})
        ll new_U = mulmod(Um, Vm, n);
        ll new_V = (mulmod(Vm, Vm, n) - 2 * Qm % n + 2*n) % n;
        ll new_Q = mulmod(Qm, Qm, n);
        Um = new_U; Vm = new_V; Qm = new_Q;

        if((k >> i) & 1){
            // Шаг: (Um, Vm) -> (U_{m+1}, V_{m+1})
            ll nu = (mulmod(P, Um, n) + Vm) % n;
            if(nu % 2 != 0) nu += n;
            nu /= 2;
            ll nv = (mulmod(D % n, Um, n) + mulmod(P, Vm, n)) % n;
            if(nv % 2 != 0) nv += n;
            nv /= 2;
            Um = nu; Vm = nv; Qm = mulmod(Qm, Q, n);
        }
    }
    return {Um, Vm, Qm};
}

// Сильный тест Лукаса
bool strong_lucas_test(ll n){
    if(n < 2) return false;
    if(n == 2) return true;
    if(n % 2 == 0) return false;

    ll D = find_D(n);
    if(D == 0) return false; // gcd(D, n) > 1

    ll P = 1;
    ll Q = (1 - D) / 4;

    // n + 1 = 2^s * d
    ll m = n + 1;
    int s = 0;
    while(m % 2 == 0){ m /= 2; s++; }

    // Вычислить U_d, V_d, Q^d
    auto [Ud, Vd, Qd] = lucas_sequence(P, Q, m, n);

    if(Ud == 0 || Vd == 0) return true;

    ll Qk = Qd;
    for(int r = 1; r < s; r++){
        // V_{2^r d} = V_{2^{r-1} d}^2 - 2 * Q^{2^{r-1} d}
        Vd = (mulmod(Vd, Vd, n) - 2 * Qk % n + 2*n) % n;
        Qk = mulmod(Qk, Qk, n);
        if(Vd == 0) return true;
    }
    return false;
}

// BPSW тест: практически безошибочный для n < 10^{30}
bool is_prime_bpsw(ll n){
    if(n < 2) return false;
    if(n == 2 || n == 3 || n == 5) return true;
    if(n % 2 == 0 || n % 3 == 0 || n % 5 == 0) return false;

    // Шаг 1: МР по основанию 2
    if(!miller_rabin_test(n, 2)) return false;

    // Шаг 2: сильный тест Лукаса
    return strong_lucas_test(n);
}

int main(){
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    // Тесты
    cout << boolalpha;
    cout << "is_prime(561) = " << is_prime(561) << "\n";      // false (Кармайкла)
    cout << "is_prime(1729) = " << is_prime(1729) << "\n";    // false
    cout << "is_prime(998244353) = " << is_prime(998244353) << "\n"; // true
    cout << "is_prime(1000000007) = " << is_prime(1000000007) << "\n"; // true

    // Числа Кармайкла
    cout << "\nЧисла Кармайкла до 10000:\n";
    for(ll n = 3; n <= 10000; n += 2)
        if(is_carmichael(n)) cout << n << " ";
    cout << "\n";
    // Ожидается: 561 1105 1729 2465 2821 6601 8911

    // Факторизация
    map<ll,int> f;
    factorize(561LL, f);
    cout << "\n561 = ";
    for(auto [p, e] : f) cout << p << "^" << e << " * ";
    cout << "\n";

    // Проверка для n = 10^18 + 9 (простое)
    ll big_prime = 1000000000000000003LL;
    cout << "is_prime(" << big_prime << ") = " << is_prime(big_prime) << "\n";

    return 0;
}

Анализ сложности:

Для $n \le 10^{18}$: is_prime работает за наносекунды (12 · 60 = 720 умножений).

---

Разбор задачи 1 (средняя): Подсчёт простых в массиве

Условие. Дан массив из $10^5$ чисел $a_i \le 10^{18}$. Найти количество простых.

Решение. Применяем is_prime к каждому элементу.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
// ... (код is_prime из выше)

int main(){
    int n; cin >> n;
    int cnt = 0;
    for(int i = 0; i < n; i++){
        ll x; cin >> x;
        if(is_prime(x)) cnt++;
    }
    cout << cnt << "\n";
    // Сложность: O(n * log^2(max_a)) = O(10^5 * 60^2) ~ 3.6 * 10^8 mulmod
    // Но mulmod быстрая (~3 нс), итого ~1 сек
}

Оптимизация: сначала проверить делимость на малые простые (2, 3, 5, ..., 997). Это отсеивает $\approx 95\%$ составных числел без дорогостоящего МР.

---

Разбор задачи 2 (сложная): Факторизация и мультипликативная функция

Условие. Дано $Q \le 10^5$ запросов. Для каждого $n \le 10^{18}$ найти $\varphi(n)$.

Решение. Факторизуем $n$ через Полларда–Ро и вычисляем $\varphi(n) = n \prod_{p|n} (1 - 1/p)$.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
// ... (код is_prime, pollard_rho, factorize)

ll euler_phi(ll n){
    map<ll,int> f;
    factorize(n, f);
    ll result = n;
    for(auto [p, e] : f){
        result -= result / p;
    }
    return result;
}

int main(){
    int Q; cin >> Q;
    while(Q--){
        ll n; cin >> n;
        cout << euler_phi(n) << "\n";
    }
}

Тест: $\varphi(561) = \varphi(3) \cdot \varphi(11) \cdot \varphi(17) = 2 \cdot 10 \cdot 16 = 320$. Проверка критерия Корселта: $(3-1)|(560), (11-1)|(560), (17-1)|(560)$: $2|560, 10|560, 16|560$ — всё верно.

Сложность: $O(Q \cdot n^{1/4} \log^2 n)$ в среднем. При $Q = 10^5$, $n = 10^{18}$: $10^5 \cdot 10^{4.5} \cdot 3600 \approx 10^{12.7}$ — слишком медленно.

Оптимизация: для каждого $n$ разбиваем: если после пробных делений до $10^6$ остаток $r = 1$, то готово; если $r$ простое (тест МР) — готово; иначе один вызов Полларда–Ро. Реально для большинства $n$ пробные деления дают один–два делителя быстро.

Практически для $Q = 10^5$ и $n \le 10^{18}$: $< 2$ сек.

---

Разбор задачи 3 (очень сложная — CF Div.1 E / ICPC уровень): Нахождение числа Кармайкла в диапазоне

Условие. Дано $T \le 1000$ запросов. Для каждого найти наименьшее число Кармайкла $\ge L$ для $L \le 10^{15}$.

Анализ. Прямой перебор всех нечётных чисел от $L$ до первого числа Кармайкла — неэффективен (расстояние между числами Кармайкла может быть велико при $L \approx 10^{15}$).

Метод: используем критерий Корселта для генерации чисел Кармайкла.

Числа Кармайкла вида $n = p_1 p_2 p_3$ с тремя простыми делителями строятся следующим образом. Зафиксируем $\lambda = \text{lcm}(p_1-1, p_2-1, p_3-1)$, тогда $n = (k\lambda + 1)$ для некоторого $k$, и $n = p_1 p_2 p_3$.

Конструктивный алгоритм Эрдёша для трёхпростых Кармайкла:

// Генерация всех чисел Кармайкла вида p*q*r до N
// за O(N^{2/3} / log^2 N) при наивном подходе

vector<ll> generate_carmichael_3primes(ll N) {
    vector<ll> result;
    // Перебираем p < q, затем находим r
    // Из (p-1)|(n-1), (q-1)|(n-1), (r-1)|(n-1)
    // n = pqr, n-1 = k*lcm(p-1,q-1,r-1)

    // Упрощение: p*(p-1) | (n-1), фиксируем p, q
    // l = lcm(p-1, q-1), n-1 = l*k для некоторого k
    // n = pqr = l*k + 1
    // r = (l*k+1)/(p*q)

    // Простые до cbrt(N)
    int limit = (int)cbrtl((long double)N) + 1;
    vector<int> primes;
    vector<bool> isp(limit+1, true);
    for(int i=2; i<=limit; i++){
        if(isp[i]) primes.push_back(i);
        for(int j=0; j<(int)primes.size()&&(ll)i*primes[j]<=limit; j++){
            isp[i*primes[j]]=false;
            if(i%primes[j]==0) break;
        }
    }

    for(int pi = 0; pi < (int)primes.size(); pi++){
        ll p = primes[pi];
        for(int qi = pi+1; qi < (int)primes.size(); qi++){
            ll q = primes[qi];
            ll L = lcm(p-1, q-1);
            // Ищем k: r = (k*L+1)/(p*q) должно быть целым и простым
            // k*L + 1 ≡ 0 (mod p*q)
            // k ≡ -1/L (mod p*q/gcd(L, p*q))
            ll pq = p * q;
            // Хотим k*L ≡ -1 (mod pq), т.е. k ≡ -L^{-1} mod pq
            // Проверяем: gcd(L, pq) должен делить 1 => gcd(L, pq) = 1
            ll g = __gcd(L, pq);
            if(g != 1) continue; // нет решений

            // k0 = обратное к L по модулю pq, умноженное на -1
            // Используем расширенный Евклид
            // ... (реализация inv_mod)
            // k ≡ k0 (mod pq), k = k0 + pq*t, t >= 0
            // r = (k*L+1)/(pq), k должно быть >= 1

            // Для краткости: наивный перебор k
            ll k = 1;
            while(true){
                ll num = k * L + 1;
                if(num > N) break;
                if(num % pq == 0){
                    ll r = num / pq;
                    if(r > q && is_prime(r)){
                        // Проверяем (r-1) | (n-1) = (k*L)
                        if(k * L % (r-1) == 0){
                            result.push_back(num);
                        }
                    }
                }
                k++;
            }
        }
    }

    sort(result.begin(), result.end());
    result.erase(unique(result.begin(), result.end()), result.end());
    return result;
}

Более эффективный подход для соревнований: предгенерировать все числа Кармайкла до $10^{15}$ (их $\approx 105000$) и отвечать на запросы бинарным поиском. Генерация через алгоритм Эрдёша занимает $O(B^{2/3} \log^2 B)$ при грамотной реализации.

// Для соревнования: предподсчёт числе Кармайкла
vector<ll> carmichael_list;

void precompute(ll limit){
    // Используем generate_carmichael_3primes (основная часть)
    // + числа с 4+ простыми делителями (их значительно меньше)
    carmichael_list = generate_carmichael_3primes(limit);
    sort(carmichael_list.begin(), carmichael_list.end());
}

ll next_carmichael(ll L){
    auto it = lower_bound(carmichael_list.begin(), carmichael_list.end(), L);
    if(it == carmichael_list.end()) return -1; // нет в диапазоне
    return *it;
}

Тест: первые числа Кармайкла: 561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341.

---

Задачи для самостоятельного решения (таблица)

---

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

1. Применение теста Ферма вместо Миллера–Рабина. Числа Кармайкла проходят тест Ферма по всем основаниям. Если вы пишете if(pow(a, n-1, n) == 1) return true — ваш код неверен для чисел Кармайкла.

2. Переполнение в mulmod. При $n > 2^{31}$ произведение $a \cdot b$ может переполнить long long. Всегда используйте (__int128)a * b % m или аналог.

3. Неправильный выбор оснований. Набор $\{2, 3\}$ не достаточен для $n \ge 1{,}373{,}654$. Конкретный контрпример: $n = 1{,}373{,}653 = 389 \cdot 3533$ проходит МР по основаниям 2 и 3, но составное. Всегда используйте проверенный набор оснований из таблицы.

4. Граничный случай $n = 2$. Некоторые реализации МР возвращают false для $n = 2$, т.к. $n - 1 = 1$ — нечётное, и $s = 0$. Добавьте явную проверку маленьких простых.

5. Алгоритм Полларда зацикливается. При неудачном выборе $c$ (например, $c = 0$ или $c = n-2$) алгоритм может зациклиться. Добавьте проверку или используйте несколько итераций с разными $c$.

6. Неверная проверка чисел Кармайкла для чисел с $p^2 | n$. Критерий Корселта требует squarefree. Если $n = 4$ или $n = 12$, то $4 | n$, и условие $(p-1) | (n-1)$ может выполняться для некоторых $p$, но $n$ не Кармайкла. Всегда проверяйте squarefree первым.

---

Совет профессионала

В реальных соревнованиях для проверки простоты одиночного числа $n \le 10^{18}$ используйте детерминированный тест МР с основаниями $\{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37\}$ — он гарантированно верен и работает за наносекунды. Не используйте случайные основания: в соревновательных условиях жюри могут специально подобрать "тяжёлые" входные данные.

Для факторизации $n \le 10^{18}$ комбинируйте: (1) пробные деления до $10^6$ (отсекает большинство случаев за $10^6 / \ln(10^6) \approx 72000$ итераций); (2) тест простоты МР для остатка; (3) Полларда–Ро для составного остатка. Такая комбинация работает $< 1$ мс для любого $n \le 10^{18}$.

---

Итог

Числа Кармайкла — составные числа, "обманывающие" тест Ферма. Их характеризует критерий Корселта: squarefree и $(p-1) | (n-1)$ для каждого простого делителя $p$.

Тест Миллера–Рабина — усиление теста Ферма через проверку квадратных корней из 1. Доля "обманных" оснований для составного $n$ не превосходит $1/4$.

Детерминированный тест с основаниями $\{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37\}$ корректен для всех $n < 3.3 \cdot 10^{24}$ — более чем достаточно для любых задач в long long и __int128.

Практические рекомендации:

1. Одиночная проверка простоты $n \le 10^{18}$: детерминированный МР, $O(\log^2 n)$.
2. Факторизация $n \le 10^{18}$: Полларда–Ро + МР, $O(n^{1/4} \log^2 n)$.
3. Никогда не используйте тест Ферма без МР — числа Кармайкла его обходят.
4. BPSW = МР(2) + Лукас — практически надёжен, используется в Python's isprime.