本文中模运算仅对 ​\mathbb{Z} 中的数定义。

前言

在组合数学中，常常会遇到特别大的答案，这些答案极有可能超出 long long甚至 int128的范围。为判断答案的正确性，题目会要求将答案取模。因模运算的性质，程序可以在运算过程中对阶段答案不断取模，保证运算结果在任何时候不超出数据类型的范围。

基本概念

模运算英文名为​\text{Modular Arithmetic}，“模”是音译的，模运算的实际是求余运算，是一种二元运算，符号为​\bmod，​a模​b记为"​a \bmod b"，运算结果为​a除以​b的余数。

形如​a \equiv b \pmod m的式子称为同余式，表示​a \bmod m = b \bmod m，也可理解为​a与​b的差为​m的倍数。同余式满足自反律（​a \equiv a）、对称律（​a \equiv b \Rightarrow b \equiv a）、传递律（​a \equiv b \equiv c \Rightarrow a \equiv c）。

基本运算性质

‍

若​a\equiv b{\pmod {m}}，​c\equiv d{\pmod {m}}，则有

\begin{align}a + c\equiv b + d\pmod m\tag{i.1.1}\\ a - c\equiv b - d\pmod m\tag{i.1.2}\\ a \times c\equiv b \times d\pmod m\tag{i.1.3}\end{align}

​\text{(i.1.3)}可推导出

a^n \equiv b^n \pmod{m} \tag{i.1.4}

由这些基本运算性质可知，若题目求答案对​m取模的值，在每一步加减乘运算之后对结果取模，最终的结果等于直接对答案取模。这样就能保证运算过程中的数和最终结果都不超出数据类型的范围。

注意：在模数下进行减法运算时，可能遇到结果为负数的情况，C++规定模运算的结果的正负性和被模数的正负性相同，所以可能导致最终答案为负，通常的做法是在减法后加上模数再模上模数。

除法和逆元

上述运算性质中并没有关于除法的性质。

事实上，​ad \equiv bd \Rightarrow a \equiv b \pmod{m}并不一定成立。小学数学告诉我们，除以一个数等于乘上这个数的倒数。但​d \in \mathbb{Z}，所以​\frac{1}{d} \in [-1, 1]，而​\text{(i.1.3)}是在​a,b,c,d \in \mathbb{Z}的情况下推导出来的，所以不能使用​\text{(i.1.3)}的性质。

如果我们能在​\mathbb{Z}中找到一个可以代替​\frac{1}{d}的数，​ad \equiv bd \Rightarrow a \equiv b \pmod{m}自然就成立了。
不妨设设该数为​d'，如果需要满足​d'd=1，显然找不到代替​\frac{1}{d}的数；但根据​\text{(i.1.3)}，我们只需要满足​d'd \equiv 1 \pmod{m}就可以了。
​d'd \equiv 1 \pmod{m}可化为​d'd+m'm=1（其中​m' \in \mathbb{Z}，​m'和​d'均为未知数）。

裴蜀定理（ ​\text{Bézout's Lemma} ）： 对于​a,b \in \mathbb{Z^+}且​a,b不全为​0，方程​ax+by=c有解的充要条件为​gcd(a,b)|c。

对于​d'd+m'm=1，因为此时​c=1，所以该方程有整数解当且仅当​gcd(d, m)=1，即​d \bot m。

所以，若​d \bot m，

ad \equiv bd \Rightarrow a \equiv b \pmod{m} \tag{i.2}

‍

现在我们考虑在程序运算过程中出现除法如何处理。

设​res为该次运算不取模的结果，​a为未取模的被除数，​a'为已取模的被除数，​d为除数，则有​\frac{a}{d}=res \Rightarrow a=res \cdot d，​a' \equiv res \cdot d \pmod{m}。
接下来考虑已知​a'和​d，求​res \bmod m。
如果​d \bot m，且​d|a'，根据我们推导出的​\text{(i.2)}，​res \equiv \frac{a'}{d} \pmod{m}。
但大多数情况下​d|a'并不成立，这时我们像前面推导​\text{(i.2)}时一样考虑将同余式两边同乘代替​\frac{1}{d}的数​d'，就得到了​a' \cdot d' \equiv res \pmod{m}。

