-
高斯消元
- 求解含有n个未知数，n个方程的多元线性方程组 O(n^3)
- 初等行变换：
- 某行乘以一个非零数
- 交换两行
- 某行加上另一行的若干倍
- 利用初等行变换将方程组化为上三角矩阵
- 解的情况：
- 完美阶梯型：唯一解
- 非完美阶梯型：
- 0 == 非0：无解
- 0 == 0：无穷解
- 步骤:
- 枚举每一列
- 找到这一列系数的绝对值最大的一行
- 将这一行与第一行交换
- 将改行的第一个数变成一(方程两边同乘某数)
- 把下面所有行的当前列的系数消成0(某行加上第一行的若干倍)
- 枚举每一列
- 代码：
const int N = 110; const double esp = 1e-6; //x < esp,则x = 0,否则 x != 0,由于浮点数精度问题,0可能是0.000001 double a[N][N]; int n; int gauss(){ int r, c; for(r = 1, c = 1; c <= n; ++ c){//遍历每一列 int t = r; for(int i = r; i <= n; ++ i) //找当前列的系数最大的行 if(fabs(a[t][c]) < fabs(a[i][c])) t = i; //记录最大行 if(fabs(a[t][c]) < esp) continue; //最大行系数为0说明该列系数均为0 for(int i = c; i <= n + 1; ++ i) swap(a[r][i], a[t][i]);//否则交换第一行与系数最大行 for(int i = n + 1; i >= c; -- i) a[r][i] /= a[r][c]; //将第一行当前列系数变为1 for(int i = r + 1; i <= n; ++ i) //将下面所有行的当前列的系数变为0 if(fabs(a[i][c]) > esp) //当前列系数非0 for(int j = n + 1; j >= c; -- j) a[i][j] -= a[r][j] * a[i][c]; r ++ ; } if(r <= n){ for(int i = r; i <= n; ++ i) if(fabs(a[i][n + 1]) > esp) return 0; //无解 return 1; //无数解 } for(int i = n; i >= 1; -- i) //从后往前求解 for(int j = i + 1; j <= n; ++ j) a[i][n + 1] -= a[i][j] * a[j][n + 1]; return 2; //唯一解 } int main(){ cin >> n; for(int i = 1; i <= n; ++ i) for(int j = 1; j <= n + 1; ++ j) cin >> a[i][j]; int t = gauss(); if(t == 1) cout << "Infinite group solutions"; else if(t == 0) cout << "No solution"; else if(t == 2){ for(int i = 1; i <= n; ++ i) printf("%.2f\n", a[i][n + 1]); } return 0; }
-
求组合数
- \(C_a^b = \frac{a(a-1)(a-2)···(a-b+1)}{b!}=\frac{a!}{b!(a-b)!}\)
- 一般直接用公式求组合数容易超时
- 故可以通过预处理的方式，用递推式求出所有可能的组合数
- 递推式：\(C_a^b = C_{a-1}^{b-1} + C_{a-1}^{b}\) \(O(N^2)\)
- 当a，b <= 2e3，询问n = 1e4次时
- 代码：
//询问10000次 //a,b <= 2000 const int N = 2e3 + 10; const int mod = 1e9 + 7; int c[N][N]; int n; void init(){ for(int i = 0; i < N; ++ i) for(int j = 0; j <= i; ++ j) if(!j) c[i][j] = 1; else c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % mod; } - 当a，b <= 1e5，询问n = 1e4次，模数p = 1e9 + 7时
- 开二维数组会爆内存，于是利用公式，预处理出阶乘以及阶乘的逆\(O(Nlog(N))\)
- 代码：
typedef long long LL; const int N = 1e5 + 10; const int mod = 1e9 + 7; int fact[N], infact[N]; int n; //快速幂 int qmi(int a, int b, int p){ int res = 1; while(b){ if(b & 1) res = (LL)res * a % p; a = (LL)a * a % p; b >>= 1; } return res; } //预处理各阶乘以及阶乘的逆 void init(){ fact[0] = infact[0] = 1; for(int i = 1; i < N; ++ i){ fact[i] = (LL)fact[i - 1] * i % mod; //递推式n! = (n - 1)! * n; infact[i] = (LL)infact[i - 1] * qmi(i, mod - 2, mod) % mod; //递推式(n!)^-1 = ((n - 1)!)^-1 * n^-1，取模意义下的逆就是逆元 } } int main(){ cin >> n; init(); while(n -- ){ int a, b; cin >> a >> b; cout << (LL)fact[a] * infact[b] % mod * infact[a - b] % mod << endl;//两个int最大值相乘不会爆longlong，但是三个相乘会爆 } return 0; }
- 当a，b <= 1e18，询问n = 2e5次，模数p不定时
- 两数相乘会爆longlong，利用卢卡斯定理预处理组合数 \(O(plog(N)log(p))\)
- 定理：\(C_a^b\equiv C_{a \mod p}^{b \mod p}·C_{a/p}^{b/p} \pmod p\)
- 代码：
typedef long long LL; int p; int n; //快速幂求逆元 int qmi(int a, int b){ int res = 1; while(b){ if(b & 1) res = (LL)res * a % p; a = (LL)a * a % p; b >>= 1; } return res; } int C(int a, int b){ //求组合数C_a^b if(a < b) return 0; //可以不要 int x = 1, y = 1; //分子分母 for(int i = 0; i < b; ++ i){ x = (LL)x * (a - i) % p; //求分子 y = (LL)y * (i + 1) % p; //求分母 } return (LL)x * qmi(y, p - 2) % p; } //卢卡斯定理 int lucas(LL a, LL b){ if(a < p && b < p) return C(a, b); //不必用卢卡斯，可以直接利用公式求出 return (LL)C(a % p, b % p) * lucas(a / p, b / p) % p; //用卢卡斯定理 } int main(){ cin >> n; while(n -- ){ LL a, b; cin >> a >> b >> p; cout << lucas(a, b) << endl; } return 0; }
- 当a，b <= 5e5，询问一次，不取模时
- 此时组合数较大，需要用高精度存储结果
- 阶乘中素数的个数：
- \(v_p(n!) = \sum_{i=1}^{\infty} \left\lfloor \frac{n}{p^i} \right\rfloor\)
- 证明：将 \(n!\) 记为 \(1\times 2\times \cdots \times p\times \cdots \times 2p\times \cdots \times \lfloor n/p\rfloor p\times \cdots \times n\) 那么其中 \(p\) 的倍数有 \(p\times 2p\times \cdots \times \lfloor n/p\rfloor p=p^{\lfloor n/p\rfloor }\lfloor n/p\rfloor !\) 然后在 \(\lfloor n/p\rfloor !\) 中继续寻找 \(p\) 的倍数即可，这是一个递归的过程
- 代码：
const int N = 5e3 + 10; int primes[N], cnt; int st[N]; int sum[N]; //线性筛，1-n中质数就是n!中的质因子 void get_primes(int a){ for(int i = 2; i <= a; ++ i){ if(!st[i]) primes[ ++ cnt] = i; for(int j = 1; primes[j] <= a / i; ++ j){ st[primes[j] * i] = 1; if(i % primes[j] == 0) break; } } } //求a!中质数p的个数 int get(int a, int p){ int res = 0; while(a){ res += a / p; a /= 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; } while(c.size() > 1 && c.back() == 0) c.pop_back(); return c; } int main(){ int a, b; cin >> a >> b; get_primes(a); for(int i = 1; i <= cnt; ++ i){ int p = primes[i]; sum[i] = get(a, p) - get(b, p) - get(a - b, p); //a!中质因子个数-b!中质因子个数-(a-b)!中质因子个数 } vector<int> res; res.push_back(1); for(int i = 1; i <= cnt; ++ i) //累乘 for(int j = 1; j <= sum[i]; ++ j) res = mul(res, primes[i]); for(int i = res.size() - 1; i >= 0; -- i)//输出 cout << res[i]; return 0; }
-
卡特兰数
- \(C_{2n}^n-C_{2n}^{n-1}=\frac{1}{n+1}·C_{2n}^n\)
- 用于解决序列排列等问题
- AcWing 889. 满足条件的01序列 - AcWing
- 代码：
#include <iostream> using namespace std; typedef long long LL; const int N = 1e5 + 10; const int p = 1e9 + 7; int n; int qmi(int a, int k){ int res = 1; while(k){ if(k & 1) res = (LL)res * a % p; a = (LL)a * a % p; k >>= 1; } return res; } int main(){ cin >> n; int res = 1; int a = 2 * n, b = n; int x = 1, y = 1; for(int i = 0; i < b; ++ i){ x = (LL)x * (a - i) % p; y = (LL)y * (i + 1) % p; } res = (LL)x * qmi(y, p - 2) % p; res = (LL)res * qmi(n + 1, p - 2) % p; cout << res; return 0; }
文章来源地址https://www.toymoban.com/news/detail-638449.html
文章来源：https://www.toymoban.com/news/detail-638449.html
到了这里，关于[数论第三节]高斯消元法/求组合数/卡特兰数的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
