模运算的世界

模运算的世界

逆元

接下来，我们来学习如何求解线性同余方程。让我们考虑如何求解线性方程 $ax\equiv b(modm)$。对于实数运算下的方程 $ax=b$ ，由于 $a$ 存在倒数，因此很容易求解。如果在 $modm$ 的运算下，也有像满足 $ay\equiv 1(modm)$ 这样的 $a$ 的倒数一样的数存在，方程就可以求解了。我们把这样的y叫做 $a$ 的逆元，记作 $a^{-1}$ 。如果能求解逆元，那么就有 $x=a^{-1}\times xa=a^{-1}\times b$ ，也就可以求出 $x$ 了。由于方程 $ax\equiv 1(modm)$ 等价于存在整数 $k$ 使得 $ax=1+mk$ ，因此稍作变形之后，可以变为求解满足 $ax-mk=1$ 的 $x$ 的问题。这个问题可以使用 $extgcd$ 求解。同时，也可以知道如果 $gcd(a，m)!=1$ ，那么逆元是不存在的。

int mod_inverse (int a, int m) {
    int x, y;
    extgcd(a, m, x, y);
    return (m + x % m) % m;
}

如果 $a$ 和 $m$ 不互素，那么$ax\equiv b(modm)$ 就等价于

$(a/gcd(a，m))x\equiv b/gcd(a，m)(modm/gcd(a，m))$

从式子中可以看出，当 $b$ 无法整除 $gcd(a，m)$ 时，原方程无解。在有解的情况下，我们有

$x\equiv (a/gcd(a，m){)}^{-1}\times (b/gcd(a，m))(modm/gcd(a，m))$

因此， $ax\equiv b(modm)$ 的解为

$x\equiv (a/gcd(a，m){)}^{-1}\times (b/gcd(a，m))+(m/gcd(a，m))\times k(modm)(0\le k<gcd(a，m))$

需要注意的是这里和实数的情况有所不同，有可能有多解，也有可能无解。

费马小定理

在 $p$ 是素数的情况下，对任意整数 $x$ 都有 $x^p\equiv x(modp)$ 。这个定理被称作费马小定理。其中如果 $x$ 无法被 $p$ 整除，我们有 $x^{p-1}\equiv 1(modp)$ 。利用这条性质，在 $p$ 是素数的情况下，就很容易求出一个数的逆元。把上面的式子变形之后得到 $a^{-1}\equiv a^{p-2}(modp)$ ，因此可以通过快速幂运算求出逆元。

在不是素数的情况下，我们也有类似的欧拉定理可以使用。假设 $m=p_1^{e_1}{p}_2^{e_2}...p_n^{e_n}$ ，那么 $m$ 的欧拉函数 $\phi (m)$ 的定义如下
$\phi (m)=m\times \Pi (p-1)/p$，（$\Pi$为连乘符号）
欧拉函数的值等于不超过 $m$ 并且和 $n$ 互素的数的个数。此时对于和加互素的 $x$ ，有 $x^{\phi (m)}\equiv 1(modm)$成立。
在m是素数时，根据定义\varphi (m)=m-1，因此欧拉定理也可以看作是费马小定理的推广。

因为 $n$ 的整数分解可以在 $O(\sqrt{n})$ 时间内完成，所以对于某一个欧拉函数也可以在 $O(\sqrt{n})$ 时间内求得。另外，我们可以在利用埃式筛法，每次发现质因子时就把它的倍数的欧拉函数乘上 $(p-1)/p$ ，这样就可以一次性求出 $1n$ 的欧拉函数值的表了。

//求欧拉函数值。O（√n）
#define MAX_N 1000
int euler_phi(int n) {
	int res = n;
	for (int i = 2; i*i <= n; i++) {
		if (n%i == 0) {
			res = res / i * (i - 1);
			for (; n%i == 0; n /= i);
		}
	}
	if (n != 1) res = res / n * (n - 1);
	return res;
}

int euler[MAX_N];

//O（MAX_N）时间筛出欧拉函数值的表
void euler_phi2() {
	for (int i = 0; i < MAX_N; i++)
		euler[i] = i;
	for (int i = 2; i < MAX_N; i++) {
		if (euler[i] == i) {
			for (int j = i; j < MAX_N; j += i)
				euler[j] = euler[j] / i * (i - 1);
		}
	}
}

线性同余方程组

详情请移步：求解线性同余方程组
下面给大家介绍一下如何求解由多条线性同余方程联立得到的线性同余方程组。用数学化的符号表示就是求解 $a_i\times x\equiv b_i(modm)(1\le i\le n)$ 这样的方程组。如果方程组有解，那么一定有无穷多解，而且解的全集一定可以写成 $x\equiv b(modm)$ 的形式，因此问题就转化为求解 $b$ 和 $m$ 。如果我们能求解方程组 $x\equiv b_1(mod{m}_1)$ ， $a\times x\equiv b_2(mod{m}_2)$ ，那么只要对方程逐个求解，对于任意有限的 $n$ 我们就都可以求出答案了。

因为 $x\equiv b(mod{m}_i)$ ，所以可以把 $x$ 写成 $x=b_1+m_1\times t$ 的形式。把它代入第二条式子，可以得到

$a(b_1+m_1\times t)\equiv b_2(mod{m}_2)$

移项整理后得到

$a\times m_1\times t\equiv b_2-a\times b_1(mod{m}_2)$

由于这只是一个一次同余方程，因此很容易求解。当 $gcd(m_2，a\times m_1)$ 无法整除 $b_2-a\times b_1$ 时原方程组无解。

