Сайт переезжает. Большинство статей уже перенесено на новую версию.
Скоро добавим автоматические переходы, но пока обновленную версию этой статьи можно найти там.

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

Часто в задачах требуется посчитать что-то по простому модулю (чаще всего \(10^9 + 7\)). Это делают для того, чтобы участникам не приходилось использовать длинную арифметику, и они могли сосредоточиться на самой задаче.

Обычные арифметические операции выполняются не сильно сложнее — просто нужно брать модули и заботиться о переполнении. Например:

c = (a + b) % mod;
c = (mod + a - b) % mod;
c = a * b % mod;

Но вот с делением возникают проблемы — мы не можем просто взять и поделить. Пример: \(\frac{8}{2} = 4\), но \(\frac{8 \% 5 = 3}{2 \% 5 = 2} \neq 4\).

Нужно найти некоторый элемент, который будет себя вести как \(\frac{1}{a} = a^{-1}\), и вместо «деления» домножать на него. Назовем такой элемент обратным.

Способ 1: бинарное возведение в степень

Если модуль \(p\) простой, то решением будет \(a^{-1} \equiv a^{p-2}\). Это следует из малой теоремы Ферма:

Теорема. \(a^p \equiv a \pmod p\) для всех \(a\), не делящихся на \(p\).

Доказательство. (для понимания несущественно, можно пропустить)

\[ \begin{aligned} a^p &= (\underbrace{1+1+\ldots+1+1}_\text{$a$ раз})^p \\ &= \sum_{x_1+x_2+\ldots+x_a = p} P(x_1, x_2, \ldots, x_a) & \text{(раскладываем по определению)} \\ &= \sum_{x_1+x_2+\ldots+x_a = p} \frac{p!}{x_1! x_2! \ldots x_a!} & \text{(какие слагаемые не делятся на $p$?)} \\ &\equiv P(p, 0, \ldots, 0) + \ldots + P(0, 0, \ldots, p) & \text{(все остальные не убьют $p$ в знаменателе)} \\ &= a \end{aligned} \]

Здесь \(P(x_1, x_2, \ldots, x_n) = \frac{k}{\prod (x_i!)}\) это мультиномиальный коеффициент — количество раз, которое элемент \(a_1^{x_1} a_2^{x_2} \ldots a_n^{x_n}\) появится при раскрытии скобки \((a_1 + a_2 + \ldots + a_n)^k\).

Теперь два раза «поделим» наш результат на \(a\).

\[ a^p \equiv a \implies a^{p-1} \equiv 1 \implies a^{p-2} \equiv a^{-1} \]

Получается, что \(a^{p-2}\) ведет себя как \(a^{-1}\), что нам по сути и нужно. Посчитать \(a^{p-2}\) можно за \(O(\log p)\) бинарным возведением в степень.

Приведем код, который позволяет считает \(C_n^k\).

int t[maxn]; // факториалы, можно предподситать простым циклом

// бинарное возведение в степень
int bp (int a, int n) {
    int res = 1;
    while (n) {
        if (n & 1) res = res * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return res;
}

// находит обратный элемент как a^(p-2)
int inv (int x) {
    return bp(x, mod-2);
}

int c (int n, int k) {
    return t[n] * inv(t[k]) % mod * inv(t[n-k]) % mod;
}

Способ 2: диофантово уравнение

Диофантовыми уравнениями называют такие штуки:

\[ ax + by = 1 \]

Требуется решить их в целых числах, то есть \(a\) и \(b\) известны, и нужно найти такие целые (возможно, отрицательные) \(x\) и \(y\), чтобы равенство выполнялось. Решают такие вещи расширенным алгоритмом Евклида. TODO: описать, как он работает.

Подставим в качестве \(a\) и \(b\) соответственно \(a\) и \(m\)

\[ ax + my = 1 \]

Одним из решений уравнения и будет \(a^{-1}\), потому что если взять уравнение по модулю \(m\), то получим

\[ ax + by = 1 \iff ax \equiv 1 \iff x \equiv a^{-1} \pmod m \]

Преимущества этого метода над возведением в степень:

Сам автор почти всегда использует возведение в степень.

Почему \(10^9+7\)?

  1. Это выражение довольно легко вбивать (1e9+7).
  2. Простое число.
  3. Достаточно большое.
  4. int не переполняется при сложении.
  5. long long не переполняется при умножении.

Кстати, \(10^9 + 9\) обладает теми же свойствами. Иногда используют и его.

Предподсчёт обратных факториалов за линейное время

Пусть нам нужно зачем-то посчитать все те же \(C_n^k\), но для больших \(n\) и \(k\), поэтому асимптотика \(O(n \log m)\) нас не устроит. Оказывается, мы можем сразу предподсчитать все обратные ко всем факториалам.

Если у нас уже написан inv, то нам не жалко потратить \(O(\log m)\) операций, посчитав \(m!^{-1}\).

После этого мы будем считать \((m-1)!^{-1}\) как \(m!^{-1} m = \frac{1}{1 \cdot 2 \cdot \ldots \cdot (m-1)}\).

int f[maxn];
f[0] = 1;
for (int i = 1; i < maxn; i++)
    f[i] = i*f[i-1] % mod;

int r[maxn];
r[maxn-1] = inv(f[maxn-1])
for (int i = maxn-1; i >= 1; i--)
    r[i-1] = r[i]*i % mod;

TODO: техника с сайта емакса.