常用算法代码模板 (4) ：数学知识

常用算法模板——数学知识：质数、约数、欧拉函数、快速幂、扩展欧几里得算法、高斯消元、组合数

---

常用算法代码模板 (1) ：基础算法
常用算法代码模板 (2) ：数据结构
常用算法代码模板 (3) ：搜索与图论
常用算法代码模板 (4) ：数学知识
算法选择——由数据范围反推算法时间复杂度

---

---

1 质数

定义：在大于 $1$ 的整数中，只包含 $1$ 和本身这两个约数的数称为质数（或素数）

1.1 试除法判定质数

时间复杂度：$O(\sqrt{n})$

bool is_prime(int x) {
    if (x < 2) return false;
    for (int i = 2; i <= x / i; i++)	// 枚举到sqrt(x)
        if (x % i == 0)
            return false;
    return true;
}

1.2 试除法分解质因数

时间复杂度：$O(\sqrt{n})$

vector<pair<int, int> > primes;	// 存储质因数及其个数

void divide(int x) {
    for (int i = 2; i <= x / i; i++)	// 枚举到sqrt(x)
        if (x % i == 0) {
            int cnt = 0;	// cnt记录质因子i的个数
            while (x % i == 0) {
                x /= i;
                cnt++;
            }
            primes.push_back({i, cnt});
        }
  
    if (x > 1) primes.push_back({x, i});	// 原理：x中只包含1个大于sqrt(x)的质因子
}

1.3 筛法求素数表

1.3.1 埃氏筛法

时间复杂度：$O(n\log \log n)$

int primes[N], len;	// 存储所有素数
bool st[N];			// st[i]标记数i是否被筛掉

/* 埃氏筛法求[2, n]上所有素数 */
void get_primes(int n) {
    for (int i = 2; i <= n; i++) {
        if (!st[i]) {	// 仅遍历未被筛去的数，且只筛它的倍数
            primes[len++] = i;
            for (int j = i + i; j <= n; j += i)	// 筛去i的倍数，朴素法遍历全部倍数
                st[j] = true;
        }
    }
}

1.3.2 线性筛法

时间复杂度：$O(n)$

核心思想：每个合数只会被其最小质因子筛掉。对于 $i$ 和素数 $P_j$ ，若 $i\mathrm{mod}{P}_j=0$ ，且 $P_j$ 是 $i$ 的最小质因子，即一定是 $P_j\cdot i$ 的最小质因子。

int primes[N], len;	// 存储所有素数
bool st[N];			// st[i]标记数i是否被筛掉

/* 线性筛法求[2, n]上所有素数 */
void get_primes(int n) {
    for (int i = 2; i <= n; i++) {
        if (!st[i]) primes[len++] = i;
        for (int j = 0; primes[j] <= n / i; j++) {	// primes[j] * i <= n
            st[primes[j] * i] = true;		// 每个合数只会被其最小质因子筛掉
            if (i % primes[j] == 0) break;	// 保证primes[j]一定是primes[j] * i的最小质因子
        }
    }
}

---

2 约数

2.1 试除法求所有约数

时间复杂度：取决于排序函数，试除的消耗是 $O(\sqrt{n})$

/* 求所有约数（去重且递增排序） */
vector<int> get_divisors(int x) {
    vector<int> res;
    for (int i = 1; i <= x / i; i++)	// 枚举到sqrt(x)
        if (x % i == 0) {	// 若i为x的约数，则x/i也是x的约数
            res.push_back(i);
            if (i != x / i) res.push_back(x / i);	// 不重复存储约数sqrt(x)
        }
    sort(res.begin(), res.end());
    return res;
}

2.2 约数个数、约束之和

设 $N=p_1^{{\alpha }_1}{p}_2^{{\alpha }_2}\cdots p_k^{{\alpha }_k}$ ，其中 $p_i$ 是试除法求得的约数（个数为 ${\alpha }_i$ ），则

• 约数个数：$\prod _{i=1}^k({\alpha }_i+1)=({\alpha }_1+1)({\alpha }_2+1)\cdots ({\alpha }_k+1)$
• 约数之和：$\prod _{i=1}^k(\sum _{j=0}^{{\alpha }_i}{p}_i^j)=(p_1^0+p_1^1+\cdots +p_1^{{\alpha }_1})\cdots (p_k^0+p_k^1+\cdots +p_k^{{\alpha }_k})$

typedef long long LL;

const int MOD = 1e9 + 7;		// 防止结果过大而溢出