//返回一个（b，m）的数对
pair<int, int>linear_congruence(const vector<int>&A, const vector<int>&B, const vector<int>&M) {
	//由于最开始没有任何限制，所以先把解设为表示所有整教的 x三0（mod 1）
	int x = 0, m = 1;
	for (int i = 0; i < A.size(); i++) {
		int a = A[i] * m, b = B[i] - A[i] * x, d = gcd(M[i], a);
		if (b%d != 0) return make_pair(0, -1);//无解
		int t = b / d * mod_inverse(a / d, M[i] / d) % (M[i] / d);
		x = x + m * t;
		m *= M[i] / d;
	}
	return make_pair(x % m, m);
}

中国剩余定理

详情请移步：中国剩余定理证明
我们假设同余方程组里所有的 $a$ 都等于 $1$ ，并且所有的 $m$ 都互素。这样，答案就一定是 $x\equiv b(mod\Pi m_i)$ 。反之，对于一个合数 $n$ ，我们假设有 $n=ab（其中a和b互素）$ 。那么如果$xmodn$ 的值确定， $xmoda$ 和 $xmodb$ 的值就都确定了。也就是说，我们有 $（xmodn）\iff （xmoda,xmodb）$ 这样一组对应关系。

换句话说，以合数 $n$ 为模数来考虑与以 $a$ 和 $b$ 为模数来考虑是等价的。这个定理叫做中国剩余定理。通过对 $n$ 进行分解，对于模合数的情况只需要考虑模 $p^k（p为素数）$ 的情况就可以了。其中，如果 $n$ 不能被任何一个完全平方数整除，那么问题就可以转化为模数为素数的情况，从而变得容易求解。中国剩余定理不是一个算法，而是可以看成在思考算法时的一个提示

例： $f(x)\equiv 0（modn）\iff f(x)\equiv 0（mod{p}^k）（p^k|n）$

n！（n的阶乘）

在计数问题中，经常需要用到 $n!$ 。在学完前面介绍的内容之后，有必要了解 $n!$ 在 $modp$ 下的一些性质。下面我们假设 $p$ 是素数， $n!=a{p}^e（a无法被p整除）$ ，并试图求解 $amodp$ 和 $e$ 。$e$ 是$n!$ 能够迭代整除 $p$ 的次数，因此可以使用下面的式子进行计算。

$n/p+n/p^2+n/p^3+\dots$

这个结论很显然，因为 $n/d$ 和不超过 $n$ 的能被 $d$ 整除的正整数的个数相等。由于只需要对于 $p^t\le n$ 的 $t$ 进行计算，因此复杂度是 $O（lo{g}_{p}n）$ 。

接下来计算 $amodp$ 。首先计算 $n!=1\times 2\times \dots \times n$ 的因数中不能被 $p$ 整除的项的积。假设 $n=10，p=3$ ，则

$n!=1\times 2\times 4\times 5\times 7\times 8\times 10\times (3\times 6\times 9)$
$1\times 2\times 4\times 5\times 7\times 8\times 10\equiv 1\times 2\times 1\times 2\times 1\times 2\times 1(modp)$

从这个例子中可以看出，不能被 $p$ 整除的项在 $modp$ 下呈周期性。因此，不能被 $p$ 整除的项的积等于 $(p-1){!}^{(n/p)}\times (nmodp)!$ 。事实上，根据威尔逊定理，我们有 $(p-1)!\equiv -1$ 。因为除了 $1$ 和 $p-1$ 之外其余的项都可以和各自的逆元相乘得到1。

接下来，计算可以被 $p$ 整除的项的积。可以被 $p$ 整除的项有 $p,2p,3p,\dots ,(n/p)p$ 。把这些项分别除以 $p$ 之后得到 $1,2,3.\dots ,n/p$ 。因此，问题的范围就由 $n$ 缩小到了 $n/p$ 。如果预处理出 $0\le n<p$ 范围中 $n!modp$ 的表，就可以在 $O(lo{g}_{p}n）$ 时间内算出答案。如果不预处理，那么复杂度是 $O(plo{g}_{p}n)$ 。

//分解n！=a p^e，返回a mod p。 O（1og，n）。
int mod_fact(int n, int p, int& e) {
	e = 0;
	if (n == 0) return 1;

	//计算p的倍数的部分
	int res = mod_fact(n / p, p, e);
	e += n / p;

	//由于（p-1）！=-1，因此（p-1）！^{(n/p)}只需要知道n/p的奇偶性就可以计算了。
	if (n / p % 2 != 0) return res * (p - fact[n % p]) % p; 
	return res * fact[n % p] % p;
}

$C_n^{k}modp$

了解了 $n!$ 在 $modp$ 下的性质之后， $C_n^{k}modp$ 也就可以计算了。首先，把 $C_n^k$ 写成 $n!$ 的积的形式。

$C_n^k=\frac{n!}{k!(n-k)!}$

设 $n!=a_1{p}_1^{e_1}，k!=a_2{p}_2^{e_2}，(n-k)!=a_3{p}_3^{e_3}$ 。从式子中可以看出，当 $e_1>e_2+e_3$ 时， $C_n^k$ 可以被 $p$ 整除， $e_1=e_2+e_3$ 时无法被 $p$ 整除。在无法整除的情况下， $C_n^k=a_1(a_2{a}_3{)}^{-1}$ 。

//求C^k_n mod p。O（1og_p n）。
int mod_comb(int n, int k, int p) {
	if (n < 0 || k < 0 || n < k) return 0;
	int e1, e2, e3;
	int a1 = mod_fact(n, p, e1), a2 = mod_fact(k, p, e2), a3 = mod_fact(n - k, p, e3);
	if (e1 > e2 + e3) return 0;
	return a1 * mod_inverse(a2 * a3 % p, p) % p;
}

扩展中国剩余定理

来自大佬的详解