高斯消元


高斯方程的作用

O(n3)O_{(n^3)} 的时间复杂度内解决一个含有n个未知量的线性方程组

\left\{ \begin{array}{**lr**} \ a_{1 1}\,x_1 + a_{1 2}\,x_2 + \; ... \; +a_{1 n}\,x_n = b_1 , \\ \ a_{2 1}\,x_1 + a_{2 2}\,x_2 + \; ... \; +a_{2 n}\,x_n = b_2 , \\ \ ..........\\ \ a_{n 1}\,x_1 + a_{n 2}\,x_2 + \; ... \; +a_{n n}\,x_n = b_n , \end{array} \right.

方程解的情况:无解, 无穷多解, 唯一解

高斯消元的步骤

将系数矩阵 (n行,n+1列) 转化为 “对角线矩阵” 的过程

为了出于精度考虑第一次消的时候要用系数绝对值最大的数置为零后去消其他行的数

核心代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
int gauss(){
int c, r;
for(c = 0, r = 0; c < n; c++){
int t = r;
for(int i = r; i < n; i++){
if((fabs(a[i][c]) > fabs(a[t][c]))) t = i;
}

if(fabs(a[t][c]) < eps) continue;

for(int i = c; i < n+1; i++)swap(a[t][i], a[r][i]);
for(int i = n; i >= c; i--) a[r][i] /= a[r][c];
for(int i = r+1; i < n; i++){
if(fabs(a[i][c]) < eps) continue;
for(int j = n; 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] > eps)) return 2; //无解
}
return 1; //无穷多解
}

for(int i = n-1; i >= 0; i--){
for(int j = i+1; j < n; j++){
a[i][n] -= a[j][n]*a[i][j];
}
}

return 0;//唯一解
}

这个算法的变量和遍历真的是相当多,需要认真理解每一个变量的含义与作用,建议上手多敲几遍。

组合数


组合数 I

原理: Cab=Ca1b+Ca1b1C_a^b = C_{a-1}^{b} + C_{a-1}^{b-1}

解释:CabC_a^b 表示从 aa 堆物品中选出 bb 间物品,这时我们假定从这 aa 件物品中取出一件标记为特殊,那么就可将所有的情况分为选这件特殊的物品与不选,这时

(1) 不选特殊则意味着要从其他 a1a-1 件物品中选出 bb 件,即 Ca1bC_{a-1}^b ,

(2) 选择则意味着要从其他 a1a-1 件物品中选出 b1b-1 件,即 Ca1bC_{a-1}^{b}

实现:

预处理每一个组合数, 再由上面的等式就可得到递推公式即可转换为核心代码

核心代码:

1
2
3
4
5
6
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] + c[i-1][j-1])%mod;
}
}

复杂度:两层循环显然 O(N2)O_{(N^2)} (N 为所求数据a, b的范围)

组合数 II

原理: $C_a^b; =; \frac{a!}{(a-b)! , *, b!} $ 和 除法逆元

记 $ fact[ i ] = i! \ % \ mod$ , infact[i]infact[i]fact[i]fact[i] 关于 modmod 的逆元

Cab  =  fact[a]  ×  infact[ab]  ×  infact[b]C_a^b \;=\; fact[a] \;\times\;infact[a-b]\;\times\;infact[b]

由于 1e9+7 是质数所以通过费马小定理用快速幂求 infact[i]infact[i] 即可

实现:

预处理出给定范围内的所有 fact[i]fact[i]infact[i]infact[i] ,然后结合上述公式即可得到

前置知识:阶乘的逆等于逆的阶乘(易证)

核心代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
fact[0] = infact[0] = 1;
for(int i = 1; i < N; i++){
fact[i] = (LL)fact[i-1] * i % mod;
infact[i] = (LL)infact[i-1] * qmi(i, mod-2, mod) % mod;
}

int n;
cin>>n;

while(n--){
int a, b;
scanf("%d%d", &a, &b);
printf("%d\n", (LL)fact[a] * infact[a-b] % mod * infact[b] % mod);
}

复杂度:预处理遍历 NN,求幂 logN\log{N} 故复杂度为 O(NlogN)O_{(N\log{N})} (N为给定的数据a, b范围)

组合数 III

原理:卢卡斯(Lucas)定理 Cab  Ca  mod  pb  mod  p × Ca/pb/p  (mod  p)C_a^b \ \equiv \ C_{a \;mod\; p}^{b\;mod\;p} \ \times \ C_{a \,/\,p}^{b\,/\,p} \ \ (mod \ \ p) 这里的除法为整除即下取整

证明可参考(简单了解):卢卡斯定理

复杂度:

