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

В предыдущей главе мы изучили функцию Мёбиуса и увидели, что она обладает особым свойством: μ(mn) = μ(m)·μ(n) при НОД(m,n) = 1. Это свойство называется мультипликативностью, и оно присуще целому классу важнейших функций теории чисел. Понимание мультипликативности — ключ к эффективным алгоритмам: зная значения функции на простых степенях, можно вычислить её на любом числе за O(log n).

Свёртка Дирихле — алгебраическая операция над арифметическими функциями — позволяет выражать одни функции через другие и формально обосновывает все «магические» тождества, которые мы видели в главе про Мёбиуса. Например:

• φ = μ * id (функция Эйлера = свёртка Мёбиуса с тождественной)
• d = 1 * 1 (число делителей = свёртка единицы с собой)
• σ = id * 1 (сумма делителей = свёртка тождественной с единицей)

В соревновательном программировании задачи на свёртку Дирихле встречаются на уровне Div.1 C–E. Понимание этой темы открывает доступ к задачам, которые без неё кажутся неподступными.

Теория

Арифметические функции

Определение. Арифметическая функция — это функция f : ℕ → ℂ (или ℤ, или ℝ — в зависимости от контекста).

Мультипликативные функции

Определение. Функция f называется мультипликативной, если:

1. f не тождественно равна нулю
2. Для любых m, n с НОД(m, n) = 1 выполняется f(mn) = f(m)·f(n)

Определение. Функция f называется полностью мультипликативной, если f(mn) = f(m)·f(n) для всех m, n (без условия взаимной простоты).

Следствие для мультипликативных функций. Если f мультипликативна, то f(1) = 1. Доказательство: f(1) = f(1·1) = f(1)·f(1), значит f(1) ∈ {0, 1}. Поскольку f не тождественно 0, и f(n) = f(1·n) = f(1)·f(n), получаем f(1) = 1.

Важнейшее следствие. Если f мультипликативна и n = p₁^{a₁}·...·p_k^{a_k}, то:

f(n) = f(p₁^{a₁}) · f(p₂^{a₂}) · ... · f(p_k^{a_k})

Это означает, что мультипликативная функция полностью определяется своими значениями на степенях простых.

Стандартные арифметические функции

Рассмотрим важнейшие примеры:

1. Тождественная единица: 1(n) = 1 для всех n.

• Полностью мультипликативна: 1(mn) = 1 = 1(m)·1(n). 2. Тождественная функция: id(n) = n.
• Полностью мультипликативна: id(mn) = mn = id(m)·id(n). 3. Нейтральный элемент: ε(n) = [n = 1].
• Мультипликативна (и полностью мультипликативна, так как ε(mn) = [mn=1] = [m=1][n=1] = ε(m)ε(n) при любых m,n ≥ 1). 4. Функция Мёбиуса: μ(n).
• Мультипликативна (доказано в главе 7). 5. Функция Эйлера: φ(n) = #{k ≤ n : НОД(k,n) = 1}.
• Мультипликативна. Доказательство через КТО: если НОД(m,n)=1, то ℤ/mn ≅ ℤ/m × ℤ/n, поэтому φ(mn) = φ(m)·φ(n).
• φ(p^k) = p^k − p^{k−1} = p^{k−1}(p−1). 6. Функция числа делителей: d(n) = Σ_{d|n} 1 = τ(n).
• d(p^k) = k + 1.
• Мультипликативна: при НОД(m,n)=1 каждый делитель mn однозначно представляется как d₁·d₂ где d₁|m, d₂|n, значит d(mn) = d(m)·d(n).

7. Функция суммы делителей: σ(n) = Σ_{d|n} d.

• σ(p^k) = 1 + p + p² + ... + p^k = (p^{k+1}−1)/(p−1).
• Мультипликативна: аналогично d(n), Σ_{d|mn} d = (Σ_{d₁|m} d₁)·(Σ_{d₂|n} d₂) = σ(m)·σ(n) при НОД(m,n)=1. 8. Обобщённые функции σ_k(n) = Σ_{d|n} d^k.
• σ₀ = d (число делителей), σ₁ = σ (сумма делителей).
• Все мультипликативны. 9. Функция Лиувилля: λ(n) = (−1)^{Ω(n)}, где Ω(n) = число простых делителей с учётом кратности.
• Полностью мультипликативна: λ(mn) = (−1)^{Ω(mn)} = (−1)^{Ω(m)+Ω(n)} = λ(m)·λ(n).
• λ(p^k) = (−1)^k.

Таблица значений основных функций:

Свёртка Дирихле

Определение. Свёртка Дирихле двух арифметических функций f и g — это функция:

(f * g)(n) = Σ_{d|n} f(d) · g(n/d)

Теорема об алгебраических свойствах. Свёртка Дирихле обладает следующими свойствами:

(a) Коммутативность: f * g = g * f.

Доказательство: (f * g)(n) = Σ_{d|n} f(d)g(n/d) = Σ_{d|n} f(n/d)g(d) = (g * f)(n), так как при замене d → n/d пробегаем те же делители.