我们称​d'为​d的逆元。因为它的作用很像​d的倒数（我们其实就是拿它在​\text{mod } m下代替倒数的），即​d^{-1}，所以我们一般也把它记作​d^{-1}。

如果一个数​d存在逆元​d^{-1}，则​d^{-1} 不唯一（裴蜀定理还有一条内容：若二元不定整数方程​ax+by=c有整数解，则该方程有无数多组整数解），一般我们求的​d^{-1}是一个特解。

前面已经证明​d^{-1}d+m'm=1有解当且仅当​d \bot m，所以​d^{-1} 仅在​d \bot m 时存在。通常情况下，题目给的模数​m都是比较大的质数，只要​d<m，那么​d \bot m。即便​m不是质数，数据也会保证​d \bot m，或者不会出现除法运算。

逆元求法

扩展欧几里得法 EXGCD

根据我们前面的推导过程，​d^{-1}是​d^{-1}d \equiv 1 \pmod{m}，即​d^{-1}d+m'm=1（​m \in \mathbb{Z}）的解。所以我们直接解这个方程就可以了。

裴蜀定理用于判断二元不定整数方程是否有整数解，而它实际也给了我们一个解二元不定整数方程的启示。因为二元不定整数方程​ax+by=c有解的充要条件为​gcd(a,b)|c，对于任意一个这样的方程，我们只要求得一组​ax+by=gcd(a,b)的解，乘上​\frac{c}{gcd(a,b)}就得到原方程的一组特解了（求出一组特解，通解也就出来了）。对于​ax+by=gcd(a,b)，其实我们在用欧几里得法求​gcd(a,b)的时候就可以一并求出了。这个用欧几里得法求二元不定整数方程的整数解的方法称为扩展欧几里得法。

回到原问题，求​d^{-1}d+m'm=1（​m \in \mathbb{Z}）的解。
该方程有解时​gcd(d,m)=1，直接求​d^{-1}d+m'm=gcd(d,m)就行了，不需要乘倍数了。

欧几里得算法通过下面的递归式求解：

\begin{cases}gcd(0,n)=n\\gcd(m,n)=gcd(n \bmod m, m) & m>0\end{cases}

接下来我们考虑在求​gcd(m,n)时如何一并求出满足​m'm+n'n=gcd(m,n)的一组整数​m',n'。

  1. 递归到​gcd(0,n)的情况，直接令​m'=0，​n'=1
  2. 递归到​gcd(m,n)的情况，
    此时要继续递归​gcd(n \bmod m, m)，我们设​r为​n \bmod m，则下一层递归为​gcd(r, m)，​gcd(r,m)也会返回一组​n',m'，为方便区分，我们把它们分别记为​\overline{r},\overline{m}。因为是递归的，返回的​\overline{r},\overline{m}是​\overline{r}r+\overline{m}m=gcd(r,m)的一组解。我们现在要用​\overline{r},\overline{m}表示当前这一层递归​m',n'。因为​gcd(r,m)=gcd(m,n)，而​\overline{r}r+\overline{m}m=gcd(r,m)，​m'm+n'n=gcd(m,n)，可建立等式​\overline{r}r+\overline{m}m=m'm+n'n。前面我们令​r=n \bmod m，现在把它化成底函数的形式，​r=n-\lfloor n/m \rfloor m，再代入等式，得到​\overline{r}(n-\lfloor n/m \rfloor m)+\overline{m}m=m'm+n'n，再整理得到​(\overline{m}-\lfloor n/m \rfloor \overline{r})m+\overline{r}n=m'm+n'n，于是​\begin{cases}m'=\overline{m}-\lfloor n/m \rfloor \overline{r} \\ n'=\overline{r}\end{cases}，这就是转移方程。

现在我们可以用这个方法解​d^{-1}d+m'm=gcd(d,m)从而求得​d^{-1}了。

代码如下：

void exgcd(int a, int b, int& x, int& y) {
  if (b == 0) {
    x = 1, y = 0;
    return;
  }
  exgcd(b, a % b, y, x);
  y -= a / b * x;
}

（代码来自OI Wiki）

​T(n)=O(\log n)（和欧几里得法求​gcd相同）

快速幂法

