常用算法模板——数学知识：质数、约数、欧拉函数、快速幂、扩展欧几里得算法、高斯消元、组合数
C++算法模板系列
目录
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 \bmod 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^k\limits_{i=1}(\alpha_i+1)=(\alpha_1+1)(\alpha_2+1)\cdots(\alpha_k+1)约数之和：
\prod^k\limits_{i=1}(\sum^{\alpha_i}\limits_{j=0}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 \bmod 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互质的数的个数称为欧拉函数，记为\phi(N)。
求法：若N=p_1^{a_1}p_2^{a_2}\cdots p_m^{a_m} ，则
\phi(N)=N\prod^m\limits_{i=1}(1-\frac1{p_i})特别地，对于质数p，有
\phi(p)=p-1欧拉定理：若a与m互质，即\gcd(a,m)=1，则
a^{\phi(m)}\equiv 1 \pmod m费马小定理：若p为质数，则
a^{\phi(p)}\equiv1\pmod p \Rightarrow a^{p-1}\equiv 1\pmod p \Rightarrow a^p\equiv a\pmod p3.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 \bmod p，时间复杂度：O(\log k)
思想：
a^k\bmod p=\prod\limits^{\log k}_{i=0}a^{2^i}\bmod p\\ a^{2^{i+1}}\bmod p=(a^{2^i}\bmod p)^2\bmod p核心：求k的二进制表示，即
a^k\bmod p=\prod\limits_{i\in\{i|k[i]=1\}} a^{2^i}\bmod p应用：求模为质数的逆元
逆元的定义：ax\equiv 1\pmod m，a,m互质，则称x为a模m的逆元，记为a^{-1}。
求法：当模为质数p时，由费马小定理得
a^{\phi(p)}\equiv1\pmod p \Rightarrow a^{p-1}\equiv 1\pmod p \Rightarrow a\cdot a^{p-2}=1\pmod 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满足ax_0+b_0y=d ①，其中d=\gcd(a,b)。原方程可化为a(x-\frac bd)+b(y+\frac ad)=d②。由①②可得通解为
\left\{\begin{matrix}
x=x_0-\frac bdk\\ y=y_0+\frac adk,
\end{matrix}\right.
,\ k\in \mathbb{Z^+}应用：求解线性同余方程。对于方程ax\equiv b\pmod m，\exist y\in\mathbb{Z^+},ax=my+b\ \xRightarrow{y'=-y} ax+my'=b，该方程有解的充要条件为\gcd(a,m)|b，此时可用扩展欧几里得算法exgcd(a, b, x, y')求得一组特解，进而求得原方程的解。特别地，求a模m的逆元即为b=1的情况。
中国剩余定理：对于两两互质的k个数m_1,m_2,\cdots,m_k，线性同余方程组
\left\{\begin{matrix}
x\equiv a_1\pmod{m_1}\\
x\equiv a_2\pmod{m_2}\\
\vdots \\
x\equiv a_k\pmod{m_k}
\end{matrix}\right.的通解为
x=a_1M_1M_1^{-1}+a_2M_2M_2^{-1}+\cdots+a_kM_kM_k^{-1}其中M=m_1m_2\cdots m_k，M_i=\frac{M}{m_i}\ (i=1,2,\cdots,k)，M_i^{-1}为M_i模m_i的逆元，可通过解M_ix\equiv1\pmod{m_i}求得。
/* 求一组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; // 这里传回特解(1, 0)
y = 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（或\text{C}(n,m)、{n \choose m}）的定义：
\text{C}_n^m=\frac{n\times(n-1)\times\cdots\times(n-m+1)}{1\times2\times\cdots\times m}=\frac{n!}{(n-m)!m!}\ (m\le n)互补性：
\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≤m≤n≤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≤m≤n≤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≤p≤m≤n，有
\text{C}_n^m \equiv \text{C}_{n \bmod p}^{m \bmod p}\text{C}_{n/p}^{m/p}\pmod p适用于较低数据量、 1≤m≤n≤10^{18},1≤p≤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 np + \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}
《C++算法模板（4）：数学知识》有1条评论