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

Разложение числа в произведение простых — одна из фундаментальных задач теории чисел. Почти всё, что связано с мультипликативными функциями, делителями, функцией Эйлера, сравнениями и НОД, в той или иной мере опирается на факторизацию. Но она же лежит в основе криптографии RSA: сложность взлома RSA эквивалентна сложности факторизации произведения двух больших простых.

В соревновательном программировании задачи на факторизацию бывают трёх уровней сложности:

• Простой уровень: подсчёт делителей, сумма делителей для $n \leq 10^7$ — достаточно решета наименьших простых делителей (SPF).
• Средний уровень: факторизация одиночных чисел до $10^{12}$ — пробное деление достаточно.
• Сложный уровень: факторизация чисел до $10^{18}$ — нужен $\rho$-алгоритм Полларда.

Последнее встречается в задачах Div.1 D/E на Codeforces и финалах ICPC. Знание этого алгоритма — редкое преимущество, выделяющее топ-участников.

Теория

Факторизация через наименьший простой делитель

Идея линейного решета из главы 2 даёт нам массив $\text{spf}[i]$ — наименьший простой делитель числа $i$. Имея этот массив, факторизация любого числа $n \leq N$ занимает $O(\log n)$ времени.

Алгоритм факторизации через SPF:

Пока n > 1:
    p = spf[n]
    Считаем степень p: пока n % p == 0, делим n на p и увеличиваем показатель
    Добавляем (p, степень) к результату

Количество шагов не превышает $\log_2 n$ (каждый шаг уменьшает число минимум вдвое), то есть $O(\log n)$.

Применение: подсчёт делителей, мультипликативных функций для всех чисел до $N$ одновременно. После построения SPF за $O(N)$ каждый запрос на факторизацию — $O(\log n)$.

Пробное деление для одиночных чисел

Для факторизации одиночного числа $n$ проверяем делимость на все числа от 2 до $\sqrt{n}$

​. Сложность $O(\sqrt{n})$. Для $n \leq 10^{12}$ это $10^6$ операций — приемлемо. Для $n \leq 10^{18}$ — $10^9$ операций, слишком медленно.

Оптимизация: проверяем 2, 3, затем числа вида $6k \pm 1$. Ускорение в $3$ раза.

Алгоритм Полларда (ρ-алгоритм)

История и идея. Джон Ми Поллард предложил алгоритм в 1975 году. Название «ρ» (ро) происходит от формы траектории: последовательность $x_0, x_1, x_2, \ldots$ в конце концов зацикливается, образуя «хвост» и «петлю» — фигуру, похожую на греческую букву ρ.

Центральная идея — парадокс дней рождения. Если случайно выбирать числа от 1 до $n$, то ожидаемое количество выборов до первого повторения — $O(\sqrt{n})$

​). Аналогично, если $p$ — делитель числа $N$, то среди случайных чисел $x_1, \ldots, x_k$ ожидается коллизия по модулю $p$ (то есть $x_i \equiv x_j \pmod{p}$, откуда $\gcd(|x_i - x_j|, N) > 1$) уже при $k = O(\sqrt{p}) = O(N^{1/4})$.

Конкретный алгоритм. Строим последовательность: $x_{i+1} = (x_i^2 + c) \bmod N$ для случайно выбранного $c \neq 0, -2$. Это псевдослучайная последовательность. По парадоксу дней рождения, через $O(\sqrt{p})$

​) шагов (где $p$ — наименьший нетривиальный делитель $N$) появится пара $x_i \equiv x_j \pmod{p}$. Тогда $\gcd(|x_i - x_j|, N)$ даёт нетривиальный делитель.

Нахождение цикла — алгоритм Флойда. Используем двух «черепах»: медленная $x$ делает один шаг, быстрая $y$ — два шага за итерацию. Когда последовательность входит в цикл по модулю $p$, быстрая черепаха нагонит медленную за $O(\sqrt{p})$

​) итераций. В этот момент $\gcd(|x - y|, N) > 1$.

