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

В предыдущих главах мы видели мощные формулы: функция Мёбиуса, свёртка Дирихле, тождества между φ, d, σ. Всё это прекрасно теоретически — но насколько быстро мы можем вычислить эти функции на практике?

Наивный подход: для каждого n от 1 до N разложить n на простые и применить мультипликативность. Разложение за O(√n), итого O(N√N) — медленно при N = 10⁶: это 10⁹ операций.

Решето Эратосфена позволяет найти простые за O(N log log N). Зная простые, можно вычислить φ(n), μ(n), d(n), σ(n) для всех n ≤ N за O(N log N) через суммирование по делителям. Для N = 10⁶ это 2·10⁷ — приемлемо.

Но существует алгоритм, называемый линейным решетом (linear sieve), который:

• Вычисляет любую мультипликативную функцию за O(N) — строго линейно
• Помечает каждое составное число ровно один раз (в отличие от решета Эратосфена)
• Попутно вычисляет несколько функций одновременно без дополнительных затрат

При N = 10⁷ разница между O(N log log N) и O(N) заметна: 5·10⁷ против 10⁷ операций. При N = 10⁸ линейное решето становится необходимостью.

Теория

Идея: каждое число помечается ровно один раз

В классическом решете Эратосфена число n = 12 = 2²·3 вычёркивается при обработке простых 2 и 3 — дважды. В решете Эратосфена с флагом «если n/p — составное, прерваться» время работы становится O(N log log N), но пометка всё равно может происходить несколько раз.

Ключевая идея линейного решета: каждое составное число n «вычёркивается» ровно один раз — при умножении на его наименьший простой делитель.

Лемма. У каждого составного числа n существует единственное разложение n = p · m, где p = наименьший простой делитель n и m = n/p. При этом наименьший простой делитель m не меньше p.

Нет! Наименьший простой делитель m может быть как равен p (если p² | n), так и больше p. Правильная инвариантность линейного решета другая:

Инвариант линейного решета: число ip (i — текущее, p — простое из списка) помечается составным, когда p = наименьший простой делитель числа ip.

Почему это верно? При переборе p:

• Пока p < наименьший простой делитель(i): p < наим.пр.дел(i) ≤ i, значит наим.пр.дел(ip) = p. Это корректная пометка.
• Когда p = наименьший простой делитель(i): наим.пр.дел(ip) = p. Пометили корректно, затем прерываемся.

Если бы мы не прерывались, то следующий простой p' > p = наим.пр.дел(i), и число ip' имеет наименьший простой делитель p, а не p'. Тогда ip' уже будет помечено при i' = ip'/p, p — когда мы дойдём до этого i'.

Таким образом, каждое составное ip помечается ровно один раз: при именно том i, для которого p = наим.пр.дел(ip).

Теорема. Линейное решето выполняется за O(N) операций.

Доказательство. Каждое составное число n = i·p помечается ровно один раз. Простых чисел до N примерно N/ln N штук. Итого N − N/ln N + N/ln N = N меток, плюс O(N) итераций внешнего цикла. Общее число операций O(N). □

Алгоритм линейного решета (базовый)

primes = []
is_composite[1..N] = false
for i = 2 to N:
    if not is_composite[i]:
        primes.add(i)
    for p in primes:
        if i * p > N: break
        is_composite[i * p] = true
        if i % p == 0: break   // ключевая строка

Инвариант при работе с мультипликативными функциями. В момент обработки числа ip:

• Если i % p ≠ 0: p — новый простой в разложении ip, не встречавшийся в i. Тогда ip = i · p, НОД(i, p) = 1, и f(ip) = f(i) · f(p).
• Если i % p = 0: p уже входит в i. Пусть i = p^e · m, НОД(p,m)=1. Тогда ip = p^{e+1} · m.

Для второго случая нужно знать e и f(p^{e+1}). Храним дополнительно:

• pe[i] = значение p^e (наименьший простой делитель в степени e)
• exp_e[i] = показатель e

При i % p == 0:

• exp_e[ip] = exp_e[i] + 1
• pe[ip] = pe[i] · p

Формулы для конкретных функций при i % p == 0:

Пусть i = p^e · m, ip = p^{e+1} · m:

Формулы при i % p ≠ 0:

Пусть ip = i · p, НОД(i, p) = 1. Тогда f(ip) = f(i) · f(p):

Вычисление произвольной мультипликативной функции

Пусть задана мультипликативная функция f, определённая через:

• f(p) для любого простого p
• f(p^k) = g(p, k) — некоторая формула

Тогда в линейном решете:

• При i % p ≠ 0: f[ip] = f[i] · f[p] (значение f[p] вычислено при пометке p как простого)
• При i % p == 0: f[ip] = f[i] / g(p, exp_e[i]) · g(p, exp_e[i]+1) (делим на текущую p-компоненту и умножаем на следующую)