(b) Ассоциативность: (f * g) * h = f * (g * h).

Доказательство: ((f * g) * h)(n) = Σ_{d|n} (f * g)(d)·h(n/d) = Σ_{d|n} Σ_{e|d} f(e)·g(d/e)·h(n/d).

Делая замену: пусть a = e, b = d/e, c = n/d (тогда abc = n): = Σ_{abc=n} f(a)·g(b)·h(c).

Аналогично (f * (g * h))(n) = Σ_{abc=n} f(a)·g(b)·h(c).

(c) Нейтральный элемент: f * ε = ε * f = f.

Доказательство: (f * ε)(n) = Σ_{d|n} f(d)·ε(n/d) = f(n)·ε(1) = f(n).

(d) Обратный элемент. Для любой f с f(1) ≠ 0 существует g такое, что f * g = ε. Функция g называется обратной по Дирихле и обозначается f^{−1}.

Доказательство существования: Определим g рекурсивно. При n=1: f(1)·g(1) = 1, значит g(1) = 1/f(1). При n>1:

0 = Σ_{d|n} f(d)·g(n/d) = f(1)·g(n) + Σ_{d|n, d<n} f(d)·g(n/d)

Отсюда g(n) = −(1/f(1)) · Σ_{d|n, d<n} f(d)·g(n/d).

Это корректная рекурсия, поскольку g(n) выражается только через g(m) с m < n.

Следствие. Арифметические функции с f(1) ≠ 0 образуют абелеву группу по свёртке Дирихле. Все арифметические функции образуют коммутативное кольцо (с поточечным сложением и свёрткой как умножением).

Мультипликативность свёртки

Теорема. Если f и g — мультипликативные функции, то f * g тоже мультипликативна.

Доказательство. Пусть НОД(m,n) = 1. Тогда:

(f * g)(mn) = Σ_{d|mn} f(d)·g(mn/d)

Ключевое наблюдение: при НОД(m,n)=1 каждый делитель d числа mn однозначно представляется как d = d₁·d₂, где d₁|m и d₂|n (и при этом НОД(d₁,d₂)=1). Поэтому:

Σ_{d|mn} f(d)·g(mn/d) = Σ_{d₁|m} Σ_{d₂|n} f(d₁d₂)·g((mn)/(d₁d₂))
                       = Σ_{d₁|m} Σ_{d₂|n} f(d₁)f(d₂)·g(m/d₁)g(n/d₂)
                       = (Σ_{d₁|m} f(d₁)g(m/d₁)) · (Σ_{d₂|n} f(d₂)g(n/d₂))
                       = (f*g)(m) · (f*g)(n)

Здесь использовали мультипликативность f и g: f(d₁d₂)=f(d₁)f(d₂) и g(m/d₁·n/d₂)=g(m/d₁)g(n/d₂).

Следствие. Обратная по Дирихле к мультипликативной функции тоже мультипликативна.

Важнейшие тождества (таблица свёрток)

Дзета-функция Дирихле (дискретная версия)

Аналог непрерывной дзета-функции Римана ζ(s) = Σ 1/n^s в дискретном случае — это функция 1 (тождественно равная 1). Умножение на 1 в свёртке Дирихле:

(f * 1)(n) = Σ_{d|n} f(d)

«суммирует» f по делителям. В аналитической теории чисел:

• Функция ζ(s) = Σ_{n≥1} 1/n^s = произведение по простым (1−1/p^s)^{−1}
• 1/ζ(s) = Σ_{n≥1} μ(n)/n^s — ряд Дирихле для μ
• ζ(s)² = Σ_{n≥1} d(n)/n^s — ряд для функции числа делителей

Дискретный аналог этих равенств:

• 1 * 1 = d (суммирование 1 по делителям даёт d(n))
• μ * 1 = ε (обращение)
• 1 * 1 = d

Формула для обратного элемента по Дирихле

Для мультипликативной f её обратная f^{−1} вычисляется через значения на степенях простых.

При n = p^k:

Σ_{j=0}^{k} f(p^j) · f^{-1}(p^{k-j}) = [k=0]

Отсюда рекурсия: f^{-1}(1) = 1, и для k≥1:

f^{-1}(p^k) = −(1/f(1)) · Σ_{j=1}^{k} f(p^j) · f^{-1}(p^{k-j})
             = −Σ_{j=1}^{k} f(p^j) · f^{-1}(p^{k-j})  (так как f(1)=1)

Пример: Найдём 1^{-1} (обратную к тождественной единице).

• 1^{-1}(1) = 1
• 1^{-1}(p) = −1 (так как f(p)·f^{-1}(1) + f(1)·f^{-1}(p) = 0 → f^{-1}(p) = −1)
• 1^{-1}(p²) = −1^{-1}(p)·1(p) − 1^{-1}(1)·1(p²) = 1 − 1 = 0
• 1^{-1}(p^k) = 0 для k ≥ 2

Итого: 1^{-1}(n) = μ(n). Что совпадает с μ * 1 = ε! ✓

