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

Факторизация больших чисел — одна из центральных задач как в криптографии (безопасность RSA), так и в соревновательном программировании. Когда число n имеет порядок 10^12–10^18, пробное деление работает за O(√n) операций, что может занять несколько секунд. Алгоритм Полларда ρ позволяет найти нетривиальный делитель за O(n^{1/4}) итераций, но и он пасует перед числами порядка 10^20 и выше.

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

1. Метод квадратичного решета (Quadratic Sieve, QS) — субэкспоненциальный алгоритм, работающий за exp(O(√(ln n · ln ln n))) операций. Один из лучших практических алгоритмов для чисел до 100 цифр.

2. Метод эллиптических кривых (ECM) Ленстры — вероятностный алгоритм, сложность которого зависит от наименьшего простого делителя p числа n, что делает его идеальным для поиска "средних" множителей (10–50 цифр).

Умение выбрать правильный инструмент — ключевой навык:

---

Теория

Идея конгруэнции квадратов

Основу большинства современных методов факторизации составляет следующая идея. Если удастся найти два числа x и y такие, что:

x² ≡ y²  (mod n),   x ≢ ±y  (mod n)

то gcd(x − y, n) даст нетривиальный делитель n. Это следует из того, что n | (x−y)(x+y), но n не делит ни x−y, ни x+y по условию.

Доказательство. Пусть d = gcd(x−y, n). Тогда d | n. Если d = n, то n | x−y, т.е. x ≡ y (mod n) — противоречие. Если d = 1, то n | x+y (из n | (x−y)(x+y) при gcd(n, x−y)=1), т.е. x ≡ −y (mod n) — тоже противоречие. Значит, 1 < d < n, что и требовалось. □

B-гладкие числа

Число называется B-гладким, если все его простые делители не превышают B. Факторная база — множество всех простых p ≤ B.

Идея QS. Выберем последовательность чисел вида Q(x) = (⌈√n⌉ + x)² − n. При x ≈ 0 эти числа малы по абсолютной величине (порядка √n · x), поэтому среди них встречается много B-гладких.

Если собрать достаточно B-гладких значений Q(x_i), каждое из которых разложено по факторной базе:

Q(x_i) = p_1^{e_{i,1}} · p_2^{e_{i,2}} · ... · p_k^{e_{i,k}}

то задача сводится к поиску подмножества, в котором все показатели степеней чётны. Тогда произведение выбранных Q(x_i) будет точным квадратом.

Алгоритм квадратичного решета (шаги)

Шаг 1: Выбор параметров. Факторная база F = {p простое : p ≤ B}. Оптимальное B ≈ exp(½ · √(ln n · ln ln n)).

Шаг 2: Просеивание. Для x от 0 до некоторого M вычисляем Q(x) = (⌈√n⌉ + x)² − n. Используем решетообразное просеивание (как решето Эратосфена): для каждого простого p ∈ F находим корни сравнения x² ≡ n (mod p) и вычитаем log(p) из журнала для всех позиций x в этих арифметических прогрессиях. Позиция считается B-гладкой, если журнал близок к нулю.

Шаг 3: Построение матрицы. Для каждого гладкого Q(x_i) строим вектор показателей степеней по модулю 2. Это вектор в GF(2)^k, где k = |F|.

Шаг 4: Линейная алгебра над GF(2). Методом Гаусса (или более быстрыми алгоритмами для разреженных матриц) находим ненулевые линейные зависимости, т.е. подмножества строк с нулевой суммой в GF(2).

Шаг 5: Восстановление делителя. Для найденного подмножества {i_1, ..., i_t}:

x = ∏ (⌈√n⌉ + x_{i_j})  (mod n)
y = √(∏ Q(x_{i_j}))      (mod n)

Тогда x² ≡ y² (mod n). Вычисляем gcd(x−y, n).

Почему QS работает быстро

Теорема (оценка вероятности гладкости). Доля B-гладких чисел среди чисел до N асимптотически равна ρ(u), где u = ln N / ln B, а ρ — функция Диксона–Канфилда. При u = 2 (т.е. B = √N), ρ(2) = 1 − ln 2 ≈ 0.307.