Важное замечание: Деление должно быть точным (целочисленным). Это гарантируется мультипликативностью: f[i] = f(p^{exp_e[i]}) · f(i / pe[i]), и f(p^{exp_e[i]}) делит f[i].

Доказательство корректности для φ

Проверим формулу φ[ip] = p · φ[i] при i % p == 0.

Пусть i = p^e · m, НОД(p,m)=1. Тогда:

• φ(i) = φ(p^e) · φ(m) = p^{e-1}(p-1) · φ(m)
• φ(ip) = φ(p^{e+1}) · φ(m) = p^e(p-1) · φ(m) = p · φ(p^e) · φ(m) = p · φ(i)

Формула верна. ✓

Доказательство корректности для d

При i = p^e · m:

• d(i) = (e+1) · d(m)
• d(ip) = (e+2) · d(m) = d(i) / (e+1) · (e+2)

Обозначим e+1 как exp_e[i] (это текущий показатель p в i, до умножения на p). Тогда:

• d[ip] = d[i] / exp_e[i] · (exp_e[i] + 1)

Но exp_e[ip] = exp_e[i] + 1, поэтому корректнее писать:

• d[ip] = d[i] / (exp_e[ip] - 1) · exp_e[ip]

Корректность для σ при i % p == 0

i = p^e · m, ip = p^{e+1} · m:

• σ(i) = σ(p^e) · σ(m)
• σ(ip) = σ(p^{e+1}) · σ(m) = σ(p^e · p + p^e + ... + 1 + p^{e+1}) ... аккуратнее:
• σ(p^{e+1}) = 1 + p + p² + ... + p^{e+1} = σ(p^e) + p^{e+1}
• σ(ip) = (σ(p^e) + p^{e+1}) · σ(m) = σ(i) · (σ(p^e) + p^{e+1}) / σ(p^e) = σ(i) / σ_pe[i] · (σ_pe[i] + pe[ip])

где σ_pe[i] = σ(p^e), pe[ip] = p^{e+1}. ✓

Сравнение с решетом Эратосфена

Практическое замечание. Из-за непредсказуемого обращения к памяти (vector primes с различными шагами) линейное решето на практике часто работает медленнее решета Эратосфена для простого нахождения простых. Его преимущество — возможность вычислять мультипликативные функции за O(N) при константах, сравнимых с O(N log log N) решетом.

Для N ≤ 10⁶ разница незначительна (оба работают за < 50ms). Для N = 10⁸ линейное решето — необходимость.

Ключевые формулы (сводка)

Инициализация:

f[1] = 1 (для всех мультипликативных функций)

Для простого i:

μ[i] = -1,  φ[i] = i-1,  d[i] = 2,  σ[i] = i+1
pe[i] = i,  exp_e[i] = 1,  σ_pe[i] = i+1

При i % p ≠ 0 (p — новый простой в ip):

μ[ip] = -μ[i]
φ[ip] = φ[i] · (p-1)
d[ip] = d[i] · 2
σ[ip] = σ[i] · (p+1)
pe[ip] = p, exp_e[ip] = 1, σ_pe[ip] = p+1

При i % p == 0 (p уже в i, показатель был e = exp_e[i]):

μ[ip] = 0
φ[ip] = φ[i] · p
d[ip] = d[i] / exp_e[i] · (exp_e[i]+1)    // = d[i]/e*(e+1)
σ_pe[ip] = σ_pe[i] + pe[i]·p              // σ(p^{e+1}) = σ(p^e) + p^{e+1}
σ[ip] = σ[i] / σ_pe[i] · σ_pe[ip]
pe[ip] = pe[i] · p, exp_e[ip] = exp_e[i]+1

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

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

// Полная реализация линейного решета с вычислением μ, φ, d, σ
struct LinearSieve {
    const int N;
    vector<int> mu, phi, d;
    vector<long long> sigma;
    vector<int> primes;
    vector<bool> is_composite;

    // Вспомогательные данные для корректного вычисления d и σ
    vector<int> pe;          // p^e: наименьший простой делитель в степени e
    vector<int> exp_e;       // e: показатель наименьшего простого делителя
    vector<long long> sigma_pe;  // σ(p^e)