Улучшение Брента. Ричард Брент в 1980 году предложил более эффективный вариант, который практически в 25-36% быстрее алгоритма Флойда. Идея: вместо алгоритма Флойда используем «степенной» цикл. Фиксируем значение $y$ и проверяем блок за блоком длиной $2^k$.

В реализации Брента мы накапливаем произведение разностей (через умножение по модулю $N$) и берём НОД от всего произведения, а не от каждой разности по отдельности. Это экономит дорогие операции $\gcd$.

Оценка сложности. Для числа $N$ с наименьшим простым делителем $p \leq \sqrt{N}$

​:

• $\rho$-алгоритм Полларда находит $p$ за ожидаемое $O(p^{1/2}) = O(N^{1/4})$ итераций.
• Каждая итерация: одно умножение и одно сложение по модулю $N$ — $O(\log N)$ битовых операций (или $O(1)$ при 64-битных числах с аппаратным умножением).
• Итого: $O(N^{1/4} \log N)$ в среднем.

Для $N \sim 10^{18}$: $N^{1/4} \approx 10^{4.5} \approx 31623$. Это около $3 \times 10^4$ итераций, то есть почти мгновенно!

Полная факторизация рекурсивно применяет алгоритм Полларда:

1. Если $N = 1$, закончить.
2. Если $N$ простое (проверить тестом Миллера-Рабина), добавить к результату.
3. Иначе найти нетривиальный делитель $d$ алгоритмом Полларда.
4. Рекурсивно факторизовать $d$ и $N/d$.

Подсчёт делителей и их суммы через факторизацию

Теорема. Если $n = p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}$, то:

1. Количество делителей: $\tau(n) = (a_1 + 1)(a_2 + 1) \cdots (a_k + 1)$

2. Сумма делителей: $\sigma(n) = \prod_{i=1}^{k} \frac{p_i^{a_i+1} - 1}{p_i - 1}$

3. Функция Эйлера: $\varphi(n) = n \prod_{p \mid n} \left(1 - \frac{1}{p}\right) = \prod_{i=1}^{k} p_i^{a_i - 1}(p_i - 1)$

Доказательство формулы для $\tau(n)$. Каждый делитель $d$ числа $n$ однозначно определяется выбором показателей: $d = p_1^{b_1} \cdots p_k^{b_k}$, где $0 \leq b_i \leq a_i$. Для $p_i^{a_i}$ есть $a_i + 1$ вариантов (от 0 до $a_i$). По правилу произведения: $\tau(n) = \prod (a_i + 1)$. $\square$

Доказательство формулы для $\sigma(n)$. Сумма делителей мультипликативна: $\sigma(p^a) = 1 + p + p^2 + \cdots + p^a = (p^{a+1}-1)/(p-1)$. По мультипликативности: $\sigma(n) = \prod \sigma(p_i^{a_i})$. $\square$

Мультипликативные функции — функции $f$ такие, что $f(1)=1$ и $f(mn)=f(m)f(n)$ при $\gcd(m,n)=1$. Все функции $\tau, \sigma, \varphi$ мультипликативны. Это ключевое свойство позволяет вычислять их через факторизацию.

Подсчёт делителей для всех чисел до $N$ (решето делителей): Можно обойтись без явной факторизации: для каждого $d$ от 1 до $N$ увеличиваем счётчик всех кратных $d$. Сложность $O(N \log N)$ (гармонический ряд).

Число делителей и задачи на делители

Теорема (Оценка $\tau(n)$). Для любого $\varepsilon > 0$: $\tau(n) = O(n^\varepsilon)$. На практике для $n \leq 10^{18}$ максимальное $\tau(n) \approx 100000$ (высококомпозитные числа).

Несколько конкретных значений:

• $\tau(720720) = 240$ (720720 = $2^4 \cdot 3^2 \cdot 5 \cdot 7 \cdot 11 \cdot 13$)
• Наибольшее $\tau(n)$ для $n \leq 10^6$: $\tau(720720) = 240$
• Наибольшее $\tau(n)$ для $n \leq 10^{18}$: около $103680$ (у числа $963761198400 \cdot \ldots$)

