渣渣本渣得知自己研究生研究方向为计算成像，而且自己可能会主要负责编程及算法实现这块的内容，于是将重心移至矩阵这块（因为老板说这块以后要用到很多），那这次先系统学习一下关于矩阵的基础常用算法
高斯消元 & 高斯-约旦消元
高斯消元
对于一个线性方程组，我们可以把它写成增广矩阵的形式（就是把等号右侧的列向量拼到系数矩阵右侧）
然后对增广矩阵不断进行初等行变换（只能进行一种变换），将其化为上三角矩阵，最后可以得到，矩阵最后一行最后一个元素的值就是对应的最后一个未知数的解，把这个解代入倒数第二行，得到倒数第二个未知数的解，以此类推，可以得到方程的所有解
时间复杂度：$O(n^3)$
举个栗子：
| x - 2y + 3z = 6
4x - 5y + 6z = 12
7x - 8y + 10z = 21
|
写成增广矩阵为：
| 1 -2 3 | 6
4 -5 6 | 12
7 -8 10 | 21
|
首先把他化为上三角矩阵，用第一行消掉之后所有行的第一个元素
PS：往往是将一列中最大的一个元素所对应的行作为消去使用的行（因为这样会减少误差，一些意会的证明见：P3389 【模板】高斯消元 题解）所以这里把第三行和第一行交换，然后用第一行消，为了写程序方便，先将首元素化为1
| 1 -8/7 10/7 | 3
0 -3/7 -2/7 | 0
0 -6/7 11/7 | 3
|
以此类推，可以得到一个上三角阵：
| 1 -8/7 10/7 | 3
0 1 -2/3 | 0
0 0 1 | 3
|
得出$x_3=3$，然后将其回代入第二个式子得到$x_2$，再回代入第一个式子得到$x_1$
代码
P3389 【模板】高斯消元法
code
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
| #include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=105;
const double eps=1e-6;
int n;
double a[maxn][maxn],ans[maxn];
inline void Gauss(){
for(int i=1;i<=n;i++){
int tp=i;
for(int j=i+1;j<=n;j++){
if(fabs(a[j][i])>fabs(a[tp][i]))
tp=j;
}
if(fabs(a[tp][i])<eps){
printf("No Solution\n");
return;
}
if(tp!=i) swap(a[tp],a[i]);
double div=a[i][i];
for(int j=i;j<=n+1;j++)
a[i][j]/=div;
for(int j=i+1;j<=n;j++){
div=a[j][i];
for(int k=i;k<=n+1;k++){
a[j][k]-=a[i][k]*div;
}
}
}
ans[n]=a[n][n+1];
for(int i=n-1;i>=1;i--){
ans[i]=a[i][n+1];
for(int j=i+1;j<=n;j++){
ans[i]-=ans[j]*a[i][j];
}
}
for(int i=1;i<=n;i++)
printf("%.2lf\n",ans[i]);
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)
for(int j=1;j<=n+1;j++)
scanf("%lf",&a[i][j]);
Gauss();
return 0;
}
|
高斯-约旦消元
操作步骤和高斯消元相同，只不过最后将增广阵消为对角阵，从而去除了繁复的回代操作，经过证明，该方法精度最高，代码也相对较短
时间复杂度：$O(n^3)$
代码
code
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
| inline void Gauss_Jordan(){
for(int i=1;i<=n;i++){
int tp=i;
for(int j=i+1;j<=n;j++){
if(fabs(a[j][i])>fabs(a[tp][i]))
tp=j;
}
if(a[tp][i]<eps){
print("No Solution\n");
return;
}
if(tp!=i) swap(a[tp],a[i]);
for(int j=1;j<=n;j++){
if(j!=i){
double div=a[j][i]/a[i][i];
for(int k=i;k<=n+1;k++){
a[j][k]-=a[i][k]*div;
}
}
}
}
for(int i=1;i<=n;i++){
printf("%.2lf\n",a[i][n+1]/a[i][i]);
}
}
|
矩阵求逆
学过线代的我们都知道，一个矩阵的逆矩阵的求法（前提是该矩阵可逆，可以通过行列式值不为零或满秩等判断），最常用的就是在它右边拼一个同大小单位阵，然后通过初等行变换把左侧矩阵化为单位阵，则此时右侧矩阵即原矩阵的逆矩阵（证明略）
对于这样一个情况，我们很容易发现这就是一个高斯-约旦消元的过程，只不过右侧不是一个列向量，而是一个同阶矩阵，那么如法炮制，很容易得到解法
代码
P4783 【模板】矩阵求逆
这个题需要注意的是有个模数，因此不需要用double存矩阵元素，除法用扩欧或费马小定理算下逆元就行了，注意出现负数，我在一开始就把他全部转正了
code
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
| #include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
int n;
ll a[405][805];
ll exgcd(ll &x,ll &y,ll a,ll b){
if(b==0){
x=1,y=0;return a;
}
ll res=exgcd(x,y,b,a%b);
ll tx=x;
x=y;
y=tx-a/b*y;
return res;
}
inline void Gauss_Jordan(){
for(int i=1;i<=n;i++){
int tp=i;
for(int j=i+1;j<=n;j++){
if(a[j][i]>a[tp][i])
tp=j;
}
if(a[tp][i]==0){
printf("No Solution\n");
return;
}
if(tp!=i) swap(a[tp],a[i]);
for(int j=1;j<=n;j++){
if(j!=i){
ll x,y;
exgcd(x,y,a[i][i],mod);
x=(x+mod)%mod;
ll div=a[j][i]*x%mod;
for(int k=i;k<=n+n;k++){
a[j][k]-=a[i][k]*div%mod;
a[j][k]=(a[j][k]+mod)%mod;
}
}
}
}
for(int i=1;i<=n;i++){
ll x,y;
exgcd(x,y,a[i][i],mod);
x=(x+mod)%mod;
for(int j=1;j<=n;j++){
printf("%lld ",a[i][j+n]*x%mod);
}
printf("\n");
}
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
scanf("%lld",&a[i][j]),
a[i][j]=(a[i][j]+mod)%mod;
for(int i=1;i<=n;i++) a[i][i+n]=1;
Gauss_Jordan();
return 0;
}
|
线性基
参考blog：
这是一个比较偏程序设计竞赛的算法（应该叫数据结构），可以用来存储一组集合的一种基，类似于线代里面的线性无关组，线性基内的元素通过异或操作能不多不少恰好表示出整个集合的所有元素
构造方式
通过不断插入集合元素来完成构造
从集合内元素最高位开始，往低位数，如果插入的数字$x$当前位为1，则判断是否该位的线性基已存在，若存在，则让$x$异或该基，看下一位，重复上述操作；否则令该位线性基为$x$，结束
| inline void insert(ll x){
for(int i=60;i>=0;i--){
if((1ll<<i)&x){
if(p[i]) x^=p[i];
else{
p[i]=x;
break;
}
}
}
}
|
经过构造过程，我们可以得到一组集合的线性基，该线性基集合中存储了从高位到低位每一位为1的所有信息，并且这些信息都是通过原集合中若干元素异或得到的，因此，可以通过线性基表示出原集合的所有元素（具体证明见上方blog）
给定集合取出任意个求最大异或和
P3812 【模板】线性基
构造好线性基后，只需要通过从高位到低位的贪心过程就可以AC此题，具体贪心流程为：
初始解为ans=0，从最高位开始遍历，如果(p[i]^ans)>ans则让ans^=p[i]，最后得到的ans就是答案
证明：显然可以发现，不断使得当前最高位变大的过程是最优的，因此只需要不断取能使答案变大的就行了
PS：^运算符优先级低于>，因此要加括号
code
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
| #include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int n;
ll p[65];
ll x;
inline void ins(ll x){
for(int i=50;i>=0;i--){
if((1ll<<i)&x){
if(p[i]) x^=p[i];
else{
p[i]=x;
break;
}
}
}
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)
scanf("%lld",&x),ins(x);
ll ans=0;
for(int i=50;i>=0;i--)
if((ans^p[i])>ans) ans^=p[i];
printf("%lld\n",ans);
return 0;
}
|