    explicit LinearSieve(int N) : N(N),
        mu(N+1, 0), phi(N+1, 0), d(N+1, 0), sigma(N+1, 0),
        is_composite(N+1, false),
        pe(N+1, 1), exp_e(N+1, 0), sigma_pe(N+1, 1)
    {
        // База: n=1
        mu[1] = 1;
        phi[1] = 1;
        d[1] = 1;
        sigma[1] = 1;
        // pe[1]=1, exp_e[1]=0, sigma_pe[1]=1 (нет простых делителей)

        for (int i = 2; i <= N; ++i) {
            if (!is_composite[i]) {
                // i — простое
                primes.push_back(i);
                mu[i] = -1;
                phi[i] = i - 1;
                d[i] = 2;
                sigma[i] = i + 1;
                pe[i] = i;
                exp_e[i] = 1;
                sigma_pe[i] = i + 1;  // σ(p^1) = p + 1
            }

            for (int p : primes) {
                // Не выходим за границу
                if ((long long)i * p > N) break;
                const int ip = i * p;
                is_composite[ip] = true;

                if (i % p == 0) {
                    // p уже присутствует в i с показателем e = exp_e[i]
                    // В ip показатель станет e+1
                    const int e_new = exp_e[i] + 1;
                    exp_e[ip] = e_new;
                    pe[ip] = pe[i] * p;  // p^{e+1}

                    // μ: если p² | ip, то μ(ip) = 0
                    mu[ip] = 0;

                    // φ: φ(p^{e+1}·m) = p · φ(p^e · m) = p · φ(i)
                    phi[ip] = phi[i] * p;

                    // d: d(p^{e+1}·m) = (e+2)·d(m) = d(p^e·m)/(e+1)·(e+2)
                    // exp_e[i] = e, e+1 = exp_e[ip]-1... нет
                    // d[i] = (e+1) · d(m), d[ip] = (e+2) · d(m)
                    // d[ip] = d[i] / (e+1) * (e+2) = d[i] / exp_e[i] * (exp_e[i]+1)
                    // но exp_e[i] — старый показатель (e), exp_e[ip] = e+1
                    // Проверим: exp_e[i]=e (показатель p в i)
                    //   d[i] = (e+1)*d(m), нужно d[ip]=(e+2)*d(m)
                    //   d[ip] = d[i] / (e+1) * (e+2) = d[i] / exp_e[i] * (exp_e[i]+1)
                    d[ip] = d[i] / exp_e[i] * (exp_e[i] + 1);

                    // σ: σ(p^{e+1}) = σ(p^e) + p^{e+1}
                    //    σ(ip) = σ(p^{e+1}) · σ(m) = (σ_pe[i] + pe[ip]) · σ(m)
                    //    σ(i)  = σ(p^e) · σ(m) = σ_pe[i] · σ(m)
                    //    σ(ip) = σ(i) / σ_pe[i] · (σ_pe[i] + pe[ip])
                    const long long new_sigma_pe = sigma_pe[i] + pe[ip];
                    sigma_pe[ip] = new_sigma_pe;
                    sigma[ip] = sigma[i] / sigma_pe[i] * new_sigma_pe;

                    break;  // Важнейшая строка! Прерываем цикл по простым.
                } else {
                    // p — новый простой делитель ip, не встречавшийся в i
                    // ip = i · p, НОД(i, p) = 1
                    exp_e[ip] = 1;
                    pe[ip] = p;

                    mu[ip] = -mu[i];
                    phi[ip] = phi[i] * (p - 1);
                    d[ip] = d[i] * 2;
                    sigma_pe[ip] = p + 1;
                    sigma[ip] = sigma[i] * (p + 1);
                }
            }
        }
    }

    // Вычисление произвольной мультипликативной функции f
    // f_prime[p] = f(p) для простых p
    // f_pk — функция, вычисляющая f(p^{k+1}) из f(p^k): принимает (p, f_pk, p^k)
    vector<long long> compute(
        const function<long long(int)>& f_prime,
        const function<long long(int, long long, int)>& f_next_pk
        // f_next_pk(p, f_pk, p^k) = f(p^{k+1})
    ) const {
        vector<long long> f(N+1, 0);
        vector<long long> f_pe(N+1, 1); // f(p^e) — значение на p-компоненте
        f[1] = 1;

        for (int i = 2; i <= N; ++i) {
            if (!is_composite[i]) {
                f[i] = f_prime(i);
                f_pe[i] = f[i];
            }
            for (int p : primes) {
                if ((long long)i * p > N) break;
                const int ip = i * p;
                if (i % p == 0) {
                    // f(ip) = f(i) / f_pe[i] * f(p^{e+1})
                    // где f(p^{e+1}) = f_next_pk(p, f_pe[i], pe[i])
                    long long new_fpe = f_next_pk(p, f_pe[i], pe[i]);
                    f_pe[ip] = new_fpe;
                    f[ip] = f[i] / f_pe[i] * new_fpe;
                    break;
                } else {
                    f_pe[ip] = f_prime(p);
                    f[ip] = f[i] * f_prime(p);
                }
            }
        }
        return f;
    }
};