Ca  mod  pb  mod  pC_{a \;mod\; p}^{b\;mod\;p}采用预处理阶乘和阶乘的逆的方法复杂度为 O(plogp)O_{(p\log{p})} ,对 Ca/pb/pC_{a \,/\,p}^{b\,/\,p} 由于除数结果可能会大于p所以每次得到的值都再次运用卢卡斯定理(类似于对a,b同时求p的进制表示)复杂度为 O(logpn)O_{(\log_{p}{n})} (n为a b的大小)故最终的结果为 O(p logpn logp)O_{(p\ \log_p{n}\ \log{p})} 注意到 logpn\log_p{n}plogpp\log{p} 的大小与 p 的大小具有相反的单调性故此时可以根据大概估计复杂度即可(如忽略 logpn\log_p{n}​ 的影响)

核心代码:

1
2
3
4
5
6
7
8
9
10
11
12
int C(int a, int b){
LL ans = 1;
for(int i = 1, j = a; i <= b; i++, j--){
ans = ans * j % p;
ans = ans * qmi(i, p-2) % p;
}
return (int)ans;
}
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;
}

组合数 IV

求出 CabC_a ^b 的具体值,结果很大需要用高精度表示

原理: $C_a^b; =; \frac{a!}{(a-b)! , *, b!} $ 提取公因数 和 高精度乘法

实现:先对 a!a!(ab)!(a-b)!b!b! 提取所有公因数及其项数,其中一种做法是类似于埃式筛法对阶乘提取所有公因数(求出1~n中所有 p1,p2,p3, ... ,pkp^1, p^2,p^3,\ ...\ ,p^k 的倍数并相加)然后将每个阶乘提取后公因数项数相减即得到组合数的最终分解公因数的结果

然后利用高精度乘法求将分解后的质因数相乘即可

简单来说就是利用分解质因数将乘除法转换为加减法和利用阶乘的规律性简化对阶乘分解质因数的过程

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#include<iostream>
#include<cstdio>
#include<vector>
using namespace std;

const int N = 5010;
int primes[N], cnt;
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;
}
}
}

int get(int n, int p){
int ans = 0;

while(n){
ans += n/p;
n /= p;
}
return ans;
}

vector<int> mul(vector<int> a, int b){
vector<int> c;
int t = 0, len_a = a.size();

for(int i = 0; i < len_a; i++){
t += a[i]*b;
c.push_back(t%10);
t /= 10;
}

while(t){
c.push_back(t%10);
t /= 10;
}

return c;
}

int main(){
int a, b;
cin>>a>>b;

get_primes(a);

for(int i = 0; i < cnt; i++){
int p = primes[i];
sum[i] = get(a, p) - get(b, p) - get(a-b, p);
}


vector<int> ans;
ans.push_back(1);

for(int i = 0; i < cnt; i++){
for(int j = 0; j < sum[i]; j++){
ans = mul(ans, primes[i]);
}
}

for(int i = ans.size()-1; ~i; i--) printf("%d", ans[i]);
puts("");

return 0;
}

卡特兰定理

描述:给定 n 个 0 和 n 个 1,它们将按照某种顺序排成长度为 2n 的序列,求它们能排列成的所有序列中,能够满足任意前缀序列中 0 的个数都不少于 1 的个数的序列有多少个。例如3个零一对可以组成 “000111, 001101, 001011, 010011, 010101“共5种。

转换:将n组10对转换组合问题转换为从(0,0)点向上向右走步到达(n, n)点的问题,0代表向右走步 1代表向上走步,所有合格的走法不能越过对角线(浅绿线)即路线(黄线)不能与红线相交。

对于所有的非法路线,找到第一个与红线相交的点然后将其交点后的路线对于红线做镜像对称(翻折)得到的终点一定会在(n-1, n+1)上,且显然这样一条经过翻折的路径与原路径是一一对应的并且是全包含的,故对应这样一条非法的路径可以用从(0, 0)到(n-1, n+1)的所有路径代替表示

故合法的走法会等于从(0, 0)到(n,n)的所有走法减去从(0, 0)到(n-1, n+1)的所有走法

ans=C2nnC2nn1ans = C_{2n}^n - C_{2n}^{n-1}

$ = \frac{(2n)!}{n! \times n!}\ \ - \ \ \frac{(2n)!}{(n-1)! \times (n+1)!} $

$= \frac{(2n)! \times(n+1) - (2n)!\times(n) }{n!\times(n+1)!} $

$= \frac{(2n)!}{n!\times(n+1)!} $

=C2nnn+1= \frac{C_{2n}^{n}}{n+1}

这个最终结果就被称为卡特兰数

实现:阶乘和阶乘逆元法求组合数

核心代码:

1
2
3
4
5
6
7
8
9
int n;
cin>>n;
int a = 2*n, b = n, ans = 1;

for(int i = 0; i < n; i++){
ans = (LL)ans * (a-i) % mod * qmi(b-i, mod-2, mod) % mod;
}

ans = (LL)ans * qmi(n+1, mod-2, mod) % mod;