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

Почти каждая задача на комбинаторику или теорию чисел в соревновательном программировании заканчивается словами «выведите ответ по модулю 10^9 + 7». Это не случайность: ответы растут экспоненциально, и хранить их целиком невозможно. Но чтобы работать с остатками правильно, нужно понимать, что сложение и умножение «переносятся» через операцию взятия остатка, а вот деление — нет, и именно здесь скрыта большая часть ошибок новичков.

Модульная арифметика — это фундамент, на котором стоит вся вычислительная теория чисел: проверка простоты, дискретное логарифмирование, быстрое преобразование Фурье (NTT), криптография RSA. Без твёрдого понимания остатков невозможно разобраться ни в одной из этих тем.

Теория

Сравнения и кольцо вычетов Z_n

Определение. Говорят, что целые числа a и b сравнимы по модулю n (записывается a ≡ b (mod n)), если n делит разность a − b, то есть существует целое k такое, что a − b = k·n.

Отношение сравнимости по модулю n является отношением эквивалентности: оно рефлексивно, симметрично и транзитивно. Множество классов эквивалентности называется кольцом вычетов и обозначается Z_n или Z/nZ.

Представители классов — числа 0, 1, 2, …, n−1. Каждое целое число принадлежно ровно одному классу: его остаток от деления на n.

Операции в Z_n:

• (a + b) mod n = ((a mod n) + (b mod n)) mod n
• (a − b) mod n = ((a mod n) − (b mod n) + n) mod n ← прибавляем n, чтобы не уйти в отрицательные числа
• (a · b) mod n = ((a mod n) · (b mod n)) mod n

Эти три формулы вытекают из простого наблюдения: если a = q₁n + r₁ и b = q₂n + r₂, то a + b = (q₁ + q₂)n + (r₁ + r₂), поэтому остаток суммы зависит только от r₁ и r₂.

Почему Z_p при простом p особенный? Потому что в нём каждый ненулевой элемент имеет мультипликативный обратный, то есть Z_p является полем. В Z_6, например, элемент 2 не имеет обратного: нет такого x, что 2x ≡ 1 (mod 6).

Обратный элемент по модулю

Определение. Число x называется обратным к a по модулю n, если a · x ≡ 1 (mod n). Обозначается a⁻¹ (mod n).

Обратный элемент существует тогда и только тогда, когда gcd(a, n) = 1.

Доказательство. Уравнение ax ≡ 1 (mod n) эквивалентно ax − ny = 1 для некоторого целого y. По теореме Безу это линейное диофантово уравнение имеет решение тогда и только тогда, когда gcd(a, n) делит 1, то есть gcd(a, n) = 1.

Существует три основных метода вычисления обратного элемента.

Метод 1: расширенный алгоритм Евклида

Расширенный алгоритм Евклида находит x, y такие, что ax + ny = gcd(a, n). Если gcd(a, n) = 1, то ax ≡ 1 (mod n), значит x и есть a⁻¹.

Рекурсивные формулы:

gcd(a, b): если b = 0, то (gcd=a, x=1, y=0)
иначе: (g, x₁, y₁) = gcd(b, a mod b)
        x = y₁,  y = x₁ − (a div b)·y₁

Доказательство корректности: на каждом шаге поддерживается инвариант gcd(a,b) = ax + by. В конце b = 0, gcd = a, x = 1, y = 0, инвариант выполнен.

Сложность: O(log min(a, n)).

Метод 2: малая теорема Ферма

Если p — простое и p ∤ a, то a^(p−1) ≡ 1 (mod p). Значит a · a^(p−2) ≡ 1 (mod p), то есть a⁻¹ ≡ a^(p−2) (mod p).

Доказательство малой теоремы Ферма. Рассмотрим множество S = {a, 2a, 3a, …, (p−1)a} по модулю p. Докажем, что это перестановка множества {1, 2, …, p−1}. Предположим, ia ≡ ja (mod p) при 1 ≤ i < j ≤ p−1. Тогда p | (j−i)a. Так как gcd(a, p) = 1 и 0 < j−i < p, это невозможно. Значит, все элементы S различны и ненулевые — перестановка. Перемножим все элементы S и всего {1, …, p−1}:

a^(p−1) · (p−1)! ≡ (p−1)! (mod p)

Так как (p−1)! взаимно просто с p, сокращаем: a^(p−1) ≡ 1 (mod p). ∎

Сложность: O(log p) — через быстрое возведение в степень.