// Тест и демонстрация
void run_tests(int N) {
    LinearSieve sieve(N);

    // Проверка μ: sum_{d|n} μ(d) = [n=1]
    bool mu_ok = true;
    for (int n = 1; n <= min(N, 1000); ++n) {
        int s = 0;
        for (int d = 1; d * d <= n; ++d) {
            if (n % d == 0) {
                s += sieve.mu[d];
                if (d != n/d) s += sieve.mu[n/d];
            }
        }
        if (s != (n == 1 ? 1 : 0)) { mu_ok = false; break; }
    }
    cout << "μ корректна: " << (mu_ok ? "ДА" : "НЕТ") << "\n";

    // Проверка φ: sum_{d|n} φ(d) = n
    bool phi_ok = true;
    for (int n = 1; n <= min(N, 1000); ++n) {
        long long s = 0;
        for (int d = 1; d * d <= n; ++d) {
            if (n % d == 0) {
                s += sieve.phi[d];
                if (d != n/d) s += sieve.phi[n/d];
            }
        }
        if (s != n) { phi_ok = false; break; }
    }
    cout << "φ корректна: " << (phi_ok ? "ДА" : "НЕТ") << "\n";

    // Проверка d: d(n) = #{делителей n}
    bool d_ok = true;
    for (int n = 1; n <= min(N, 1000); ++n) {
        int cnt = 0;
        for (int d = 1; d * d <= n; ++d) {
            if (n % d == 0) { cnt++; if (d != n/d) cnt++; }
        }
        if (cnt != sieve.d[n]) { d_ok = false; break; }
    }
    cout << "d корректна: " << (d_ok ? "ДА" : "НЕТ") << "\n";

    // Проверка σ: σ(n) = sum делителей n
    bool sigma_ok = true;
    for (int n = 1; n <= min(N, 1000); ++n) {
        long long s = 0;
        for (int d = 1; d * d <= n; ++d) {
            if (n % d == 0) { s += d; if (d != n/d) s += n/d; }
        }
        if (s != sieve.sigma[n]) { sigma_ok = false; break; }
    }
    cout << "σ корректна: " << (sigma_ok ? "ДА" : "НЕТ") << "\n";
}

// Вычисление σ_k(n) = sum_{d|n} d^k через линейное решето
// Используем обобщённую функцию compute
vector<long long> compute_sigma_k(const LinearSieve& sieve, int k) {
    // f(p) = p^k + 1 (σ_k(p) = 1^k + p^k = 1 + p^k)
    // f(p^e) = 1 + p^k + p^{2k} + ... + p^{ek} = (p^{k(e+1)} - 1)/(p^k - 1)
    // f(p^{e+1}) = f(p^e) + p^{k(e+1)}
    // Нам нужно знать p^{ke}, хранить в f_pe
    return sieve.compute(
        [k](int p) -> long long {
            // f(p) = 1 + p^k
            long long pk = 1;
            for (int i = 0; i < k; ++i) pk *= p;
            return 1 + pk;
        },
        [k](int p, long long fpe, int pe_val) -> long long {
            // fpe = σ_k(p^e) = 1 + p^k + ... + p^{ek}
            // f(p^{e+1}) = fpe + p^{k(e+1)}
            // pe_val = p^e, p^{k(e+1)} = pe_val^k * p^k
            long long pk = 1;
            for (int i = 0; i < k; ++i) pk *= p;
            long long pek = 1;
            long long pv = pe_val;
            for (int i = 0; i < k; ++i) pek *= pv;
            return fpe + pek * pk;
        }
    );
}

// Замер времени: линейное решето vs обычное
#include <chrono>

void benchmark(int N) {
    using namespace chrono;

    // Линейное решето
    auto t1 = high_resolution_clock::now();
    LinearSieve sieve(N);
    auto t2 = high_resolution_clock::now();
    auto linear_time = duration_cast<milliseconds>(t2 - t1).count();

    // Обычное O(N log N) вычисление φ
    t1 = high_resolution_clock::now();
    vector<long long> phi_slow(N+1, 0);
    for (int d = 1; d <= N; ++d)
        for (int n = d; n <= N; n += d)
            phi_slow[n] += sieve.mu[d] * (n/d);
    t2 = high_resolution_clock::now();
    auto slow_time = duration_cast<milliseconds>(t2 - t1).count();

    cout << "N=" << N << "\n";
    cout << "  Линейное решето (mu,phi,d,sigma): " << linear_time << " ms\n";
    cout << "  O(N log N) через μ: " << slow_time << " ms\n";
}

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

    const int N = 1'000'000;

    cout << "=== Тесты корректности (N=" << N << ") ===\n";
    run_tests(N);

    cout << "\n=== Таблица значений (n=1..20) ===\n";
    LinearSieve small(20);
    cout << setw(4) << "n"
         << setw(5) << "mu"
         << setw(5) << "phi"
         << setw(5) << "d"
         << setw(8) << "sigma"
         << setw(6) << "pe"
         << setw(6) << "exp_e" << "\n";
    for (int i = 1; i <= 20; ++i) {
        cout << setw(4) << i
             << setw(5) << small.mu[i]
             << setw(5) << small.phi[i]
             << setw(5) << small.d[i]
             << setw(8) << small.sigma[i]
             << setw(6) << small.pe[i]
             << setw(6) << small.exp_e[i] << "\n";
    }

    cout << "\n=== Бенчмарк ===\n";
    benchmark(N);

    return 0;
}