Знание этой оценки важно: если задача требует перебора всех делителей числа, то при $n \leq 10^{18}$ это не более $10^5$ делителей — перебор допустим после факторизации.

Ключевые формулы

╔═══════════════════════════════════════════════════════════════════╗
║  СВОДКА: ФАКТОРИЗАЦИЯ                                             ║
╠═══════════════════════════════════════════════════════════════════╣
║  n = p₁^a₁ · p₂^a₂ · ... · pₖ^aₖ                                   ║
║  τ(n) = (a₁+1)(a₂+1)···(aₖ+1)           кол-во делителей          ║
║  σ(n) = ∏ (pᵢ^(aᵢ+1) - 1) / (pᵢ - 1)   сумма делителей            ║
║  φ(n) = n · ∏ (1 − 1/pᵢ)                функция Эйлера            ║
║  Пробное деление: O(√n)                                           ║
║  Решето SPF + факторизация: O(N) + O(log n)                       ║
║  Алгоритм Полларда (Брент): O(n^{1/4} log n) ожидаемое            ║
║  Тест простоты МР: O(k log² n), k=12 для детерминированного       ║
╚═══════════════════════════════════════════════════════════════════╝

Реализация на C++20

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

// ============================================================
// 1. Тест простоты Миллера-Рабина (из главы 2, повторяем)
//    Необходим как субпроцедура для алгоритма Полларда
// ============================================================
ll mulmod(ll a, ll b, ll m) {
    return (__int128)a * b % m;
}

ll powmod(ll base, ll exp, ll mod) {
    ll result = 1;
    base %= mod;
    if (base < 0) base += mod;
    while (exp > 0) {
        if (exp & 1) result = mulmod(result, base, mod);
        base = mulmod(base, base, mod);
        exp >>= 1;
    }
    return result;
}

bool miller_rabin_round(ll n, ll a) {
    if (n % a == 0) return n == a;
    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;
}

bool is_prime(ll n) {
    if (n < 2) return false;
    if (n == 2 || n == 3) return true;
    if (n % 2 == 0 || n % 3 == 0) return false;
    for (ll a : {2LL, 3LL, 5LL, 7LL, 11LL, 13LL, 17LL, 19LL,
                 23LL, 29LL, 31LL, 37LL}) {
        if (a >= n) continue;
        if (!miller_rabin_round(n, a)) return false;
    }
    return true;
}

// ============================================================
// 2. Факторизация через SPF (наименьший простой делитель)
//    Предподсчёт за O(N), факторизация одного числа за O(log n)
// ============================================================
const int SPF_MAX = 1 << 20; // 10^6, настраивайте под задачу
vector<int> spf; // spf[i] = наименьший простой делитель i

void build_spf(int n = SPF_MAX) {
    spf.resize(n + 1);
    iota(spf.begin(), spf.end(), 0); // spf[i] = i изначально

    for (int i = 2; (ll)i * i <= n; i++) {
        if (spf[i] == i) { // i — простое (не было перезаписано)
            for (int j = i * i; j <= n; j += i) {
                if (spf[j] == j) // ещё не встречали делитель
                    spf[j] = i;
            }
        }
    }
}

// Факторизация n <= SPF_MAX за O(log n)
map<int, int> factorize_spf(int n) {
    map<int, int> factors;
    while (n > 1) {
        int p = spf[n];
        while (n % p == 0) {
            factors[p]++;
            n /= p;
        }
    }
    return factors;
}

// ============================================================
// 3. Пробное деление для одиночных чисел
//    Работает для n <= ~10^12 (sqrt(n) ~= 10^6 итераций)
// ============================================================
map<ll, int> factorize_trial(ll n) {
    map<ll, int> factors;

    // Делители 2 и 3 — отдельно
    for (ll p : {2LL, 3LL}) {
        while (n % p == 0) {
            factors[p]++;
            n /= p;
        }
    }

    // Числа вида 6k±1, от 5 до sqrt(n)
    for (ll p = 5; p * p <= n; p += 6) {
        for (ll q : {p, p + 2}) {
            while (n % q == 0) {
                factors[q]++;
                n /= q;
            }
        }
    }

    if (n > 1) factors[n]++; // n — оставшийся простой делитель
    return factors;
}