Для QS нам нужно собрать O(|F|) ≈ O(B) гладких значений. Каждый Q(x_i) по величине ≈ M√n. Вероятность гладкости ≈ ρ(ln(M√n)/ln(B)). При оптимальном выборе B и M суммарное время просеивания и линейной алгебры даёт субэкспоненциальную сложность.

Метод эллиптических кривых Ленстры (ECM)

Эллиптическая кривая над Z/nZ задаётся уравнением Вейерштрасса:

y² = x³ + ax + b  (mod n),   4a³ + 27b² ≢ 0  (mod n)

Множество точек кривой вместе с "точкой на бесконечности" O образует абелеву группу с законом сложения, задаваемым геометрически (хорда-касательная).

Закон сложения. Для точек P = (x_1, y_1) и Q = (x_2, y_2) на кривой над полем:

• Если x_1 = x_2 и y_1 = −y_2, то P + Q = O
• Иначе λ = (y_2 − y_1)/(x_2 − x_1) (или (3x_1² + a)/(2y_1) при P = Q)
• x_3 = λ² − x_1 − x_2, y_3 = λ(x_1 − x_3) − y_1

Идея ECM. Работаем над Z/nZ, но n составное. Пусть p — неизвестный простой делитель n. Группа точек E(F_p) имеет порядок m_p ≈ p ± 2√p (теорема Хассе). Если m_p имеет только малые простые делители (т.е. m_p B-гладкое), то:

k = lcm(1, 2, ..., B) = ∏_{p_i ≤ B} p_i^{⌊log_{p_i} B⌋}

При вычислении k·P над Z/nZ на каком-то шаге знаменатель (x_2 − x_1 или 2y) окажется кратным p, но не n. Тогда gcd(знаменатель, n) даст нетривиальный делитель.

Алгоритм ECM:

1. Выбрать случайные a и точку P = (x_0, y_0) над Z/nZ. Вычислить b = y_0² − x_0³ − ax_0 (mod n).
2. Вычислить Q = k · P методом удвоения-сложения.
3. На каждом шаге: при вычислении обратного по модулю n вызвать gcd с n. Если gcd ≠ 1 и gcd ≠ n, найден делитель.
4. Если не найдено, повторить с другой кривой.

Почему мы меняем кривые? Если m_p не B-гладкое, алгоритм неудачен. Но у разных кривых разные m_p, и вероятность того, что хотя бы одна из T кривых даст гладкий порядок, быстро растёт с T.

Сравнение методов

Pollard ρ (O(p^{1/2})): эффективен для малых p (до ~20 знаков). Прост в реализации.

ECM (exp(O(√(ln p · ln ln p)))): сложность зависит от p — наименьшего делителя. Идеален для задач, где у n есть делитель 20–50 знаков.

QS (exp(O(√(ln n · ln ln n)))): сложность зависит от n — числа целиком. Лучше для "трудных" случаев — произведения двух простых сопоставимого размера.

---

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

---

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

Pollard ρ с улучшением Брента

#include <bits/stdc++.h>
using namespace std;
using u64 = unsigned long long;
using u128 = __uint128_t;

// Умножение по модулю без переполнения
u64 mulmod(u64 a, u64 b, u64 m) {
    return (__uint128_t)a * b % m;
}

// Возведение в степень по модулю
u64 powmod(u64 a, u64 b, u64 m) {
    u64 res = 1;
    a %= m;
    while (b > 0) {
        if (b & 1) res = mulmod(res, a, m);
        a = mulmod(a, a, m);
        b >>= 1;
    }
    return res;
}

// Тест Миллера-Рабина (детерминированный для n < 3,317,044,064,679,887,385,961,981)
bool is_prime(u64 n) {
    if (n < 2) return false;
    if (n == 2 || n == 3 || n == 5 || n == 7) return true;
    if (n % 2 == 0 || n % 3 == 0) return false;
    u64 d = n - 1;
    int r = 0;
    while (d % 2 == 0) { d /= 2; r++; }
    // Базисы, достаточные для n < 3.3 * 10^24
    for (u64 a : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
        if (a >= n) continue;
        u64 x = powmod(a, d, n);
        if (x == 1 || x == n - 1) continue;
        bool composite = true;
        for (int i = 0; i < r - 1; i++) {
            x = mulmod(x, x, n);
            if (x == n - 1) { composite = false; break; }
        }
        if (composite) return false;
    }
    return true;
}