Вычисление нескольких функций одновременно

Ключевое преимущество линейного решета — возможность вычислять μ, φ, d, σ и другие функции в одном проходе за O(N), без увеличения суммарной сложности. Достаточно добавить соответствующие строки в тело решета.

// Минимальная версия: вычисление mu, phi, d за один проход
vector<int> mu(N+1, 0), phi(N+1, 0), d(N+1, 0), exp_e(N+1, 0);
vector<int> primes;
vector<bool> is_composite(N+1, false);

mu[1] = phi[1] = d[1] = 1;

for (int i = 2; i <= N; ++i) {
    if (!is_composite[i]) {
        primes.push_back(i);
        mu[i] = -1; phi[i] = i-1; d[i] = 2; exp_e[i] = 1;
    }
    for (int p : primes) {
        if ((long long)i*p > N) break;
        int ip = i*p;
        is_composite[ip] = true;
        if (i % p == 0) {
            mu[ip] = 0;
            phi[ip] = phi[i] * p;
            exp_e[ip] = exp_e[i] + 1;
            d[ip] = d[i] / exp_e[i] * exp_e[ip];
            break;
        }
        mu[ip] = -mu[i];
        phi[ip] = phi[i] * (p-1);
        exp_e[ip] = 1;
        d[ip] = d[i] * 2;
    }
}

Разбор задачи 1 (простая): Сумма φ(n) для n от 1 до N

Задача. Дано N (1 ≤ N ≤ 10⁷). Найдите Σ_{n=1}^{N} φ(n) по модулю 10⁹+7.

Анализ. Тривиальная задача, демонстрирующая применение линейного решета.

Из тождества id = φ * 1 следует, что Σ_{n=1}^{N} φ(n) = (N(N+1)/2 + [N≥2]·... ) — нет, это другое. Просто сумма φ.

Есть формула через число взаимно простых пар: Σ_{n=1}^{N} φ(n) = (количество пар (a,n) с 1≤a≤n≤N и НОД(a,n)=1) = (2·Π(N) - 1), где Π(N) — количество взаимно простых пар (a,b) с 1≤a<b≤N. Но проще просто посчитать линейным решетом.

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

const int MOD = 1e9 + 7;

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

    int N;
    cin >> N;

    vector<int> phi(N+1, 0), exp_e(N+1, 0);
    vector<int> primes;
    vector<bool> is_composite(N+1, false);
    phi[1] = 1;

    for (int i = 2; i <= N; ++i) {
        if (!is_composite[i]) {
            primes.push_back(i);
            phi[i] = i - 1;
            exp_e[i] = 1;
        }
        for (int p : primes) {
            if ((long long)i * p > N) break;
            int ip = i * p;
            is_composite[ip] = true;
            if (i % p == 0) {
                phi[ip] = phi[i] * p;
                exp_e[ip] = exp_e[i] + 1;
                break;
            }
            phi[ip] = phi[i] * (p - 1);
            exp_e[ip] = 1;
        }
    }

    long long ans = 0;
    for (int n = 1; n <= N; ++n)
        ans = (ans + phi[n]) % MOD;

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

Сложность: O(N) — существенно быстрее, чем O(N log N) через суммирование.

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

Задача. Дано N (1 ≤ N ≤ 10⁶) и k (0 ≤ k ≤ 20). Для каждого числа от 1 до N определить ω(n) — количество различных простых делителей. Найдите Σ_{n=1}^{N} [ω(n) = k].

Анализ. ω(n) = Σ_{p|n} 1 — количество различных простых делителей. Это аддитивная (не мультипликативная!) функция. Но её можно вычислить через линейное решето.

Через свёртку: ω = Σ_{n: p|n} 1 для каждого простого p, что суммируется за O(N log log N) наивно.

Через линейное решето: при обработке числа ip:

• Если i % p == 0: ω(ip) = ω(i) (p уже учтён)
• Если i % p != 0: ω(ip) = ω(i) + 1 (p — новый простой)