// ============================================================
// 4. ρ-алгоритм Полларда с улучшением Брента
//    Находит НЕТРИВИАЛЬНЫЙ делитель n (не обязательно простой).
//    Возвращает n, если делитель не найден (перезапустить с другим c).
//    Ожидаемое время: O(n^{1/4} log n)
// ============================================================
ll pollard_rho(ll n, ll c) {
    if (n % 2 == 0) return 2;

    // Стартовое значение
    ll x = rand() % (n - 2) + 2;
    ll y = x;
    ll d = 1;

    // Функция итерации: f(x) = (x^2 + c) mod n
    auto f = [&](ll x) {
        return (mulmod(x, x, n) + c) % n;
    };

    // Алгоритм Брента
    ll ys, r = 1, q = 1;
    ll product; // будем накапливать произведение для batch-НОД

    do {
        // y — «быстрый» указатель, x — «медленный» (фиксируется)
        // x остаётся фиксированным на протяжении цикла длиной r
        ll saved_x = x; // сохраняем x для возврата при неудаче
        for (ll i = 0; i < r; i++) y = f(y); // двигаем y на r шагов

        ll k = 0;
        while (k < r && d == 1) {
            ys = y;
            // Накапливаем произведение |x - y| mod n
            // Берём НОД только раз на M шагов (batch) — ускорение
            product = 1;
            ll M = min((ll)128, r - k); // размер пакета
            for (ll i = 0; i < M; i++) {
                y = f(y);
                product = mulmod(product, abs(x - y), n);
            }
            d = __gcd(product, n);
            k += M;
        }
        r *= 2;
    } while (d == 1);

    // Если d == n, попробуем пошаговый откат
    if (d == n) {
        d = 1;
        y = ys;
        while (d == 1) {
            y = f(y);
            d = __gcd(abs(x - y), n);
        }
    }

    return d == n ? n : d; // n означает неудачу — нужно перезапустить
}

// ============================================================
// 5. Полная факторизация за O(n^{1/4} log n)
//    Рекурсивно раскладывает n через Полларда + тест Миллера-Рабина
// ============================================================
void factorize_pollard(ll n, map<ll, int>& factors) {
    if (n == 1) return;

    if (is_prime(n)) {
        factors[n]++;
        return;
    }

    // Ищем нетривиальный делитель
    ll d = n;
    for (ll c = 1; d == n || d == 1; c++) {
        // Перебираем разные c до успеха
        // В среднем нужно O(1) попыток
        d = pollard_rho(n, c);
    }

    // Рекурсивно факторизуем оба делителя
    factorize_pollard(d, factors);
    factorize_pollard(n / d, factors);
}

map<ll, int> factorize_fast(ll n) {
    map<ll, int> factors;

    // Сначала выносим маленькие делители (ускоряет Полларда)
    for (ll p : {2LL, 3LL, 5LL, 7LL, 11LL, 13LL}) {
        while (n % p == 0) {
            factors[p]++;
            n /= p;
        }
    }

    if (n > 1) factorize_pollard(n, factors);
    return factors;
}

// ============================================================
// 6. Вычисление функций через факторизацию
// ============================================================

// Количество делителей τ(n)
ll count_divisors(const map<ll, int>& factors) {
    ll result = 1;
    for (auto [p, e] : factors)
        result *= (e + 1);
    return result;
}

// Сумма делителей σ(n)
ll sum_divisors(const map<ll, int>& factors) {
    ll result = 1;
    for (auto [p, e] : factors) {
        ll pk = 1, sum = 1;
        for (int i = 0; i < e; i++) {
            pk *= p;
            sum += pk;
        }
        result *= sum;
    }
    return result;
}

// Функция Эйлера φ(n)
ll euler_phi(const map<ll, int>& factors) {
    ll result = 1;
    for (auto [p, e] : factors) {
        result *= (p - 1);         // p^1 - p^0 = p-1
        for (int i = 1; i < e; i++)
            result *= p;           // p^(e-1)
    }
    return result;
}