费马小定理（ ​\text{Fermat's Little Theorem} ） ：若​p为质数，​\gcd(a, p) = 1，则​a^{p - 1} \equiv 1 \pmod{p}。

若模数​m为质数，且​d \bot m，就有​d^{m-1} \equiv 1 \pmod{m}。
联立方程​\begin{cases}d^{m-1} \equiv 1 \\ d^{-1}d \equiv 1 \end{cases}\pmod{m}，得​d^{-1}d \equiv d^{m-1} \pmod{m}，从而​d^{-1} \equiv d^{m-2}，用快速幂求就可以了。

注意：只有在模数​m为质数的情况下才能用快速幂求逆元。

代码就是快速幂的代码：

int qpow(long long a, int b) {
  int ans = 1;
  a = (a % p + p) % p;
  for (; b; b >>= 1) {
    if (b & 1) ans = (a * ans) % p;
    a = (a * a) % p;
  }
  return ans;
}
int GetInv(long long d, int m){
	return qpow(d,m-2);
}

​T(n)=O(\log n)（快速幂的时间复杂度）

线性求逆元

该方法可在​O(n)的时间复杂度下求出​d \in \mathbb{Z} \cap [1,n]对应的​d^{-1}。

线性算法通常是递推。

警告：因为我不知道OI Wiki上那个方法是怎么想到的，所以以下求递推式的过程是我瞎编的，谨慎阅读！

设​f(x)=x^{-1}，有​f(1) \equiv 1 \pmod{m}，因为​\forall m \in \mathbb{Z}，​1 \times 1 \equiv 1 \pmod{m}，这是递推的基础。

这个转移方程不好求。我们要设法找到​f(i)与​f(j)（​j<i）的关系。显然​f(i)的值与​i,j,f(j),m有关。 因为​j<i，在计算​f(i)之前​f(j)就已经被计算出来，所以只要确定​j就行了。我们设​f(i) \equiv R(i,j,m) \pmod{m}（​R是一个我瞎编的NB的函数，假设我们已经确定​j，它将​i,j,f(j),m以某种我们猜不到的方式进行运算，得到一个关于模​m与​f(i)同余的结果），那么​f(i)-R(i,j,m) \equiv 0 \pmod{m}，​f(i)-R(i,j,m)一定为​m的倍数。我们尝试将​m以某种NB的方式拆成一堆与​f(i),i,j,f(j),m有关的东西的和。首先我们现在只知道​i，所以我们先尝试强行把​m拆成和​i有关的东西，我们大胆地瞎猜合理地猜想​m可分解为​ki+b这种形式，我们现在只能确定​i，所以​k,b也应该与​i有关，又因为我们要求的东西与同余有关，所以我们大胆地把​m拆成​\lfloor m/i \rfloor i + m \bmod i，则​\begin{cases}k=\lfloor m/i \rfloor \\ b=m \bmod i\end{cases}，这样保证了​b<i，所以​f(b)也是知道的，后面可以用。显然​ki+b \equiv 0 \pmod{m}，我们尝试让​i^{-1}出现，将同余式两边同乘​i^{-1}，则​k+bi^{-1} \equiv 0 \pmod{m}，然后发现​i^{-1}还有​b这个系数，就再乘上​b^{-1}把​b消掉，​b^{-1}就是​f(b)，是已经知道的，可以放心用，然后就有​kf(b)+i^{-1} \equiv 0 \pmod{m}，其中​k,f(b)都是知道的，把​\begin{cases}k=\lfloor m/i \rfloor \\ b=m \bmod i\end{cases}代入，再把​kf(b)这一项移到右边，得到​i^{-1} \equiv -\lfloor m/i \rfloor f(m \bmod i) \pmod{m}，即​f(i) \equiv -\lfloor m/i \rfloor f(m \bmod i) \pmod{m}。是的，就这样神奇地出来了。为什么？因为我照着OI Wiki抄的我们的猜想很合理。

最终递推式如下：

\begin{cases}f(1) \equiv 1 \\ f(i) \equiv -\lfloor m/i \rfloor f(m \bmod i) \pmod{m} & i>1\end{cases}

代码如下：

void GetInv(){
	inv[1] = 1;
	for (int i = 2; i <= n; ++i) {
  		inv[i] = (long long)(p - p / i) * inv[p % i] % p;
	}
}

Reference 参考

  • 《具体数学》
  • OI Wiki
  • 维基百科

‍