Это ровно O(N) операций!

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

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

    int N, K;
    cin >> N >> K;

    vector<int> omega(N+1, 0);  // ω(n) — число различных простых делителей
    vector<int> primes;
    vector<bool> is_composite(N+1, false);

    // omega[1] = 0 (нет простых делителей)
    for (int i = 2; i <= N; ++i) {
        if (!is_composite[i]) {
            primes.push_back(i);
            omega[i] = 1;  // само является простым
        }
        for (int p : primes) {
            if ((long long)i * p > N) break;
            int ip = i * p;
            is_composite[ip] = true;
            if (i % p == 0) {
                omega[ip] = omega[i];      // p уже в i
                break;
            }
            omega[ip] = omega[i] + 1;     // p — новый простой
        }
    }

    // Попутно вычислим Ω(n) — число простых с учётом кратности
    vector<int> Omega(N+1, 0);
    for (int i = 2; i <= N; ++i) {
        if (!is_composite[i]) {
            // Omega уже 0, обновим при пометке через прямой цикл
        }
    }
    // Для Ω проще через отдельный цикл:
    for (int p = 2; p <= N; ++p) {
        if (!is_composite[p]) {
            for (long long pk = p; pk <= N; pk *= p)
                for (long long n = pk; n <= N; n += pk)
                    Omega[n]++;
        }
    }

    long long ans = 0;
    for (int n = 1; n <= N; ++n)
        if (omega[n] == K) ans++;

    cout << ans << "\n";

    // Бонус: гистограмма распределения ω(n)
    cout << "\nРаспределение ω(n) для n в [1," << N << "]:\n";
    vector<int> hist(21, 0);
    for (int n = 1; n <= N; ++n)
        if (omega[n] <= 20) hist[omega[n]]++;
    for (int k = 0; k <= 7; ++k)
        cout << "ω(n)=" << k << ": " << hist[k] << " чисел\n";

    return 0;
}

Замечание. Через линейное решето мы вычислили ω(n) для всех n за O(N). Отдельное вычисление Ω(n) потребовало O(N log log N). Можно совместить и Ω в одном решете:

• При i % p == 0: Ω[ip] = Ω[i] + 1 (p повторяется)
• При i % p != 0: Ω[ip] = Ω[i] + 1 (p новый, но в обоих случаях +1!)

Итого: Ω[ip] = Ω[i] + 1 в обоих случаях. Но тогда решето теряет смысл — Ω[ip] = Ω[i] + 1, что можно вычислить и при пометке. Функция Лиувилля λ(n) = (−1)^{Ω(n)} тогда: λ[ip] = −λ[i] всегда.

Разбор задачи 3 (сложная: ICPC/CF Div.1 D уровень): Минимальное произведение с НОД=1

Задача. Дан массив a₁, a₂, ..., a_m (1 ≤ aᵢ ≤ N ≤ 10⁶, m ≤ 10⁶). Найдите количество пар индексов (i, j) с i < j таких, что НОД(aᵢ, aⱼ) = 1. Ответ по модулю 10⁹+7.

Анализ. Обозначим cnt[v] = количество элементов массива, равных v. Тогда:

Ответ = Σ_{v1 < v2} cnt[v1] · cnt[v2] · [НОД(v1,v2)=1]
       + Σ_{v} C(cnt[v], 2) · [НОД(v,v)=1]
       = (общая формула для пар)

Пусть f(d) = (Σ_{d|v} cnt[v])² = (число элементов, делящихся на d)². Тогда:

Σ_{d} μ(d) · f(d) = Σ_{d} μ(d) · (Σ_{d|a} cnt[a])²

Это даёт сумму cnt[a]·cnt[b] по всем парам (a,b) с d|НОД(a,b), взвешенно на μ(d). По формуле Мёбиуса это сумма cnt[a]·cnt[b] по парам с НОД(a,b)=1.

Но нам нужны пары, а не упорядоченные пары. Аккуратно:

Г = Σ_{d=1}^{N} μ(d) · (Σ_{k=1}^{⌊N/d⌋} cnt[dk])²

Это Σ_{i,j} cnt[i]·cnt[j]·[НОД(i,j)=1] = 2·(ответ) + (количество v с cnt[v]>0 и v=1).

Алгоритм:

1. Вычислить cnt[v] за O(m).
2. Вычислить g[d] = Σ_{k} cnt[dk] за O(N log N) (суммирование по делителям).
3. Вычислить Г = Σ_d μ(d)·g[d]² за O(N).
4. Ответ = (Г - [1 в массиве]) / 2.

Ключевое ускорение с линейным решетом: μ нужно вычислить за O(N), иначе будет TLE при N = 10⁶ из-за медленной факторизации.

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