// Получить все делители n (отсортированные)
// Работает для n <= 10^18, если τ(n) разумное (< 10^5)
vector<ll> get_divisors(const map<ll, int>& factors) {
    vector<ll> divs = {1};
    for (auto [p, e] : factors) {
        int sz = divs.size();
        ll pe = 1;
        for (int i = 0; i < e; i++) {
            pe *= p;
            for (int j = 0; j < sz; j++)
                divs.push_back(divs[j] * pe);
        }
    }
    sort(divs.begin(), divs.end());
    return divs;
}

// ============================================================
// 7. Решето делителей — подсчёт τ(i) для всех i <= N за O(N log N)
// ============================================================
vector<int> divisor_count_sieve(int n) {
    vector<int> cnt(n + 1, 0);
    for (int d = 1; d <= n; d++)
        for (int m = d; m <= n; m += d)
            cnt[m]++;
    return cnt;
}

// Более быстрая версия через SPF и мультипликативность
// τ(p^a * m) = (a+1) * τ(m) если gcd(p^a, m) = 1
// Используем: τ(i) = τ(spf[i]^k) * τ(i / spf[i]^k)

// ============================================================
// Пример использования
// ============================================================
int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    srand(42); // для воспроизводимости

    // Тест пробного деления
    auto f1 = factorize_trial(360);
    cout << "360 = ";
    for (auto [p, e] : f1) cout << p << "^" << e << " * ";
    cout << "\n"; // 2^3 * 3^2 * 5^1

    // Функции от 360
    cout << "τ(360) = " << count_divisors(f1) << "\n"; // 24
    cout << "σ(360) = " << sum_divisors(f1) << "\n";   // 1170
    cout << "φ(360) = " << euler_phi(f1) << "\n";      // 96

    // Все делители 360
    auto divs = get_divisors(f1);
    cout << "Делители 360: ";
    for (ll d : divs) cout << d << " ";
    cout << "\n";

    // Факторизация большого числа алгоритмом Полларда
    ll big = 1234567891011121314LL; // специально составное
    auto f2 = factorize_fast(big);
    cout << big << " = ";
    for (auto [p, e] : f2) cout << p << "^" << e << " * ";
    cout << "\n";

    // Проверка: произведение должно совпасть
    ll check = 1;
    for (auto [p, e] : f2) {
        ll pe = 1;
        for (int i = 0; i < e; i++) pe *= p;
        check *= pe;
    }
    cout << "Проверка: " << check << " == " << big << ": "
         << (check == big ? "OK" : "ОШИБКА") << "\n";

    // SPF-решето
    build_spf();
    auto f3 = factorize_spf(720720);
    cout << "720720 = ";
    for (auto [p, e] : f3) cout << p << "^" << e << " * ";
    cout << "\n"; // 2^4 * 3^2 * 5 * 7 * 11 * 13
    cout << "τ(720720) = " << count_divisors(f3) << "\n"; // 240

    // Решето подсчёта делителей
    auto cnt = divisor_count_sieve(20);
    cout << "τ(i) для i=1..20: ";
    for (int i = 1; i <= 20; i++) cout << cnt[i] << " ";
    cout << "\n";

    return 0;
}

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

Разбор задачи 1 — Число делителей (простая задача)

Постановка. Дано $T$ чисел $n_i$ ($1 \leq n_i \leq 10^6$). Для каждого числа вывести количество его делителей.

Идея. Предподсчёт с помощью SPF или решета делителей. SPF строим за $O(N)$, затем каждый запрос — за $O(\log n)$.

#include <bits/stdc++.h>
using namespace std;

const int MAXN = 1e6 + 5;
int spf[MAXN];

void build_spf() {
    iota(spf, spf + MAXN, 0);
    for (int i = 2; (long long)i * i < MAXN; i++)
        if (spf[i] == i) // i простое
            for (int j = i * i; j < MAXN; j += i)
                if (spf[j] == j)
                    spf[j] = i;
}