Вычисление свёртки за O(N log N)

Для массивов f[1..N] и g[1..N] вычислить h = f * g за O(N log N):

for d = 1 to N:
    for k = 1 to N/d:
        h[d*k] += f[d] * g[k]

Количество операций: Σ_{d=1}^{N} N/d ≈ N·ln N = O(N log N).

Более эффективный способ через быстрое преобразование Дирихле (аналог БПФ, но в базисе простых), но это выходит за рамки главы.

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

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

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

// Вычисление мультипликативных функций через линейное решето
// Возвращает: mu, phi, d (число делителей), sigma (сумма делителей)
struct MultiplicativeFunctions {
    int N;
    vector<int> mu, phi, d, smallest_prime, prime_power_exponent;
    vector<long long> sigma;
    vector<int> primes;
    vector<bool> is_composite;

    explicit MultiplicativeFunctions(int N) : N(N),
        mu(N+1), phi(N+1), d(N+1), sigma(N+1),
        smallest_prime(N+1, 0), prime_power_exponent(N+1, 0),
        is_composite(N+1, false)
    {
        mu[1] = phi[1] = d[1] = sigma[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;
                sigma[i] = i + 1;
                smallest_prime[i] = i;
                prime_power_exponent[i] = 1;  // i = i^1, показатель = 1
            }

            for (int p : primes) {
                if ((long long)i * p > N) break;
                int ip = i * p;
                is_composite[ip] = true;
                smallest_prime[ip] = p;

                if (i % p == 0) {
                    // p делит i, p² делит ip
                    mu[ip] = 0;
                    phi[ip] = phi[i] * p;

                    // Для d и sigma нужно знать показатель p в i
                    int e = prime_power_exponent[i] + 1;
                    prime_power_exponent[ip] = e;

                    // i = A * p^(e-1), ip = A * p^e
                    // d(ip) = d(A) * (e+1) / (e) * (e) = d(A) * (e+1) / e * ...
                    // Проще: ip/p^e = i/p^(e-1) = A (часть без p)
                    // Но для решета аккуратнее:
                    // d(ip) = d(i/p^(e-1)) * (e+1) -- нет, нужна часть без p
                    // Используем хитрость: prime_power[i] = p^(e-1) * что-то
                    // Проще хранить d[i / (p^e / p)] * (e+1)
                    // Здесь ip = i*p, i делится на p с показателем e-1
                    // Пусть A = i / p^(e-1), тогда d(A) * e = d(i)
                    // d(ip) = d(A) * (e+1)
                    // d(A) = d(i) / e
                    d[ip] = d[i] / e * (e + 1);

                    // sigma(p^e) = 1 + p + ... + p^e = sigma(p^(e-1)) + p^e
                    // sigma(ip) = sigma(A) * sigma(p^e)
                    // sigma(i) = sigma(A) * sigma(p^(e-1))
                    // sigma(ip) / sigma(i) = sigma(p^e) / sigma(p^(e-1))
                    // sigma(p^e) = sigma(p^(e-1)) * (p^e + sigma(p^(e-1)) - sigma(p^(e-1)) + p^e)
                    // Аккуратнее: sigma(p^e) = sigma(p^{e-1}) + p^e (геом. прогрессия)
                    // Пусть sp = sigma(p^{e-1}) — нам нужно sigma(p^e)
                    // sigma(p^e) = (p^{e+1} - 1)/(p-1) -- точная формула, но нужна рекурсия
                    // sigma(ip) = (sigma(i) / sigma_pk) * sigma_pk_plus_one
                    // где sigma_pk = sigma(p^{e-1}), sigma_pk_plus_one = sigma(p^e)
                    // Это сложно без дополнительного хранения.
                    // Используем другой подход для sigma:
                    sigma[ip] = sigma[i] + sigma[i] / (/* sigma(p^{e-1}) */ 1) * 1;
                    // (будет заменено в отдельной реализации)
                    break;
                } else {
                    // p не делит i, p — новый простой в ip
                    mu[ip] = -mu[i];
                    phi[ip] = phi[i] * (p - 1);
                    d[ip] = d[i] * 2;
                    sigma[ip] = sigma[i] * (p + 1);
                    prime_power_exponent[ip] = 1;
                }
            }
        }
    }
};