const long long MOD = 1e9 + 7;

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

    int N, m;
    cin >> N >> m;

    // Линейное решето для μ
    vector<int> mu(N+1, 0);
    vector<int> primes;
    vector<bool> is_composite(N+1, false);
    mu[1] = 1;
    for (int i = 2; i <= N; ++i) {
        if (!is_composite[i]) {
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int p : primes) {
            if ((long long)i*p > N) break;
            is_composite[i*p] = true;
            if (i % p == 0) { mu[i*p] = 0; break; }
            mu[i*p] = -mu[i];
        }
    }

    // Подсчёт cnt[v]
    vector<long long> cnt(N+1, 0);
    for (int i = 0; i < m; ++i) {
        int x; cin >> x;
        cnt[x]++;
    }

    // g[d] = sum_{k: dk<=N} cnt[dk] = число элементов, делящихся на d
    vector<long long> g(N+1, 0);
    for (int d = 1; d <= N; ++d) {
        for (int k = d; k <= N; k += d) {
            g[d] += cnt[k];
        }
    }

    // Gamma = sum_d mu[d] * g[d]^2
    // = sum_{i,j} cnt[i]*cnt[j]*[gcd(i,j)=1] (упорядоченные)
    long long Gamma = 0;
    for (int d = 1; d <= N; ++d) {
        if (mu[d] == 0 || g[d] == 0) continue;
        long long gd = g[d] % MOD;
        long long term = mu[d] * (gd * gd % MOD) % MOD;
        Gamma = (Gamma + term + MOD) % MOD;
    }

    // Gamma = 2*answer + [НОД(v,v)=1 для каких-то v, т.е. v=1]
    // Вклад пары (v,v): cnt[v]^2 включает cnt[v]^2 = cnt[v]*(cnt[v]-1) + cnt[v]
    // Упорядоченные пары (i,j) с i!=j: Gamma - cnt[1]^2 (вклад пары (1,1)??)
    // Нет: Gamma = sum_{i,j} cnt[i]*cnt[j]*[gcd=1] включает i=j
    // При i=j: НОД(i,i)=i=1 только при i=1, вклад = cnt[1]^2
    // Ответ (пары i<j) = (Gamma - cnt[1]) / 2
    // (cnt[1] = число пар (a,a) где a=1, т.е. cnt[1]^2 упорядоченных включает
    //  cnt[1]*(cnt[1]-1) пар (разные индексы) + cnt[1] пар (i=i))
    // Ответ = (Gamma - cnt[1]) / 2

    long long cnt1_mod = cnt[1] % MOD;
    long long ans = (Gamma - cnt1_mod + MOD) % MOD;
    // Деление на 2 по модулю: умножить на обратный элемент
    // 2^{-1} mod (10^9+7) = (10^9+8)/2 = 500000004
    const long long inv2 = 500000004LL;
    ans = ans * inv2 % MOD;

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

Сложность:

• Решето μ: O(N)
• Подсчёт g: O(N log N)
• Подсчёт Γ: O(N)
• Итого: O(N log N)

Без линейного решета (используя факторизацию каждого числа для μ) шаг вычисления μ занял бы O(N√N) — недопустимо медленно.

Важность линейного решета: Для N = 10⁶:

• O(N√N) = 10¹² операций — невозможно
• O(N log log N) решето + O(N log N) суммирование = ~5·10⁷ — быстро
• O(N) линейное + O(N log N) — примерно так же, но линейное решето даёт μ без overhead

Проверка: Если массив = [1, 2, 3], N=3:

• cnt[1]=1, cnt[2]=1, cnt[3]=1
• g[1]=3, g[2]=1, g[3]=1
• Γ = μ(1)·9 + μ(2)·1 + μ(3)·1 = 9 − 1 − 1 = 7
• ans = (7 − 1) / 2 = 3

Пары (i<j): (1,2)→НОД=1✓, (1,3)→НОД=1✓, (2,3)→НОД=1✓. Итого 3. ✓

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

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

1. Неверный порядок вычислений в решете. Значения pe[ip], exp_e[ip], sigma_pe[ip] должны быть обновлены до использования в формулах для d[ip] и sigma[ip]. Следите за порядком присваиваний.

// НЕВЕРНО:
d[ip] = d[i] / exp_e[i] * (exp_e[i] + 1);
exp_e[ip] = exp_e[i] + 1;  // обновили после использования — OK в этом случае, но опасно

// ВЕРНО и явно:
const int e_old = exp_e[i];       // сохраняем старый показатель
const int e_new = e_old + 1;
exp_e[ip] = e_new;
d[ip] = d[i] / e_old * e_new;

2. Пропущенный break после i % p == 0. Это самая критичная ошибка в линейном решете. Без break некоторые числа будут помечены несколько раз, и значения функций будут неверны.

if (i % p == 0) {
    // ... обновляем функцию ...
    break;  // ОБЯЗАТЕЛЬНО!
}