Метод работает только для простых модулей!

Метод 3: линейное решето обратных

Если нам нужны обратные для всех чисел от 1 до N, можно воспользоваться рекуррентным соотношением. Запишем p = (p div i) · i + (p mod i), то есть:

(p div i) · i + (p mod i) ≡ 0 (mod p)

Обозначим q = p div i, r = p mod i. Тогда qi + r ≡ 0 (mod p), умножим на i⁻¹ · r⁻¹:

q · r⁻¹ + i⁻¹ ≡ 0 (mod p)

i⁻¹ ≡ −q · r⁻¹ ≡ −(p div i) · inv[p mod i] (mod p)

Поскольку p mod i < i, inv[p mod i] уже вычислен! Инициализируем inv[1] = 1 и итеративно считаем остальные за O(1) каждый.

Сложность: O(N) для всех чисел 1…N.

Теорема Эйлера

Теорема Эйлера. Если gcd(a, n) = 1, то a^φ(n) ≡ 1 (mod n), где φ(n) — функция Эйлера (количество натуральных чисел, не превосходящих n и взаимно простых с n).

Доказательство. Рассмотрим группу единиц Z_n* — множество классов вычетов {r₁, r₂, …, r_φ(n)}, взаимно простых с n. Умножим каждый элемент на a (gcd(a,n)=1): множество {ar₁, ar₂, …, ar_φ(n)} mod n является перестановкой Z_n* (по аналогии с доказательством теоремы Ферма). Перемножая:

a^φ(n) · r₁r₂…r_φ(n) ≡ r₁r₂…r_φ(n) (mod n)

Произведение r₁…r_φ(n) взаимно просто с n (как произведение элементов, взаимно простых с n), поэтому можно сократить: a^φ(n) ≡ 1 (mod n). ∎

Малая теорема Ферма — частный случай при n = p простом, ибо φ(p) = p − 1.

Быстрое возведение в степень

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

a^n = a^(b_k·2^k + … + b_1·2 + b_0) = ∏ a^(b_i · 2^i)

Вычисляем последовательно a, a², a⁴, a⁸, … (квадратируя), и перемножаем только те степени, где бит установлен. Всего O(log n) умножений.

Доказательство корректности по индукции: при n=0 результат 1. При n>0: если n чётное, a^n = (a^(n/2))²; если нечётное, a^n = a · a^(n−1). Рекурсия работает корректно на всех шагах.

Класс ModInt

Инкапсуляция арифметики по модулю в класс позволяет писать читаемый код без постоянного % MOD и защищает от забытых взятий остатка. Для соревнований критично, чтобы класс был header-only, без виртуальных методов, и компилировался с -O2 без предупреждений.

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

---

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

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