int x;
vector<pair<int, int> > primes;	// 用4.1(2) divide()返回的数组，存储约数p_k及其个数a_k

/* 求约数个数 */
LL cnt = 1;
for (auto prime : primes) cnt = cnt * (prime.second + 1) % MOD;

/* 求约数之和 */
LL sum = 1;
for (auto prime: primes) {
    int p = prime.first, a = prime.second;	// 约数p与指数a
    LL t = 1;	// 记录p^0+...+p^a
    while (a--) t = (t * p + 1) % MOD;		// 秦九韶算法
    sum = sum * t % MOD;
}

2.3 欧几里得算法

整除的性质：若 $d\mid a$ ，$d\mid b$ ，则 $d\mid ax+by$

欧几里得算法（辗转相除法）求最大公约数：$\gcd (a,b)=\gcd (b,a\mathrm{mod}b)$

时间复杂度：$O(\log n)$

int gcd(int a, int b) {
    return b ? gcd(b, a % b) : a;	// gcd(a, 0) = a
}

2.4 求最小公倍数

$\text{lcm}(a,b)=\frac{ab}{\gcd (a,b)}$

int lcm(int a, int b) {
	return a / gcd(a, b) * b;
}

---

3 欧拉函数

定义：$1$ ~ $N$ 中与 $N$ 互质的数的个数称为欧拉函数，记为 $\varphi (N)$ 。

求法：若$N=p_1^{a_1}{p}_2^{a_2}\cdots p_m^{a_m}$ ，则 $\varphi (N)=N\prod _{i=1}^m(1-\frac{1}{p_i})$ 。特别地，对于质数 $p$ ，$\varphi (p)=p-1$ 。

欧拉定理：若 $a$ 与 $m$ 互质，即 $\gcd (a,m)=1$ ，则 $a^{\varphi (m)}\equiv 1(\mathrm{mod}m)$ 。

费马小定理：若 $p$ 为质数，则 $a^{\varphi (p)}\equiv 1(\mathrm{mod}p)\Rightarrow a^{p-1}\equiv 1(\mathrm{mod}p)\Rightarrow a^p\equiv a(\mathrm{mod}p)$ 。

3.1 求欧拉函数

时间复杂度：$O(\sqrt{n})$

int phi(int x) {
	int res = x;
	for (int i = 2; i <= x / i; i++)	// 试除法分解质因数
		if (x % i == 0) {
			res = res / i * (i - 1);	// 化简(1 - 1 / i)所得
			while (x % i == 0) x /= i;
		}

	if (x > 1) res = res / x * (x - 1);
	return res;
}

3.2 筛法求欧拉函数表

时间复杂度：$O(n)$

int primes[N], len;	// 存储所有素数
int euler[N];		// euler[x]存储x的欧拉函数
bool st[N];			// st[x]存储x是否被筛掉

/* 线性筛法求[1, n]上所有数的欧拉函数 */
void get_eulers(int n) {
	euler[1] = 1;	// 规定1与任何数互质
	for (int i = 2; i <= n; i++) {
		if (!st[i]) {
			primes[len++] = i;
			euler[i] = i - 1;	// 若i为质数，则phi(i)=i-1（1~i-1均与i互质）
		}
		for (int j = 0; primes[j] <= n / i; j++) {	// p_j * i <= n
			st[primes[j] * i] = true;
			if (i % primes[j] == 0) {	// 若p_j是i的最小质因子，则一定是p_j * i的最小质因子
				euler[primes[j] * i] = euler[i] * primes[j];	// 因此phi(p_j * i)的\prod部分与phi(i)完全相同
				break;
			}
			euler[primes[j] * i] = euler[i] * (primes[j] - 1);	// 否则phi(p_j * i) = p_j * phi(i) * (1 - 1 / p_j)
		}
	}
}

---

4 快速幂

求 $a^k\mathrm{mod}p$，时间复杂度：$O(\log k)$

思想：$a^k\mathrm{mod}p=\prod _{i=0}^{\log k}{a}^{2^i}\mathrm{mod}p$ ，$a^{2^{i+1}}\mathrm{mod}p=(a^{2^i}\mathrm{mod}p{)}^2\mathrm{mod}p$ 。

核心：求 $k$ 的二进制表示，即 a k   m o d   p = ∏ i ∈ { i ∣ k [ i ] = 1 } a 2 i   m o d   p a^k\bmod p=\prod\limits_{i\in\{i|k[i]=1\}} a^{2^i}\bmod p akmodp=i∈{i∣k[i]=1}∏​a2imodp 。

