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

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

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

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

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

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

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

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

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

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

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

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

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

$$a^p\equiv a⟹a^{p-1}\equiv 1⟹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⟺ax\equiv 1⟺x\equiv a^{-1}(\mathrm{mod}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 \dots \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: техника с сайта емакса.