// Алгоритм Брента (модификация Полларда ρ)
u64 brent(u64 n, u64 x0 = 2, u64 c = 1) {
    u64 x = x0, y = x0, d = 1;
    u64 m = 128;
    u64 q = 1, ys, r = 1;
    while (d == 1) {
        y = x;
        for (u64 i = 0; i < r; i++)
            x = (mulmod(x, x, n) + c) % n;
        u64 k = 0;
        while (k < r && d == 1) {
            ys = y;
            for (u64 i = 0; i < min(m, r - k); i++) {
                y = (mulmod(y, y, n) + c) % n;
                q = mulmod(q, x > y ? x - y : y - x, n);
            }
            d = __gcd(q, n);
            k += m;
        }
        r *= 2;
    }
    if (d == n) {
        while (true) {
            ys = (mulmod(ys, ys, n) + c) % n;
            d = __gcd(x > ys ? x - ys : ys - x, n);
            if (d > 1) break;
        }
    }
    return d == n ? 0 : d;
}

// Полная факторизация числа n
map<u64, int> factorize(u64 n) {
    map<u64, int> factors;
    if (n == 1) return factors;
    // Очередь для обработки
    queue<u64> q;
    q.push(n);
    while (!q.empty()) {
        u64 x = q.front(); q.pop();
        if (x == 1) continue;
        if (is_prime(x)) {
            factors[x]++;
            continue;
        }
        // Пробуем найти делитель
        u64 d = x;
        u64 c = 1;
        while (d == x) {
            d = brent(x, 2, c++);
        }
        q.push(d);
        q.push(x / d);
    }
    return factors;
}

Анализ сложности Pollard ρ. Ожидаемое число итераций — O(p^{1/2}), где p — наименьший простой делитель n. Итерация занимает O(log n) из-за mulmod. Итого: O(n^{1/4} · log n) для числа n = p·q, p ≈ q.

ECM: упрощённая реализация

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

// Точка на эллиптической кривой y² = x³ + ax + b (mod n)
struct Point {
    u64 x, y;
    bool inf = false; // точка на бесконечности
};

struct Curve {
    u64 a, n; // кривая y² = x³ + ax + b, работаем mod n
};

// Обратный элемент: возвращает {inv, gcd}
pair<u64,u64> extended_gcd_mod(u64 a, u64 m) {
    // Расширенный алгоритм Евклида
    long long old_r = a, r = m;
    long long old_s = 1, s = 0;
    while (r != 0) {
        long long q = old_r / r;
        tie(old_r, r) = make_pair(r, old_r - q * r);
        tie(old_s, s) = make_pair(s, old_s - q * s);
    }
    if (old_r < 0) old_r += m; // нормализация
    return {(u64)((old_s % (long long)m + m) % m), (u64)old_r};
}

u64 mulmod(u64 a, u64 b, u64 m) {
    return (__uint128_t)a * b % m;
}

// Сложение точек на кривой над Z/nZ
// Возвращает {результат, нетривиальный делитель n или 0}
pair<Point, u64> add_points(Point P, Point Q, const Curve& c) {
    if (P.inf) return {Q, 0};
    if (Q.inf) return {P, 0};
    u64 n = c.n;
    u64 lam_num, lam_den;
    if (P.x == Q.x) {
        if ((P.y + Q.y) % n == 0) return {Point{0, 0, true}, 0}; // P + (-P) = O
        // Касательная: λ = (3x² + a) / (2y)
        lam_num = (mulmod(3, mulmod(P.x, P.x, n), n) + c.a) % n;
        lam_den = 2 * P.y % n;
    } else {
        // Хорда: λ = (y2 - y1) / (x2 - x1)
        lam_num = (Q.y >= P.y) ? (Q.y - P.y) : (n - P.y + Q.y);
        lam_den = (Q.x >= P.x) ? (Q.x - P.x) : (n - P.x + Q.x);
    }
    auto [inv, g] = extended_gcd_mod(lam_den, n);
    if (g != 1) {
        // Нашли нетривиальный делитель!
        return {{}, (g != n) ? g : 0};
    }
    u64 lam = mulmod(lam_num, inv, n);
    u64 x3 = (mulmod(lam, lam, n) - P.x % n - Q.x % n + 2 * n) % n;
    u64 y3 = (mulmod(lam, (P.x >= x3 ? P.x - x3 : n - x3 + P.x), n) - P.y + n) % n;
    return {Point{x3, y3}, 0};
}