3. Переполнение pe[ip] = pe[i] * p. При больших N (10⁷ и выше) и простых p ~ √N значение pe[i] может достигать p^k ~ 10^{3.5}. Обычно не проблема при N ≤ 10⁷, но при N = 10⁸ нужен long long.

4. Деление на ноль в sigma[ip] / sigma_pe[i]. При i=1: sigma_pe[1] = 1, деление безопасно. Но если sigma_pe[i] = 0 по ошибке инициализации — получим деление на ноль. Инициализируйте sigma_pe все значения единицей.

5. Неверная инициализация exp_e[1]. exp_e[1] = 0 (у 1 нет простых делителей). Если случайно обратиться к этому значению в формуле d[ip] / exp_e[i] при i=1, получим деление на ноль. Но i=1 не попадёт в цикл (он начинается с i=2), поэтому проблем нет — если только вы не используете exp_e[1] явно.

6. Неверное обобщение на полностью мультипликативные функции. Для полностью мультипликативной f (например, id или 1): f(ip) = f(i)·f(p) ВСЕГДА, без различия случаев. Если написать if (i % p == 0) { f[ip] = f[i] * f[p]; break; }, то получите break и не посетите другие простые. Для полностью мультипликативных функций код корректен, но надо убедиться, что f[ip] = f[i] * f[p] верно и при p|i.

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

На практике выбирать между решетом Эратосфена и линейным нужно с умом:

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

• Для одной мультипликативной функции (например, только φ): можно использовать решето Эратосфена с модификацией. Для φ: инициализировать phi[i]=i, затем для каждого простого p: для всех кратных n умножить phi[n] на (p-1)/p. Это O(N log log N) с хорошей кэш-локальностью.

• Для нескольких функций одновременно или при N ≥ 10⁷: линейное решето.

// Профессиональный шаблон: линейное решето в 30 строк
// вычисляет mu, phi, d, sigma — все сразу
void linear_sieve(int N,
    vector<int>& mu, vector<int>& phi, vector<int>& d, vector<long long>& sigma)
{
    mu.assign(N+1, 0); phi.assign(N+1, 0);
    d.assign(N+1, 0); sigma.assign(N+1, 0);
    vector<int> primes, pe(N+1, 1), ee(N+1, 0);
    vector<long long> spe(N+1, 1);
    vector<bool> comp(N+1, false);

    mu[1]=phi[1]=d[1]=sigma[1]=1;
    for (int i = 2; i <= N; ++i) {
        if (!comp[i]) {
            primes.push_back(i);
            mu[i]=-1; phi[i]=i-1; d[i]=2; sigma[i]=i+1;
            pe[i]=i; ee[i]=1; spe[i]=i+1;
        }
        for (int p : primes) {
            if ((long long)i*p > N) break;
            int ip=i*p; comp[ip]=true;
            if (i%p==0) {
                mu[ip]=0; phi[ip]=phi[i]*p;
                ee[ip]=ee[i]+1; pe[ip]=pe[i]*p;
                long long ns=spe[i]+pe[ip]; spe[ip]=ns;
                d[ip]=d[i]/ee[i]*ee[ip];
                sigma[ip]=sigma[i]/spe[i]*ns;
                break;
            }
            mu[ip]=-mu[i]; phi[ip]=phi[i]*(p-1);
            ee[ip]=1; pe[ip]=p; spe[ip]=p+1;
            d[ip]=d[i]*2; sigma[ip]=sigma[i]*(p+1);
        }
    }
}

Второй совет: при реализации задач на соревнованиях, если N ≤ 3·10⁶, линейное решето должно работать в пределах 1 секунды. Для N = 10⁷ может потребоваться оптимизация: избегайте vector<bool> (используйте vector<char>) и старайтесь минимизировать промахи кэша, группируя вычисления.

Итог

Линейное решето — венец изучения базового числового решета. Основные выводы:

1. Идея: Каждое составное число помечается ровно один раз — как наименьший простой делитель × некоторое число. Это достигается ключевым break при p|i.
2. Сложность O(N): Строго линейная, каждое составное посещается один раз. На практике быстрее O(N log log N) при больших N.
3. Мультипликативные функции: Два случая в решете — p|i и p∤i — соответствуют двум правилам мультипликативности: рекурсия по p-компоненте и f(ip)=f(i)f(p).
4. Несколько функций одновременно: μ, φ, d, σ вычисляются в одном проходе без дополнительных накладных расходов.
5. pe, exp_e, sigma_pe: Вспомогательные массивы для корректного вычисления d и σ. Их правильное обновление — главная сложность реализации.
6. Практическое применение: Линейное решето необходимо при N ≥ 10⁷ или когда нужно несколько мультипликативных функций. При N ≤ 10⁶ можно использовать более простое решето Эратосфена.