应用：求模为质数的逆元（逆元的定义：$ax\equiv 1(\mathrm{mod}m)$ ，$a,m$ 互质，则称 $x$ 为 $a$ 模 $m$ 的逆元，记为 $a^{-1}$ ）。当模为质数 $p$ 时，由费马小定理得 $a^{\varphi (p)}\equiv 1(\mathrm{mod}p)\Rightarrow a^{p-1}\equiv 1(\mathrm{mod}p)\Rightarrow a\cdot a^{p-2}=1(\mathrm{mod}p)$，故此时 $a^{-1}=a^{p-2}$ ，可用快速幂求得。

typedef long long LL;

/* 求a^k mod p */
LL qpow(int a, int k, int p) {
	LL res = 1, t = a;	// t记录a^2^i，其中i>=0，表示逻辑上当前迭代至k的第i位
	while (k) {
		if (k & 1) res = res * t % p;	// 当k末位(k[i])为1时，结果乘上a^2^i mod p
		t = t * t % p;	// 更新操作 t <- a^2^(i+1) mod p = (a^2^i mod p)^2 mod p
		k >>= 1;		// k去掉当前末位，使得逻辑上i++
	}
	return res;
}

---

5 扩展欧几里得算法

裴蜀定理：对于任意正整数 $a,b$ ，存在非零整数 $x,y$ ，使得 $ax+by=\gcd (a,b)$ 。

求通解：设特解 $x_0,y_0$ 满足 $a{x}_0+b_{0}y=d$ ①，其中 $d=\gcd (a,b)$ 。原方程可化为 $a(x-\frac{b}{d})+b(y+\frac{a}{d})=d$ ②。由①②可得通解为 $x=x_0-\frac{b}{d}k,y=y_0+\frac{a}{d}k,k\in {\mathbb{Z}}^{+}$ 。