// Умножение точки на скаляр k методом удвоения
// Возвращает {результат, делитель или 0}
pair<Point, u64> scalar_mult(u64 k, Point P, const Curve& c) {
    Point result{0, 0, true};
    Point base = P;
    while (k > 0) {
        if (k & 1) {
            auto [r, d] = add_points(result, base, c);
            if (d) return {{}, d};
            result = r;
        }
        auto [r, d] = add_points(base, base, c);
        if (d) return {{}, d};
        base = r;
        k >>= 1;
    }
    return {result, 0};
}

// ECM: попытка найти делитель n
// B1 — граница первой фазы
u64 ecm_once(u64 n, u64 B1, mt19937_64& rng) {
    // Случайная кривая и точка
    u64 x = rng() % n;
    u64 y = rng() % n;
    u64 a = rng() % n;
    u64 b = (mulmod(y, y, n) - mulmod(x, mulmod(x, x, n), n) - mulmod(a, x, n) % n + 3 * n) % n;
    // Проверка невырожденности: gcd(4a³ + 27b², n)
    u64 disc = (mulmod(4, mulmod(a, mulmod(a, a, n), n), n) +
                mulmod(27, mulmod(b, b, n), n)) % n;
    u64 g = __gcd(disc, n);
    if (g != 1 && g != n) return g;
    if (g == n) return 0;

    Curve curve{a, n};
    Point P{x, y};

    // Первая фаза: умножаем P на k = ∏ p^e ≤ B1
    // Генерируем простые до B1 решетом
    vector<bool> sieve(B1 + 1, true);
    sieve[0] = sieve[1] = false;
    for (u64 i = 2; i * i <= B1; i++)
        if (sieve[i])
            for (u64 j = i*i; j <= B1; j += i)
                sieve[j] = false;

    for (u64 p = 2; p <= B1; p++) {
        if (!sieve[p]) continue;
        u64 pk = p;
        while (pk <= B1) pk *= p;
        pk /= p;
        auto [new_P, d] = scalar_mult(pk, P, curve);
        if (d) return d;
        P = new_P;
    }
    return 0;
}

// Полный ECM с несколькими попытками
u64 ecm_factor(u64 n, int attempts = 100, u64 B1 = 50000) {
    mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());
    for (int i = 0; i < attempts; i++) {
        u64 d = ecm_once(n, B1, rng);
        if (d && d != n) return d;
    }
    return 0; // не нашли
}

Анализ сложности ECM. При границе B1 вероятность успеха на одной кривой ≈ exp(−√(ln p · ln ln p) / √2), где p — наименьший делитель n. Чтобы найти делитель с вероятностью 1/2, нужно T ≈ exp(√(ln p · ln ln p) · √2 / 2) кривых. Каждая кривая требует O(B1 · log B1) операций. Оптимальное B1 = exp(√(ln p · ln ln p) / 2), что даёт общую сложность exp(√2 · √(ln p · ln ln p)).

---

Разбор задачи 1 (средняя)

Задача. Дано число n ≤ 10^12. Найти количество делителей. Тест: n = 963761198400.

Решение. Число умещается в 64 бита, делители малы. Применяем Pollard ρ.

963761198400 = 2^7 · 3^4 · 5^2 · 7 · 11 · 13 · 17 · 19

Количество делителей: (7+1)(4+1)(2+1)(1+1)^5 = 8·5·3·32 = 3840.

Код применяет factorize(n) из раздела выше, затем:

long long num_divisors(long long n) {
    auto f = factorize(n);
    long long res = 1;
    for (auto [p, e] : f) res *= (e + 1);
    return res;
}

Время: доли миллисекунды.

---