// Аккуратная реализация линейного решета для mu, phi, d
// sigma требует дополнительного хранения sigma_pk
struct LinearSieve {
    int N;
    vector<int> mu, phi, d;
    vector<long long> sigma;
    vector<int> primes;
    vector<bool> is_composite;
    // pe[i] = p^e где p = наименьший простой делитель i, e = его показатель
    // Нужно для правильного вычисления d и sigma
    vector<int> pe;      // p^e (мощность наименьшего простого в n)
    vector<int> exp_pe;  // показатель e наименьшего простого
    vector<long long> sigma_pe; // sigma(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_pe(N+1,0), sigma_pe(N+1,1)
    {
        mu[1] = phi[1] = d[1] = sigma[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;
                sigma[i] = i + 1;
                pe[i] = i;       // p^1 = p
                exp_pe[i] = 1;
                sigma_pe[i] = i + 1;  // sigma(p^1) = p+1
            }

            for (int p : primes) {
                if ((long long)i * p > N) break;
                long long ip = (long long)i * p;
                is_composite[ip] = true;

                if (i % p == 0) {
                    // p | i, p^(e+1) || ip
                    int e = exp_pe[i] + 1;
                    exp_pe[ip] = e;
                    pe[ip] = pe[i] * p;  // p^e

                    // d(ip) = d(i / p^(e-1)) * (e+1)
                    // d(i) = d(i/p^(e-1)) * e
                    // d(ip) = d(i) / e * (e+1)
                    d[ip] = d[i] / (e - 1) * e;

                    // sigma(p^e) = sigma(p^{e-1}) + p^e = sigma_pe[i] + pe[i]*p
                    long long sp = sigma_pe[i] + (long long)pe[ip];
                    sigma_pe[ip] = sp;
                    // sigma(ip) = sigma(i / p^{e-1}) * sigma(p^e)
                    // sigma(i) = sigma(i / p^{e-1}) * sigma(p^{e-1})
                    // sigma(ip) = sigma(i) / sigma_pe[i] * sp
                    sigma[ip] = sigma[i] / sigma_pe[i] * sp;

                    mu[ip] = 0;
                    phi[ip] = phi[i] * p;
                    break;
                } else {
                    // p не делит i
                    exp_pe[ip] = 1;
                    pe[ip] = p;

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

// Свёртка Дирихле двух функций за O(N log N)
vector<long long> dirichlet_convolution(
    const vector<long long>& f,
    const vector<long long>& g)
{
    int N = (int)f.size() - 1;
    vector<long long> h(N + 1, 0);
    for (int d = 1; d <= N; ++d) {
        if (f[d] == 0) continue;
        for (int k = 1; k * d <= N; ++k) {
            h[d * k] += f[d] * g[k];
        }
    }
    return h;
}

// Инверсия Мёбиуса: g = f * mu, т.е. если f = h * 1, то g = h
vector<long long> mobius_inversion(const vector<long long>& f, const vector<int>& mu) {
    int N = (int)f.size() - 1;
    vector<long long> g(N + 1, 0);
    for (int d = 1; d <= N; ++d) {
        if (mu[d] == 0) continue;
        for (int k = 1; (long long)d * k <= N; ++k) {
            g[d * k] += (long long)mu[d] * f[k];
        }
    }
    return g;
}

// Вычисление обратного по Дирихле для мультипликативной функции
// Используем рекурсию: f_inv[n] = -sum_{d|n, d<n} f[d] * f_inv[n/d]
vector<long long> dirichlet_inverse(const vector<long long>& f) {
    int N = (int)f.size() - 1;
    assert(f[1] == 1 || f[1] == -1);  // f(1) должна быть +-1

    vector<long long> inv(N + 1, 0);
    inv[1] = 1;  // предполагаем f[1] = 1

    for (int n = 2; n <= N; ++n) {
        for (int d = 1; d < n; ++d) {
            if (n % d == 0) {
                inv[n] -= f[n / d] * inv[d];
            }
        }
        // inv[n] /= f[1] = 1
    }
    return inv;
}

// Демонстрация тождеств
void verify_identities(int N) {
    LinearSieve sieve(N);
    auto& mu = sieve.mu;
    auto& phi = sieve.phi;
    auto& d_func = sieve.d;
    auto& sigma = sieve.sigma;

    // Преобразуем в long long для свёртки
    vector<long long> one(N+1, 1), id_func(N+1), mu_ll(N+1), phi_ll(N+1);
    for (int i = 1; i <= N; ++i) {
        id_func[i] = i;
        mu_ll[i] = mu[i];
        phi_ll[i] = phi[i];
    }

    // Проверка: 1 * 1 = d
    auto d_check = dirichlet_convolution(one, one);
    bool ok1 = true;
    for (int i = 1; i <= N; ++i)
        if (d_check[i] != d_func[i]) { ok1 = false; break; }
    cout << "1*1 = d: " << (ok1 ? "OK" : "FAIL") << "\n";

    // Проверка: id * 1 = sigma
    auto sigma_check = dirichlet_convolution(id_func, one);
    bool ok2 = true;
    for (int i = 1; i <= N; ++i)
        if (sigma_check[i] != sigma[i]) { ok2 = false; break; }
    cout << "id*1 = sigma: " << (ok2 ? "OK" : "FAIL") << "\n";

    // Проверка: mu * id = phi
    auto phi_check = dirichlet_convolution(mu_ll, id_func);
    bool ok3 = true;
    for (int i = 1; i <= N; ++i)
        if (phi_check[i] != phi_ll[i]) { ok3 = false; break; }
    cout << "mu*id = phi: " << (ok3 ? "OK" : "FAIL") << "\n";

    // Проверка: phi * 1 = id
    auto id_check = dirichlet_convolution(phi_ll, one);
    bool ok4 = true;
    for (int i = 1; i <= N; ++i)
        if (id_check[i] != id_func[i]) { ok4 = false; break; }
    cout << "phi*1 = id: " << (ok4 ? "OK" : "FAIL") << "\n";

    // Проверка: mu * 1 = epsilon
    auto eps_check = dirichlet_convolution(mu_ll, one);
    bool ok5 = (eps_check[1] == 1);
    for (int i = 2; i <= N; ++i)
        if (eps_check[i] != 0) { ok5 = false; break; }
    cout << "mu*1 = eps: " << (ok5 ? "OK" : "FAIL") << "\n";
}

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

    const int N = 1000;

    cout << "Проверка тождеств для N=" << N << ":\n";
    verify_identities(N);

    // Таблица значений
    LinearSieve sieve(20);
    cout << "\nТаблица мультипликативных функций:\n";
    cout << setw(4) << "n"
         << setw(6) << "mu"
         << setw(6) << "phi"
         << setw(6) << "d"
         << setw(8) << "sigma" << "\n";
    for (int i = 1; i <= 20; ++i) {
        cout << setw(4) << i
             << setw(6) << sieve.mu[i]
             << setw(6) << sieve.phi[i]
             << setw(6) << sieve.d[i]
             << setw(8) << sieve.sigma[i] << "\n";
    }

    return 0;
}

Разбор задачи 1 (простая): Количество делителей и сумма делителей

Задача. Даны числа N и M (1 ≤ N, M ≤ 10⁶). Для каждого числа от 1 до N найдите: (a) количество делителей d(n) (b) сумму делителей σ(n) (c) сумму Σ_{n=1}^{N} d(n) · σ(n) по модулю 10⁹+7.

Анализ. d и σ — мультипликативные функции, их можно вычислить:

1. Через линейное решето за O(N).
2. Через суммирование по делителям за O(N log N).

Для пункта (c) после вычисления массивов d[i] и σ[i] ответ считается за O(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> d(N+1, 1), sigma(N+1, 1);
    vector<int> pe(N+1, 1);    // p^e для наименьшего простого делителя
    vector<int> exp_e(N+1, 0); // показатель этого простого
    vector<long long> sigma_pe(N+1, 1); // sigma(p^e)

    vector<int> primes;
    vector<bool> is_composite(N+1, false);

    d[1] = sigma[1] = 1;

    for (int i = 2; i <= N; ++i) {
        if (!is_composite[i]) {
            primes.push_back(i);
            d[i] = 2;
            sigma[i] = i + 1;
            pe[i] = i;
            exp_e[i] = 1;
            sigma_pe[i] = 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) {
                int e = exp_e[i] + 1;
                exp_e[ip] = e;
                pe[ip] = pe[i] * p;
                // d(ip) = d(i/p^{e-1}) * e = d(i) / (e-1) * e
                d[ip] = d[i] / (e-1) * e;
                // sigma(p^e) = sigma(p^{e-1}) + p^e
                long long sp = sigma_pe[i] + (long long)pe[ip];
                sigma_pe[ip] = sp;
                // sigma(ip) = sigma(i) / sigma_pe[i] * sp
                sigma[ip] = (long long)sigma[i] / sigma_pe[i] * sp;
                break;
            } else {
                exp_e[ip] = 1;
                pe[ip] = p;
                d[ip] = d[i] * 2;
                sigma_pe[ip] = p + 1;
                sigma[ip] = (long long)sigma[i] * (p + 1);
            }
        }
    }

    // Подсчёт суммы d(n) * sigma(n) mod MOD
    long long ans = 0;
    for (int n = 1; n <= N; ++n) {
        ans = (ans + (long long)d[n] % MOD * (sigma[n] % MOD)) % MOD;
    }
    cout << ans << "\n";

    return 0;
}

Замечание по реализации. В линейном решете для d и σ ключевой момент — аккуратное отслеживание показателя наименьшего простого. Формулы:

• d(p^e · m) = (e+1) · d(m) при НОД(p,m)=1 → при добавлении нового p: d[ip] = d[i]*2
• d(p^e · m) при p|i: e увеличивается на 1, d[ip] = d[i]/(e_old)*(e_old+1)

Разбор задачи 2 (средняя): Сумма φ(НОД(i,j)) для всех пар

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

Анализ. Воспользуемся тождеством φ * 1 = id (то есть Σ_{d|n} φ(d) = n):

Σ_{i,j} φ(НОД(i,j)) = Σ_{d=1}^{N} φ(d) · #{(i,j): НОД(i,j)=d}

Число пар (i,j) с НОД(i,j)=d: это пары (i',j') с i=d·i', j=d·j', НОД(i',j')=1, и 1≤i',j'≤⌊N/d⌋. Обозначим C(m) = #{(a,b): 1≤a,b≤m, НОД(a,b)=1} = Σ_{k=1}^{m} μ(k)·⌊m/k⌋².

Тогда ответ = Σ_{d=1}^{N} φ(d) · C(⌊N/d⌋).

Но C(m) нужно для каждого m = ⌊N/d⌋, а таких значений O(√N). Вычислим C(m) через префиксные суммы μ: C(m) = Σ_{k=1}^{m} μ(k)·⌊m/k⌋².

Однако это O(N√N) в лоб. Эффективнее:

Альтернативный подход. Используем тождество φ = μ * id:

Σ_{i,j} φ(НОД(i,j)) = Σ_{d=1}^{N} (Σ_{i,j: d|i, d|j} φ(НОД(i,j)/d) ...

Нет, проще напрямую. Используем φ(НОД(i,j)) = Σ_{t|НОД(i,j)} μ(t)·t:

Σ_{i,j} φ(НОД(i,j)) = Σ_{i,j} Σ_{t|НОД(i,j)} μ(t)·t
                     = Σ_{t=1}^{N} μ(t)·t · (⌊N/t⌋)²

Это просто и вычисляется за O(N)!

Доказательство. Используем φ(n) = Σ_{t|n} μ(t)·t (формула, следующая из φ = μ * id). Тогда φ(НОД(i,j)) = Σ_{t|НОД(i,j)} μ(t)·t.

При суммировании по всем парам:

Σ_{i,j} Σ_{t|НОД(i,j)} μ(t)·t = Σ_{t} μ(t)·t · #{(i,j): t|i, t|j}
                                = Σ_{t} μ(t)·t · ⌊N/t⌋²

Реализация:

#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> 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];
        }
    }