// ─────────────────────────────────────────────────────────────────
// Расширенный алгоритм Евклида
// Возвращает gcd(a, b), заполняет x, y такие, что a*x + b*y = gcd
// ─────────────────────────────────────────────────────────────────
long long extended_gcd(long long a, long long b, long long &x, long long &y) {
    if (b == 0) { x = 1; y = 0; return a; }
    long long x1, y1;
    long long g = extended_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

// Обратный элемент через расш. Евклид; требует gcd(a, mod) == 1
// Сложность: O(log mod)
long long inv_ext_gcd(long long a, long long mod) {
    long long x, y;
    long long g = extended_gcd(a % mod, mod, x, y);
    assert(g == 1);  // gcd должен быть 1
    return (x % mod + mod) % mod;
}

// ─────────────────────────────────────────────────────────────────
// Быстрое возведение в степень по модулю
// Сложность: O(log e)
// ─────────────────────────────────────────────────────────────────
long long power(long long base, long long exp, long long mod) {
    long long result = 1;
    base %= mod;
    if (base < 0) base += mod;  // корректируем отрицательное основание
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

// Обратный элемент через малую теорему Ферма (только для простого mod!)
// Сложность: O(log mod)
long long inv_ferma(long long a, long long p) {
    return power(a, p - 2, p);
}

// ─────────────────────────────────────────────────────────────────
// Линейное решето обратных 1..N по простому модулю p
// Сложность: O(N) времени, O(N) памяти
// ─────────────────────────────────────────────────────────────────
vector<long long> linear_inverse_sieve(int N, long long p) {
    vector<long long> inv(N + 1, 0);
    inv[1] = 1;
    for (int i = 2; i <= N; ++i) {
        // inv[i] = -(p / i) * inv[p % i] % p
        inv[i] = (p - (p / i) * inv[p % i] % p) % p;
    }
    return inv;
}

// ─────────────────────────────────────────────────────────────────
// Класс ModInt — шаблон с модулем как параметром
// Использование: using Mint = ModInt<1'000'000'007LL>;
// ─────────────────────────────────────────────────────────────────
template<long long MOD>
struct ModInt {
    static_assert(MOD >= 2, "MOD must be at least 2");
    long long v;

    ModInt() : v(0) {}
    ModInt(long long x) : v(x % MOD) { if (v < 0) v += MOD; }

    // Арифметические операторы
    ModInt operator+(const ModInt& o) const { return ModInt(v + o.v); }
    ModInt operator-(const ModInt& o) const { return ModInt(v - o.v + MOD); }
    ModInt operator*(const ModInt& o) const { return ModInt(v * o.v); }
    ModInt operator/(const ModInt& o) const { return *this * o.inv(); }

    ModInt& operator+=(const ModInt& o) { return *this = *this + o; }
    ModInt& operator-=(const ModInt& o) { return *this = *this - o; }
    ModInt& operator*=(const ModInt& o) { return *this = *this * o; }
    ModInt& operator/=(const ModInt& o) { return *this = *this / o; }

    bool operator==(const ModInt& o) const { return v == o.v; }
    bool operator!=(const ModInt& o) const { return v != o.v; }

    // Возведение в степень за O(log e)
    ModInt pow(long long e) const {
        ModInt result(1), base(v);
        if (e < 0) { base = base.inv(); e = -e; }
        while (e > 0) {
            if (e & 1) result *= base;
            base *= base;
            e >>= 1;
        }
        return result;
    }

    // Обратный элемент (MOD должен быть простым для метода Ферма)
    // Если нужен для составного модуля — заменить на extended_gcd
    ModInt inv() const {
        return pow(MOD - 2);
    }

    friend ostream& operator<<(ostream& os, const ModInt& m) {
        return os << m.v;
    }
    friend istream& operator>>(istream& is, ModInt& m) {
        long long x; is >> x; m = ModInt(x); return is;
    }
};

// Удобный псевдоним для соревнований
using Mint = ModInt<1'000'000'007LL>;

// ─────────────────────────────────────────────────────────────────
// Факториалы и обратные к ним — для комбинаций
// ─────────────────────────────────────────────────────────────────
constexpr int MAXN = 2'000'001;
Mint fact[MAXN], inv_fact[MAXN];

void precompute_factorials() {
    fact[0] = 1;
    for (int i = 1; i < MAXN; ++i) fact[i] = fact[i-1] * i;
    inv_fact[MAXN-1] = fact[MAXN-1].inv();
    for (int i = MAXN-2; i >= 0; --i) inv_fact[i] = inv_fact[i+1] * (i+1);
}

Mint C(int n, int k) {
    if (k < 0 || k > n) return 0;
    return fact[n] * inv_fact[k] * inv_fact[n-k];
}

// ─────────────────────────────────────────────────────────────────
// Демонстрация (main для проверки корректности)
// ─────────────────────────────────────────────────────────────────
int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    const long long MOD = 1'000'000'007LL;

    // Проверка inv_ext_gcd
    long long a = 123456789LL;
    long long inv_a = inv_ext_gcd(a, MOD);
    assert(a * inv_a % MOD == 1);

    // Проверка power / inv_ferma
    assert(inv_ferma(a, MOD) == inv_a);

    // Проверка линейного решета
    auto inv = linear_inverse_sieve(10, MOD);
    for (int i = 1; i <= 10; ++i) {
        assert((long long)i * inv[i] % MOD == 1);
    }

    // Проверка ModInt
    Mint x = 3, y = 7;
    cout << x + y << "\n";  // 10
    cout << x * y << "\n";  // 21
    cout << x.inv() << "\n";  // 333333336 (3^{-1} mod 10^9+7)

    // Проверка биномиальных коэффициентов
    precompute_factorials();
    cout << C(10, 3) << "\n";  // 120
    cout << C(50, 25) << "\n"; // 126410606777 mod 10^9+7

    return 0;
}

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

Разбор задачи 1 — Биномиальный коэффициент по модулю простого числа

Условие. Вычислить C(n, k) mod p, где p = 10^9 + 7, 0 ≤ k ≤ n ≤ 2·10^6.

Разбор. Прямое вычисление C(n, k) = n! / (k! · (n−k)!) даёт астрономические числа. Ключевое наблюдение: так как p простое и p > n, в знаменателе нет кратных p — можно пользоваться обратными элементами.

Алгоритм:

1. Предвычислить fact[0..n] за O(n).
2. Предвычислить inv_fact[0..n]: сначала inv_fact[n] = fact[n]^(p−2) за O(log p), затем inv_fact[i] = inv_fact[i+1] * (i+1) за O(1) каждый (используя тождество (i!)⁻¹ = ((i+1)!)⁻¹ · (i+1)).
3. C(n, k) = fact[n] * inv_fact[k] * inv_fact[n−k] mod p.

Почему трюк с inv_fact работает? Из (i+1)! = i! · (i+1) следует inv((i+1)!) = inv(i!) · inv(i+1), откуда inv(i!) = inv((i+1)!) · (i+1). Мы вычисляем массив справа налево, начиная от inv(n!), которое получаем одним возведением в степень.

precompute_factorials();  // из блока кода выше
cout << C(2000000, 1000000) << "\n";  // моментально

Итог: O(n) предобработка, O(1) каждый запрос. Это стандартный шаблон, который используется в 80% задач на комбинаторику.

Разбор задачи 2 — Сумма геометрической прогрессии по модулю

Условие. Дано a, r, n, m (1 ≤ n ≤ 10^18, m — простое). Найти S = a + ar + ar² + … + ar^(n−1) mod m.

Разбор. Формула S = a · (r^n − 1) / (r − 1) при r ≠ 1. Но n может быть до 10^18, и r − 1 может не иметь обратного... нет, m простое и r ≢ 1 (mod m), поэтому r−1 взаимно просто с m (если r−1 ≢ 0 mod m). Случай r ≡ 1 (mod m): S = n · a mod m.

Более элегантный подход — рекуррентный метод «быстрого суммирования»:

Пусть f(n) = 1 + r + r² + … + r^(n−1).

• Если n чётное: f(n) = f(n/2) · (1 + r^(n/2))
• Если n нечётное: f(n) = f(n−1) + r^(n−1)

Это работает за O(log² n) (или O(log n) с предвычислением степеней).

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

long long MOD;

long long power(long long base, long long exp, long long mod) {
    long long result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

// Сумма 1 + r + r^2 + ... + r^(n-1) mod MOD
// Используем разделяй-и-властвуй: O(log^2 n)
long long geo_sum(long long r, long long n) {
    if (n == 0) return 0;
    if (n == 1) return 1 % MOD;
    if (n % 2 == 0) {
        // f(n) = f(n/2) * (1 + r^(n/2))
        long long half = geo_sum(r, n / 2);
        long long rh = power(r, n / 2, MOD);
        return half % MOD * ((1 + rh) % MOD) % MOD;
    } else {
        // f(n) = f(n-1) + r^(n-1)
        return (geo_sum(r, n - 1) + power(r, n - 1, MOD)) % MOD;
    }
}

int main() {
    long long a, r, n;
    cin >> a >> r >> n >> MOD;
    // S = a * geo_sum(r, n)
    cout << a % MOD * geo_sum(r % MOD, n) % MOD << "\n";
    return 0;
}

Ловушки:

• При n = 0 ответ 0.
• При r = 0: S = a (только первый член).
• При r ≡ 1 (mod m): S = n * a mod m — нужна отдельная ветка.
• Не забываем брать r % MOD перед передачей в функцию.

Разбор задачи 3 — Дискретный логарифм (Baby-step Giant-step)

Условие. Дано a, b, p (p простое). Найти наименьшее неотрицательное x такое, что a^x ≡ b (mod p), или сообщить, что решения нет. (Условие уровня Codeforces Div.1 D / ICPC Region)

Разбор. Метод «шаг младенца — шаг великана» (BSGS).

Пусть m = ⌈√p⌉. Запишем x = i·m − j, где 0 ≤ i, j ≤ m. Тогда:

a^(im − j) ≡ b (mod p) (a^m)^i ≡ b · a^j (mod p)

Алгоритм:

1. Baby steps: вычислить b · a^j mod p для j = 0, 1, …, m. Сохранить в хеш-таблицу: значение → j.
2. Giant steps: для i = 1, 2, …, m вычислить (a^m)^i mod p. Если встретилось значение из таблицы, нашли x = i·m − j.

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

long long power(long long base, long long exp, long long mod) {
    long long result = 1; base %= mod;
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod; exp >>= 1;
    }
    return result;
}

// Возвращает наименьший x >= 0: a^x == b (mod p), или -1
// Требует: p простое, gcd(a, p) = 1
// Сложность: O(sqrt(p) * log(p))
long long bsgs(long long a, long long b, long long p) {
    a %= p; b %= p;
    if (a == 0) return (b == 0) ? 1 : (b == 1 ? 0 : -1);

    long long m = (long long)ceil(sqrt((double)p)) + 1;

    // Baby steps: сохраняем b * a^j -> j
    unordered_map<long long, long long> table;
    long long cur = b;
    for (long long j = 0; j <= m; ++j) {
        table[cur] = j;
        cur = cur * a % p;
    }

    // Giant steps: ищем (a^m)^i в таблице
    long long am = power(a, m, p);
    cur = am;  // i = 1
    for (long long i = 1; i <= m + 1; ++i) {
        auto it = table.find(cur);
        if (it != table.end()) {
            long long x = i * m - it->second;
            if (x >= 0) return x;
        }
        cur = cur * am % p;
    }
    return -1;
}

int main() {
    long long a, b, p;
    cin >> a >> b >> p;
    long long ans = bsgs(a, b, p);
    if (ans == -1) cout << "No solution\n";
    else cout << ans << "\n";
    return 0;
}

Почему сложность O(√p)? Мы делаем m ≈ √p baby-шагов и m giant-шагов, каждый — O(1) при хеш-таблице (амортизированно) плюс O(log p) для power.

Подводные камни:

• При b = 1 и x = 0: нужно явно проверить, иначе цикл может пропустить.
• unordered_map имеет худший случай O(n) — для гарантии можно использовать хеш с солью или robin-hood hashing.
• Обобщённый BSGS (когда gcd(a, p) > 1) требует дополнительного этапа «обнуления» кратных gcd — это материал следующего уровня.

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

---

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

1. Отрицательные остатки после вычитания.

// Неверно:
long long res = (a - b) % MOD;  // может быть отрицательным!
// Верно:
long long res = ((a - b) % MOD + MOD) % MOD;

2. Переполнение при умножении.

// При MOD ~ 10^9 и a, b < MOD: a * b ~ 10^18, влезает в long long (max ~9.2 * 10^18)
// Но при MOD ~ 10^18 (например, в задачах с 128-битной арифметикой):
__int128 res = (__int128)a * b % MOD;  // используем __int128

3. Применение Ферма к составному модулю. Если модуль не простое число, a^(p−2) не является обратным! Нужен расширенный Евклид. Всегда проверяйте условие.

4. Деление нуля при вычислении обратного. inv(0) не существует. Если вы вычисляете C(n, k) и n или k равны нулю, нужна отдельная проверка.

5. Забытый % MOD в промежуточных вычислениях. Даже если все входные данные < MOD, произведение двух таких чисел уже > MOD² ≈ 10^18. После каждого умножения берите остаток.

6. Неправильная инициализация решета. В линейном решете inv[0] и inv[1] инициализируются отдельно; запуск цикла с i=0 даёт деление на ноль.

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

На соревнованиях держите в шаблоне готовый ModInt и предвычисленные факториалы. Большинство команд ICPC теряют время именно на написание/отладку модульной арифметики под давлением.

Второй практический совет: всегда тестируйте решение на a = 0, a = MOD−1, n = 0, n = 1 — это граничные случаи, где чаще всего появляются WA (Wrong Answer).

Третье наблюдение: если задача говорит «выведите ответ по модулю M», и M не дано простым, проверьте, можно ли применить КТО (Китайскую теорему об остатках — см. следующую главу), разбив M на простые множители.

Итог

Мы построили полную арифметику в кольце вычетов Z_n: сложение, вычитание, умножение и деление (через обратный элемент). Три метода нахождения обратного — расширенный Евклид (O(log n), универсальный), теорема Ферма (O(log p), только для простых), линейное решето (O(N) для всех до N) — покрывают все практические случаи.

Теорема Эйлера обобщает малую теорему Ферма и является ключевым инструментом для работы со степенями при больших показателях. Быстрое возведение в степень за O(log e) — один из базовых алгоритмических примитивов, который встречается во всей теории чисел.

Класс ModInt и предвычисленные факториалы с обратными — стандартный соревновательный шаблон, ускоряющий написание чистого и корректного кода.