Разбор задачи 2 (сложная)

Задача. Дано n ≤ 10^18. Разложить на простые множители. Гарантируется, что n = p · q, где p и q — простые. Найти p и q.

Решение. Это типичное RSA-подобное число. Каждый из p, q ≈ 10^9. Pollard ρ с оптимизацией Брента справится быстро: O(p^{1/2}) ≈ O(10^{4.5}) ≈ 30000 итераций.

int main() {
    u64 n;
    cin >> n;
    if (is_prime(n)) {
        cout << n << " " << 1 << "\n";
        return 0;
    }
    auto f = factorize(n);
    for (auto [p, e] : f)
        cout << p << "^" << e << " ";
    cout << "\n";
}

Нюанс: при n = p² нужно корректно обработать случай brent возвращает p. В коде factorize это учтено — при d == x меняем c.

---

Разбор задачи 3 (ICPC WF / CF Div.1 E)

Задача (CF 385E / похожая структура). Дан граф с n вершинами. Каждой вершине v назначен вес w(v). Запрос: для вершин из подтерева, посчитать количество пар (u, v) с gcd(w(u), w(v)) > 1.

Задача сводится к: для массива весов подтерева подсчитать пары с нетривиальным gcd. Используем включение-исключение по простым делителям.

Ключевая часть: быстрая факторизация всех w(v). Веса могут достигать 10^12. Применяем Pollard ρ для каждого — O(n · n^{1/4} · log n) — что при n = 10^5 составляет ~10^7 операций.

// Для каждого числа получаем список различных простых делителей
vector<u64> prime_factors(u64 n) {
    auto f = factorize(n);
    vector<u64> res;
    for (auto [p, e] : f) res.push_back(p);
    return res;
}

// Мёбиус + включение-исключение для подсчёта пар с gcd > 1
// Для массива a: ответ = C(n,2) - количество пар с gcd = 1
// Количество пар с gcd = 1: sum_{d=1}^{max} mu(d) * C(cnt[d], 2)
// где cnt[d] = количество элементов, кратных d

---

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

---

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

1. Переполнение при умножении. Для n ~ 10^18 нельзя использовать a * b % n — переполнение u64. Обязательно использовать __uint128_t или __int128.

2. Pollard ρ зацикливается. Если brent вернул n (а не делитель), нужно менять параметр c. Классическая ошибка — не обрабатывать этот случай.

3. Не проверять простоту. Перед тем как запускать ECM/Pollard на числе, всегда проверяйте is_prime. Запуск Полларда на простом числе — бесконечный цикл.

4. ECM: вырожденная кривая. При случайном выборе a, b нужно проверять дискриминант 4a³ + 27b² ≢ 0 (mod p). На практике gcd(дискриминант, n) иногда сам даёт делитель.

5. Неверные пределы просеивания в QS. При выборе M слишком малым не наберём достаточно гладких чисел, слишком большим — много лишней работы.

---

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

Для задач соревновательного программирования реализовывать полный QS обычно не требуется — числа редко превышают 10^18. Достаточно Полларда с Брентом + Miller-Rabin. Однако ECM в упрощённом варианте иногда полезен для разложения чисел 10^15–10^18 с "неудобными" делителями. Хорошая точка разрыва: если Pollard не нашёл делитель за 10^6 итераций, попробуйте ECM с B1 = 1000 на 50 кривых. Если и это не помогает, число скорее всего простое (проверьте Miller-Rabin ещё раз).

Для самых серьёзных случаев (n до 10^60) в CP иногда нужен QS. Рекомендую изучить реализацию msieve или yafu — они демонстрируют все описанные техники с оптимизациями.

---

Итог

Мы изучили два мощных метода факторизации: QS и ECM. Метод квадратичного решета основан на построении конгруэнций квадратов через B-гладкие числа и линейную алгебру над GF(2). ECM Ленстры использует случайные эллиптические кривые и зависит от наименьшего делителя — что делает его идеальным для "средних" множителей. Для практических задач в CP достаточен Pollard ρ с Брентом, но знание ECM открывает дорогу к задачам более высокого уровня. Главный принцип остаётся неизменным: любая современная факторизация опирается на конгруэнцию квадратов.