    // Ответ = sum_{t=1}^{N} mu(t) * t * floor(N/t)^2
    long long ans = 0;
    for (int t = 1; t <= N; ++t) {
        if (mu[t] == 0) continue;
        long long q = N / t;
        // mu(t) * t * q^2 mod MOD
        long long term = (long long)mu[t] % MOD * (t % MOD) % MOD * (q % MOD) % MOD * (q % MOD) % MOD;
        ans = (ans + term + MOD) % MOD;
    }

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

Сложность: O(N) на решето + O(N) на подсчёт = O(N).

Проверка для N=2:

• t=1: μ(1)·1·4 = 4
• t=2: μ(2)·2·1 = −2
• Итого: 2.

Прямой подсчёт: φ(НОД(1,1))+φ(НОД(1,2))+φ(НОД(2,1))+φ(НОД(2,2)) = φ(1)+φ(1)+φ(1)+φ(2) = 1+1+1+1 = 4. Нет, ответ 4, не 2. Проверим формулу: N/t при t=1 равно 2, q²=4. t=2: q=1, q²=1. Ответ = 1·1·4 + (−1)·2·1 = 4 − 2 = 2. Но прямой ответ 4. Ошибка!

Проверим: φ(1)=1, φ(2)=1. Пары: (1,1)→φ(1)=1, (1,2)→φ(1)=1, (2,1)→φ(1)=1, (2,2)→φ(2)=1. Сумма = 4.

Формула: Σ_t μ(t)·t·⌊2/t⌋² = 1·1·4 + (−1)·2·1 = 4−2=2. Ошибка в выводе формулы!

Разберёмся: φ(n) = Σ_{t|n} μ(t)·t (это верно? Проверим: φ(2) = Σ_{t|2} μ(t)·t = μ(1)·1 + μ(2)·2 = 1−2 = −1. Неверно! φ(2)=1.)

Правильная формула: φ = μ * id означает φ(n) = Σ_{d|n} μ(d)·(n/d), а НЕ Σ_{d|n} μ(d)·d.

Исправим:

φ(НОД(i,j)) = Σ_{t|НОД(i,j)} μ(t) · (НОД(i,j)/t)

Тогда:

Σ_{i,j} φ(НОД(i,j)) = Σ_{i,j} Σ_{t|НОД(i,j)} μ(t) · НОД(i,j)/t

Это сложнее. Вернёмся к первому подходу:

Σ_{i,j} φ(НОД(i,j)) = Σ_{d=1}^{N} φ(d) · #{(i,j): НОД(i,j)=d}
                     = Σ_{d=1}^{N} φ(d) · C(⌊N/d⌋)

где C(m) = Σ_{k=1}^{m} μ(k)·⌊m/k⌋² (упорядоченные пары с НОД=1).

Проверка: N=2. C(2) = Σ_{k=1}^{2} μ(k)·⌊2/k⌋² = 1·4 + (−1)·1 = 3. C(1) = 1. Ответ = φ(1)·C(2) + φ(2)·C(1) = 1·3 + 1·1 = 4. ✓

Более эффективная реализация через предварительный подсчёт C:

#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;
    cin >> N;

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