// Подсчёт делителей через SPF
int count_divs(int n) {
    int result = 1;
    while (n > 1) {
        int p = spf[n], e = 0;
        while (n % p == 0) { n /= p; e++; }
        result *= (e + 1);
    }
    return result;
}

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

    build_spf();

    // Кэшируем результаты, если запросов много
    vector<int> tau(MAXN, 0);
    for (int i = 1; i < MAXN; i++)
        tau[i] = count_divs(i);

    int T;
    cin >> T;
    while (T--) {
        int n;
        cin >> n;
        cout << tau[n] << "\n";
    }
    return 0;
}

Альтернатива: если запросов мало, можно не кэшировать, вызывать count_divs(n) напрямую.

Замечание. Решето делителей (for d: for m=d,2d,...: cnt[m]++) за $O(N \log N)$ тоже даёт правильный ответ, но чуть медленнее для $N = 10^6$. Разница незначительна, но SPF дополнительно даёт факторизацию.

Разбор задачи 2 — Сумма мультипликативной функции (средняя задача)

Постановка. Дано $n$ ($1 \leq n \leq 10^7$). Вычислить: $S = \sum_{i=1}^{n} \varphi(i)$ Ответ вывести по модулю $10^9 + 7$.

Идея. Функция Эйлера мультипликативна. Её можно вычислить для всех чисел до $n$ аналогично решету.

Формула: для числа $p^k$ (простого степени): $\varphi(p^k) = p^{k-1}(p-1)$.

Линейное решето Эйлера: ключевое свойство линейного решета — при умножении $i \cdot p$ на простое $p$:

• Если $p \mid i$: $\varphi(i \cdot p) = \varphi(i) \cdot p$
• Если $p \nmid i$: $\varphi(i \cdot p) = \varphi(i) \cdot (p - 1)$

#include <bits/stdc++.h>
using namespace std;

const int MOD = 1e9 + 7;
const int MAXN = 1e7 + 5;

int phi_arr[MAXN];
int primes[700000]; // примерно столько простых до 10^7
int prime_cnt = 0;
bool not_prime[MAXN];

void compute_phi(int n) {
    phi_arr[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!not_prime[i]) {
            // i — простое
            primes[prime_cnt++] = i;
            phi_arr[i] = i - 1; // φ(p) = p - 1
        }
        for (int j = 0; j < prime_cnt; j++) {
            int p = primes[j];
            if ((long long)i * p > n) break;
            not_prime[i * p] = true;
            if (i % p == 0) {
                // p | i, значит p^2 | i*p
                // φ(i*p) = φ(i) * p
                phi_arr[i * p] = phi_arr[i] * p;
                break; // важно: прерываем, чтобы сохранить линейность
            } else {
                // gcd(p, i) = 1
                // φ(i*p) = φ(i) * φ(p) = φ(i) * (p-1)
                phi_arr[i * p] = phi_arr[i] * (p - 1);
            }
        }
    }
}

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

    int n;
    cin >> n;

    compute_phi(n);

    long long S = 0;
    for (int i = 1; i <= n; i++)
        S = (S + phi_arr[i]) % MOD;

    cout << S << "\n";
    return 0;
}

Анализ. Линейное решето работает за $O(N)$ времени. Сумма вычисляется за ещё один проход $O(N)$. Память: $O(N)$.

Дополнительная оптимизация: на больших $N$ (до $10^8$) используют решето с битовыми масками и блочную обработку для лучшей работы кэша процессора.

Формула для суммы Эйлера: известно, что $\sum_{i=1}^{n} \varphi(i) = \frac{1}{2}\left(1 + \sum_{i=1}^{n} \mu(i) \lfloor n/i \rfloor^2\right)$, что позволяет вычислить ответ за $O(n^{2/3})$ с помощью решета Люси Хедж (Lucy_Hedgehog). Но это уже тема продвинутых математических трюков.

Разбор задачи 3 — Сложная задача уровня Div.1 E

Постановка (CF 1033G, адаптация). Дано $n$ чисел $a_1, \ldots, a_n$ ($2 \leq a_i \leq 10^{18}$). Для каждого числа найдите его наибольший простой делитель. Выведите сумму всех наибольших простых делителей.

Алгоритм:

1. Для каждого $a_i$: применить тест Миллера-Рабина.
2. Если простое — оно само является своим наибольшим делителем.
3. Если составное — рекурсивно факторизовать через алгоритм Полларда.
4. Выбрать максимальный простой делитель.

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

// [mulmod, powmod, is_prime — как выше]
ll mulmod(ll a, ll b, ll m) { return (__int128)a * b % m; }

ll powmod(ll base, ll exp, ll mod) {
    ll result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = mulmod(result, base, mod);
        base = mulmod(base, base, mod);
        exp >>= 1;
    }
    return result;
}

bool miller_round(ll n, ll a) {
    if (n % a == 0) return n == a;
    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;
}

bool is_prime(ll n) {
    if (n < 2) return false;
    for (ll a : {2LL, 3LL, 5LL, 7LL, 11LL, 13LL, 17LL, 19LL, 23LL, 29LL, 31LL, 37LL})
        if (a < n && !miller_round(n, a)) return false;
    return true;
}

// Алгоритм Полларда-Брента
mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());

ll pollard_brent(ll n) {
    if (n % 2 == 0) return 2;
    if (n % 3 == 0) return 3;

    while (true) {
        // Случайные начальные значения
        ll x = rng() % (n - 2) + 2;
        ll c = rng() % (n - 1) + 1;
        ll y = x, d = 1;

        auto f = [&](ll v) { return (mulmod(v, v, n) + c) % n; };

        ll r = 1, q = 1, ys;
        do {
            ll saved = y;
            for (ll i = 0; i < r; i++) y = f(y);
            ll k = 0;
            while (k < r && d == 1) {
                ys = y;
                ll M = min((ll)128, r - k);
                for (ll i = 0; i < M; i++) {
                    y = f(y);
                    q = mulmod(q, abs(x - y), n);
                }
                d = __gcd(q, n);
                k += M;
            }
            x = saved; // Брент: x не двигается, а перешагивает через r
            // Исправление: в алгоритме Брента x обновляется
            // Используем стандартную схему:
            x = y; // упрощённый Брент для корректности
            r *= 2;
        } while (d == 1);

        if (d == n) {
            // Откат к последней сохранённой точке
            d = 1;
            y = ys;
            while (d == 1) {
                y = f(y);
                d = __gcd(abs(x - y), n);
            }
        }
        if (d != n) return d; // нашли нетривиальный делитель
        // иначе повторяем с новыми случайными значениями
    }
}

// Рекурсивная факторизация
ll max_prime_factor(ll n) {
    if (n == 1) return 1;
    if (is_prime(n)) return n;

    ll d = pollard_brent(n);
    return max(max_prime_factor(d), max_prime_factor(n / d));
}

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

    int T;
    cin >> T;

    ll total = 0;
    while (T--) {
        ll n;
        cin >> n;
        total += max_prime_factor(n);
    }
    cout << total << "\n";
    return 0;
}

Анализ производительности. Для чисел $n \leq 10^{18}$ с $n = 10^4$ запросами:

• В худшем случае (полупростые числа $n = p \cdot q$, $p \approx q \approx \sqrt{n} \approx 10^9$): $O(n^{1/4}) \approx 10^{4.5}$ итераций на число.
• Итого: $\sim 10^{4.5} \times 10^4 = 10^{8.5}$ операций — на грани 2 секунд.
• На практике алгоритм Полларда значительно быстрее из-за случайного характера: среднее число итераций существенно меньше теоретического максимума.

Оптимизации для конкурса:

1. Сначала отделить маленькие простые делители (2, 3, 5, 7, ...) пробным делением.
2. Использовать mt19937_64 вместо rand() для лучшего случайного поведения.
3. Размер пакета M = 128 — хороший баланс между частотой вычисления НОД и накоплением произведения.

Пример сложного теста: $n = 2^{60}$ — алгоритм Полларда сразу найдёт делитель 2. Или $n = p^2$ для большого простого $p$ — нужно найти $p$. Или $n = p \cdot q$ с $p \approx q$ — «сбалансированное полупростое», наихудший случай.

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

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