应用：求解线性同余方程。对于方程 $ax\equiv b(\mathrm{mod}m)$ ，$\exists y\in {\mathbb{Z}}^{+},ax=my+b\stackrel{y^{\prime }=-y}{\Rightarrow }ax+m{y}^{\prime }=b$，该方程有解的充要条件为 $\gcd (a,m)|b$ ，此时可用扩展欧几里得算法exgcd(a, b, x, y')求得一组特解，进而求得原方程的解。特别地，求 $a$ 模 $m$ 的逆元即为 $b=1$ 的情况。

中国剩余定理：对于两两互质的 $k$ 个数 $m_1,m_2,\cdots ,m_k$ ，线性同余方程组 $x\equiv a_i(\mathrm{mod}{m}_i)$ 的通解为 $x=a_1{M}_1{M}_1^{-1}+a_2{M}_2{M}_2^{-1}+\cdots +a_k{M}_k{M}_k^{-1}$，其中 $M=m_1{m}_2\cdots m_k,M_i=\frac{M}{m_i},$ $M_i^{-1}$ 为 $M_i$ 模 $m_i$ 的逆元（可通过解 $M_{i}x\equiv 1(\mathrm{mod}{m}_i)$ 求得）$,i=1,2,\cdots ,k$ 。

/* 求一组x, y特解，满足a*x + b*y = gcd(a, b)。函数返回最大公约数 */
int exgcd(int a, int b, int &x, int &y) {	// x, y用引用型
	if (!b) {	// gcd(a, 0) = a，此时a*x + 0*y = a的通解为(1, 任何数)
		x = 1, y = 0;	// 这里传回特解(1, 0)
		return a;
	}
	int d = exgcd(b, a % b, y, x);	// 将x、y翻转(便于对比条件)，使得b*y + (a mod b)*x = gcd(b, a mod b) = gcd(a, b)
	y -= (a / b) * x;	// a mod b = a - [a/b]*b，代入上式化简得a*x + b*(y - [a/b]*x) = d，对比条件可知y的变化量！
	return d;
}

---

6 高斯消元

化增广矩阵为最简行阶梯形矩阵，解 $n$ 元线性方程组，时间复杂度：$O(n^3)$

const double eps = 1e-6;

int n;
double a[N][N];	// n*(n+1)的增广矩阵 a[0 ... n-1][0 ... n]

/* 0：有唯一解（此时将增广矩阵化为最简行阶梯形矩阵），1：有无穷多组解，2：无解 */
int gauss() {
	int c, r;	// 枚举的列、行（同时也记录实际方程个数）
	for (c = 0, r = 0; c < n; c++) {	// 枚举每一列c，最终化为行阶梯形矩阵
		int t = r;
		for (int i = r; i < n; i++)	// 寻找绝对值最大的行，记录于t
			if (fabs(a[i][c]) > fabs(a[t][c]))
				t = i;
				
		if (fabs(a[t][c]) < eps) continue;	// 若为0则无需消元，跳至下一列

		for (int i = c; i <= n; i++) swap(a[t][i], a[r][i]);	// 将绝对值最大的行t换到最顶端的当前行r
		for (int i = n; i >= c; i--) a[r][i] /= a[r][c];	// 将当前行r同除以该行首a[r][c]，使得行r首非零元变成1
		for (int i = r + 1; i < n; i++)	// 用当前行r将行r首非零元该列下方所有元素消成0
			if (fabs(a[i][c]) > eps)	// 若为0则无需再遍历操作该行，节省时间
				for (int j = n; j >= c; j--)	// 行i同减行r各列同列元a[r][j]乘以行i首非零元a[i][c]（为了消之为0）
					a[i][j] -= a[r][j] * a[i][c];

		r++;	// 完成本列c消元操作后才跳至下一行
	}

	if (r < n) {	// 若化简后的方程个数小于n（最后n-r行系数阵部分全为0，即0行），则有无穷多组解或无解
		for (int i = r; i < n; i++)	// 若某行左边为0而右边非0，则直接判无解
			if (fabs(a[i][n]) > eps)
				return 2;	// 无解
		return 1;	// 有无穷多组解
	}

	for (int i = n - 1; i >= 0; i--)	// 有唯一解则化为最简行阶梯形矩阵，列n所存的即为解
		for (int j = i + 1; j < n; j++)	// 同化行阶梯形操作，用各行首非零元将其列上上全部元素消为0
			a[i][n] -= a[i][j] * a[j][n];
	return 0;	// 有唯一解
}

---

7 组合数

组合数的定义： ${\text{C}}_n^m=\frac{n\times (n-1)\times \cdots \times (n-m+1)}{1\times 2\times \cdots \times m}=\frac{n!}{(n-m)!m!}(m\le n)$ ，或写作 $\text{C}(n,m)$、$\left(\genfrac{}{}{0ex}{}{n}{m}\right)$

互补性：${\text{C}}_n^m={\text{C}}_n^{n-m}$

递推式：${\text{C}}_n^m={\text{C}}_{n-1}^m+{\text{C}}_{n-1}^{m-1}$

7.1 递推法求组合数

适用于处理 $10^5$ 量级的数据量、 $1\le m\le n\le 2000$ 的情况，时间复杂度：$O(n^2)$

const int MOD = 1e9 + 7;	// 防止结果过大溢出

int c[N][N];	// c[i][j]即为C(i, j)，表示从i个不同元素中取j个的方案数

/* 计算C(0, 0) ~ C(N-1, N-1) */
for (int i = 0; i < N; i++)
	for (int j = 0; j <= i; j++)
		if (!j) c[i][j] = 1;	// 规定“取0个/不取”算作只有1种方案
		else c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % MOD;

7.2 通过逆处理逆元的方式求组合数

${\text{C}}_n^m=n!\cdot ((n-m)!{)}^{-1}\cdot (m!{)}^{-1}$ ，其中$((n-m)!{)}^{-1},(m!{)}^{-1}$ 分别为 $(n-m)!,m!$ 模 $p$ 的逆元，$p$ 为质数。由费马小定理可知逆元 $a^{-1}=a^{p-2}$ 。

适用于处理 $10^4$ 量级的数据量、$1\le m\le n\le 10^5$ 的情况，时间复杂度：$O(n\log n)$

typedef long long LL;		// 预处理时临时防爆int

const int MOD = 1e9 + 7;	// 防止结果过大溢出

int fact[N];	// fact[i]存储i的阶乘再取模
int infact[N];	// infact[i]存储fact[i]的模为质数的逆元再取模

/* 快速幂模板，用来求模为质数的逆元 */
LL qpow(int a, int k, int p) {
	LL res = 1, t = a;
	while (k) {
		if (k & 1) res = res * t % p;
		t = t * t % p;
		k >>= 1;
	}
	return res;
}

/* 预处理阶乘的余数fact[]和阶乘逆元的余数infact[] */
fact[0] = infact[0] = 1;	// 0! = 1，1的任何逆元为1
for (int i = 1; i < N; i++) {	// 递推求解
	fact[i] = (LL)fact[i - 1] * i % MOD;
	infact[i] = (LL)infact[i - 1] * qpow(i, MOD - 2, MOD) % MOD;
}

/* C(a, b)的值 */
int C(int a, int b) {
	return fact[a] * infact[a - b] * infact[b] % MOD;
}

7.3 Lucas定理

Lucas定理：若 $p$ 为质数，则对于任意整数 $1\le p\le m\le n$ ，有 ${\text{C}}_n^m\equiv {\text{C}}_{n\mathrm{mod}p}^{m\mathrm{mod}p}{\text{C}}_{n/p}^{m/p}(\mathrm{mod}p)$

适用于较低数据量、 $1\le m\le n\le 10^{18},1\le p\le 10^5$ 的情况，时间复杂度：$O({\log }_{p}n\cdot p\log p)$

typedef long long LL;

/* 快速幂模板，用来求模为质数的逆元 */
LL qpow(int a, int k, int p) {
	LL res = 1, t = a;
	while (k) {
		if (k & 1) res = res * t % p;
		t = t * t % p;
		k >>= 1;
	}
	return res;
}

/* 求int型数的组合数C(a, b) */
int C(int a, int b, int p) {
	if (a < b) return 0;

	LL x = 1, y = 1;	// 由最初的定义式，x为分子，y为分母
	for (int i = a, j = 1; j <= b; i--, j++) {
		x = x * i % p;	// x = a! / (a-b)!
		y = y * j % p;	// y = b!
	}

	return x * qpow(y, p - 2, p) % p;	// 结果即为x * (y的逆元) mod p
}

/* 通过Lucas定理求long long型数的组合数lucas(a, b) */
int lucas(LL a, LL b, int p) {
	if (a < p && b < p) return C(a, b, p);	// a, b都小于p时用逆元法即可
	return (LL)C(a % p, b % p, p) * lucas(a / p, b / p, p) % p;
}

7.4 分解质因数法求组合数

不取模，求出组合数的真实值。

1. 筛法求出范围内的所有质数
2. 通过公式 ${\text{C}}_n^m=\frac{n!}{(n-m)!m!}$ 求出每个质因子的次数。$n!$ 中 $p$ 的次数是 $\frac{n}{p}+\frac{n}{p^2}+\frac{n}{p^3}+\cdots$
3. 用高精度乘法将所有质因子相乘

int n;
int primes[N], len;	// 存储所有质数
int sum[N];			// 存储每个质数的次数
bool st[N];			// 存储每个数是否已被筛掉

/* 线性筛法求素数 */
void get_primes(int n) {
	for (int i = 2; i <= n; i++) {
		if (!st[i]) primes[cnt++] = i;
		for (int j = 0; primes[j] <= n / i; j++) {
			st[primes[j] * i] = true;
			if (i % primes[j] == 0) break;
		}
	}
}

/* 求n!中质因子p的次数：n / p + n / p^2 + n / p^3 + ... */
int get(int n, int p) {
	int res = 0;
	while (n) {
		res += n / p;
		n /= p;
	}
	return res;
}

/* 高精度乘低精度模板 */
vector<int> mul(vector<int> &A, int b) {
    vector<int> C;
    int t = 0;
    for (int i = 0; i < A.size() || t; i++) {
        if (i < A.size()) t += A[i] * b;
        C.push_back(t % 10);
        t /= 10;
    }
  
    // if (t) C.push_back(t);
    while (C.size() > 1 && C.back() == 0) C.pop_back();
    return C;
}

/* 预处理范围n以内所有质因子，存储每个质因子的个数 */
int a, b;	// 所求组合数C(a, b)的上下标

get_primes(a);
for (int i = 0; i < len; i++) {	// 求每个质因子的次数
	int p = primes[i];
	sum[i] = get(a, p) - get(a - b, p) - get(b, p);	// C(a, b) = a! / ((a - b)! * b!)
}

/* 用高精度乘法将所有质因子相乘，存于数组 */
vector<int> res;
res.push_back(1);	// 初值为1
for (int i = 0; i < len; i++)
	for (int j = 0; j < sum[i]; j++)	// 乘sum[i]次primes[i]
		res = mul(res, primes[i]);

7.5 卡特兰数

卡特兰数：$C_n=\frac{{\text{C}}_{2n}^n}{n+1}$