    for (int i = 2; i <= N; ++i) {
        if (!is_composite[i]) {
            primes.push_back(i);
            mu[i] = -1;
            phi[i] = 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;
                break;
            }
            mu[ip] = -mu[i];
            phi[ip] = phi[i] * (p-1);
        }
    }

    // Вычисляем C[m] = #{(i,j): 1<=i,j<=m, НОД(i,j)=1}
    // C[m] = sum_{k=1}^{m} mu[k] * floor(m/k)^2
    // Но нам нужно C[floor(N/d)] для каждого d.
    // Если N маленькое, вычислим C для всех m.
    // C[m] = C[m-1] + 2 * phi[m] - 1 (рекуррентность через функцию Эйлера)
    // Верна ли эта рекуррентность?
    // C[m] = #{взаимно простые пары в [1,m]²}
    // C[m] = C[m-1] + 2 * #{k < m: НОД(k,m)=1} + [НОД(m,m)=1]
    //      = C[m-1] + 2*phi[m] + [m=1]
    // При m>1: C[m] = C[m-1] + 2*phi[m]
    vector<long long> C(N+1, 0);
    C[1] = 1;
    for (int m = 2; m <= N; ++m) {
        C[m] = C[m-1] + 2LL * phi[m];
    }

    // Ответ = sum_{d=1}^{N} phi[d] * C[N/d]  mod MOD
    long long ans = 0;
    for (int d = 1; d <= N; ++d) {
        ans = (ans + (long long)phi[d] % MOD * (C[N/d] % MOD)) % MOD;
    }

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