Ошибка 1: Бесконечный цикл в алгоритме Полларда.

// Если pollard_rho вернул n (неудача), нужно перезапустить!
// Неправильно:
ll d = pollard_rho(n, 1); // без проверки, что d != n

// Правильно:
ll d = n;
for (ll c = 1; d == n; c++)
    d = pollard_rho(n, c);

Ошибка 2: Переполнение в abs(x - y) при long long.

// x и y в диапазоне [0, n), разность может быть отрицательной
// abs() для long long: std::abs() или явное условие
ll diff = abs(x - y); // ПРАВИЛЬНО в C++ при #include <cstdlib>
// Или безопаснее: (x > y ? x - y : y - x)

Ошибка 3: Неправильная инициализация phi_arr[1].

// φ(1) = 1 по определению
phi_arr[1] = 1; // не забываем!
// Линейное решето начинается с i = 2, phi[1] нужно задать вручную

Ошибка 4: Целочисленное переполнение в sum_divisors для больших n.

// σ(n) может быть значительно больше n
// Для n = 2^60 ≈ 10^18: σ(2^60) = 2^61 - 1 ≈ 2*10^18 > LLONG_MAX
// Всегда проверяйте, не нужен ли __int128 или модульная арифметика

Ошибка 5: Не проверяем is_prime перед вызовом Полларда.

// pollard_rho НЕ работает с простыми числами (вернёт само число бесконечно)
// Обязательно:
if (is_prime(n)) { factors[n]++; return; }
// и только потом вызываем Полларда

Ошибка 6: Неправильный SPF для чисел-квадратов.

// spf[4] = 2, spf[9] = 3 — правильно.
// spf[1] = 1 — не является простым!
// Цикл факторизации должен иметь условие n > 1, а не n >= 1
while (n > 1) { // правильно
    int p = spf[n];
    ...
}

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

Наиболее распространённая ловушка с алгоритмом Полларда — случай $d = n$ («неудача»). Это происходит примерно в 5-10% запусков. Правильная стратегия — перезапуск с другим значением константы $c$ в функции итерации $f(x) = x^2 + c$.

Значения $c = 0$ и $c = -2$ нужно пропускать: при $c = 0$ последовательность вырождается ($0 \to 0$), при $c = -2$ — тоже.

Для максимальной скорости в конкурсных условиях:

// Используйте случайный c при каждом перезапуске:
mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());

ll find_factor(ll n) {
    while (true) {
        ll c = rng() % (n - 2) + 1; // c в [1, n-2]
        if (c == n - 2) c = 1; // избегаем c = -2 mod n
        ll d = pollard_rho_with_c(n, c);
        if (d != n && d != 1) return d;
    }
}

Для чисел специального вида (например, $n = 2^k - 1$ — числа Мерсенна) алгоритм Полларда может работать медленно. В таких случаях помогает предварительное пробное деление на первые несколько тысяч простых.

Полезный трюк для задач с несколькими числами: сортируйте числа по убыванию. Простые числа обрабатываются очень быстро (один тест МР), и часто половина чисел окажется простыми.

Итог

Факторизация — мост между общими числами и их простыми «атомами». Мы изучили три уровня сложности инструментов:

1. SPF-решето ($O(N)$ предобработка, $O(\log n)$ запрос) — стандарт для задач с $n \leq 10^7$ и множественными запросами. Даёт факторизацию и все мультипликативные функции.
2. Пробное деление ($O(\sqrt{n})$) — достаточно для $n \leq 10^{12}$. Просто в реализации и надёжно.
3. Алгоритм Полларда ($O(n^{1/4} \log n)$ ожидаемое) — необходим для $n \leq 10^{18}$. Требует теста Миллера-Рабина как субпроцедуры. При правильной реализации факторизует $10^{18}$ менее чем за миллисекунду.

Знание формул $\tau(n) = \prod (a_i + 1)$, $\sigma(n) = \prod (p_i^{a_i+1}-1)/(p_i-1)$ и $\varphi(n) = n \prod (1-1/p_i)$ позволяет мгновенно решать задачи на делители после факторизации.