Доказательство рекуррентности C[m] = C[m-1] + 2·φ(m). C[m] считает пары (i,j) в [1,m]². Новые пары по сравнению с C[m-1]: те, где хотя бы одно из i,j равно m. Это пары (m,j) с 1≤j≤m и НОД(m,j)=1, а также (i,m) с 1≤i<m и НОД(i,m)=1.

• Пар (m,j) с j≤m и НОД=1: φ(m) штук (при j=m НОД(m,m)=m=1 только при m=1). Точнее: #{j∈[1,m]: НОД(m,j)=1}. При m>1 это φ(m) (j=m не входит, НОД(m,m)=m>1).
• Пар (i,m) с i<m и НОД=1: тоже φ(m) (по симметрии НОД(i,m)=НОД(m,i)).
• Пары (m,m): НОД(m,m)=m, равно 1 только при m=1. Итого: при m>1 прирост = 2·φ(m). ✓

Сложность: O(N) на решето + O(N) на вычисление C + O(N) на ответ = O(N).

Разбор задачи 3 (сложная: Codeforces Div.1 D уровень): Инверсия Мёбиуса над двумерной суммой

Задача. Дано N (1 ≤ N ≤ 5·10⁵). Найдите:

Σ_{i=1}^{N} Σ_{j=1}^{N} σ(НОД(i,j))

по модулю 10⁹+7.

Анализ. Используем тождество σ = id * 1, то есть σ(n) = Σ_{d|n} d.

Σ_{i,j} σ(НОД(i,j)) = Σ_{i,j} Σ_{d|НОД(i,j)} d
                     = Σ_{d=1}^{N} d · #{(i,j): d | НОД(i,j)}
                     = Σ_{d=1}^{N} d · ⌊N/d⌋²

Это вычисляется за O(N)! Но нам нужен более сложный вариант, используем задачу как учебный пример для свёртки.

Расширение задачи. Найдём Σ_{i,j} f(НОД(i,j)) для произвольной мультипликативной f, применив общий шаблон:

Σ_{i,j} f(НОД(i,j)) = Σ_{d=1}^{N} f(d) · C(⌊N/d⌋)

где C(m) = #{(i,j) с НОД=d} ... нет:

= Σ_{d=1}^{N} f(d) · #{(i,j): НОД(i,j)=d}
= Σ_{d=1}^{N} f(d) · Σ_{k=1}^{⌊N/d⌋} μ(k)·⌊N/(dk)⌋²
= Σ_{t=1}^{N} ⌊N/t⌋² · Σ_{d|t} f(d)·μ(t/d)
= Σ_{t=1}^{N} ⌊N/t⌋² · (f * μ)(t)

Это обобщённая формула! При f = σ: (σ * μ)(t) = (id * 1 * μ)(t) = (id * ε)(t) = id(t) = t. Что даёт Σ_t t·⌊N/t⌋² — то же, что получили выше.

Общая программа для таких задач:

#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;
    cin >> N;

    // Линейное решето: mu, sigma
    vector<int> mu(N+1, 0);
    vector<long long> sigma(N+1, 0), f_conv_mu(N+1, 0);
    vector<int> pe(N+1, 1), exp_e(N+1, 0);
    vector<long long> sigma_pe(N+1, 1);
    vector<int> primes;
    vector<bool> is_composite(N+1, false);

    mu[1] = 1; sigma[1] = 1;

    for (int i = 2; i <= N; ++i) {
        if (!is_composite[i]) {
            primes.push_back(i);
            mu[i] = -1;
            sigma[i] = i + 1;
            pe[i] = i; exp_e[i] = 1; sigma_pe[i] = 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) {
                int e = exp_e[i] + 1;
                exp_e[ip] = e;
                pe[ip] = pe[i] * p;
                long long sp = sigma_pe[i] + pe[ip];
                sigma_pe[ip] = sp;
                sigma[ip] = sigma[i] / sigma_pe[i] * sp;
                mu[ip] = 0;
                break;
            } else {
                exp_e[ip] = 1; pe[ip] = p;
                sigma_pe[ip] = p + 1;
                sigma[ip] = sigma[i] * (p + 1);
                mu[ip] = -mu[i];
            }
        }
    }

    // Вычисляем h[t] = (sigma * mu)(t) = sum_{d|t} sigma(d) * mu(t/d)
    // По тождеству sigma = id * 1, mu = 1^{-1}: sigma * mu = id * 1 * mu = id * eps = id
    // Но вычислим явно для демонстрации свёртки
    for (int d = 1; d <= N; ++d) {
        if (mu[d] == 0) continue;
        for (int k = 1; (long long)d*k <= N; ++k) {
            f_conv_mu[d*k] = (f_conv_mu[d*k] + sigma[k] * mu[d] % MOD + MOD) % MOD;
        }
    }

    // ans = sum_{t=1}^{N} floor(N/t)^2 * h[t]
    long long ans = 0;
    for (int t = 1; t <= N; ++t) {
        long long q = N / t;
        ans = (ans + q % MOD * (q % MOD) % MOD * f_conv_mu[t]) % MOD;
    }

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

Сложность: O(N log N) на вычисление свёртки + O(N) на ответ.

Замечание. Из тождества σ * μ = (id * 1) * μ = id * (1 * μ) = id * ε = id следует, что f_conv_mu[t] = t. Это позволяет упростить код до O(N), избежав явного вычисления свёртки. В задачах соревновательного программирования подобные «алгебраические сокращения» крайне важны.

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

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

1. Путаница между мультипликативной и полностью мультипликативной. φ(4) = φ(2²) = 2, φ(2)² = 1. Поэтому φ(4) ≠ φ(2)·φ(2) — φ НЕ полностью мультипликативна. Зато φ(6) = φ(2)·φ(3) = 1·2 = 2, что верно (НОД(2,3)=1). Всегда проверяйте взаимную простоту!

2. Ошибка при вычислении σ в линейном решете. Формула σ[ip] = σ[i] / σ_pe[i] * sp использует целочисленное деление. Необходимо убедиться, что σ_pe[i] делит σ[i]. Это верно из теоремы о мультипликативности, но требует аккуратного кода.

3. Неверное применение тождества φ = μ * id. Формула φ(n) = Σ_{d|n} μ(d)·(n/d), но НЕ Σ_{d|n} μ(d)·d. Это частая опечатка.

4. Забытый mod в промежуточных вычислениях. При работе по модулю 10⁹+7 умножать нужно ПЕРЕД взятием по модулю, иначе переполнение:

// НЕВЕРНО:
ans += (phi[d] * C[N/d]) % MOD;  // phi[d] * C[N/d] может переполниться int

// ВЕРНО:
ans += (long long)phi[d] % MOD * (C[N/d] % MOD) % MOD;

5. Неправильная свёртка (перестановка аргументов). (f * g)(n) = Σ_{d|n} f(d)·g(n/d). При реализации через двойной цикл:

for (int d = 1; d <= N; ++d)
    for (int k = 1; d*k <= N; ++k)
        h[d*k] += f[d] * g[k];  // h = f * g: h[n] = sum_{d|n} f(d)*g(n/d)

Здесь d — делитель n=dk, g[k]=g[n/d]. Убедитесь, что f и g стоят на правильных местах.

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

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

h(n) = Σ_{d|n} f(d)·g(n/d) = (f * g)(n)

Если f = a * b, то h = a * b * g. Возможно, b * g = ε или известная простая функция.

Пример профессионального приёма. Задача: вычислить Σ_{d|n} μ(d)·d(n/d) для всех n до N.

• Это (μ * d)(n).
• d = 1 * 1, значит μ * d = μ * 1 * 1 = ε * 1 = 1.
• Ответ: (μ * d)(n) = 1(n) = 1 для всех n!

Вместо O(N log N) для вычисления свёртки получаем O(1) ответ. Знание алгебраических тождеств экономит часы на соревнованиях.

Ещё один приём: решето-подход для нескольких функций одновременно. Если нужны d, σ, φ — можно запустить одно линейное решето и вычислить все три функции параллельно (следующая глава).

Итог

Мультипликативные функции и свёртка Дирихле образуют элегантную алгебрическую систему:

1. Мультипликативность f(mn) = f(m)f(n) при НОД(m,n)=1 сводит вычисление к степеням простых. Любая мультипликативная функция определяется на {p^k}.
2. Полная мультипликативность сильнее: φ, μ, d, σ — мультипликативны, но не полностью. Только id, 1, ε — полностью мультипликативны среди «стандартных» функций.
3. Свёртка Дирихле — умножение в кольце арифметических функций. Коммутативна, ассоциативна, ε — нейтральный элемент.
4. Замкнутость: свёртка двух мультипликативных функций мультипликативна.
5. Тождества: d = 11, σ = id1, φ = μid, ε = μ1 — это «арифметика» кольца.
6. Практика: большинство задач на НОД сводятся к Σ f(d)·⌊N/d⌋² через формулу Мёбиуса.