ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

【笔记】线性代数

2022-08-27 13:02:19  阅读:141  来源: 互联网

标签:node frac int ll 矩阵 笔记 线性代数 ans


矩阵乘法

首先给出矩阵乘法的代数意义:

结合一个具体的例子来理解:

设答案矩阵为 \(ans\) 。根据公式:

\(ans_{1,1}\) 是由 \(a\) 矩阵的第一行与 \(b\) 矩阵的第一列逐位相乘并求和得到的。

\(ans_{1,2}\) 是由 \(a\) 矩阵的第一行与 \(b\) 矩阵的第二列逐位相乘并求和得到的。

\(ans_{2,1}\) 是由 \(a\) 矩阵的第二行与 \(b\) 矩阵的第一列逐位相乘并求和得到的。

\(ans_{2,2}\) 是由 \(a\) 矩阵的第二行与 \(b\) 矩阵的第二列逐位相乘并求和得到的。

由此,可以简单理解,矩阵乘法得出的结果,\(ans_{i,j}\) 是由 \(a\) 矩阵的第 \(i\) 行与 \(b\) 矩阵的第 \(j\) 列逐位相乘并求和得到的。同时,也可以得知,一般情况下,要做矩阵乘法,\(a\) 矩阵的列数应该和 \(b\) 矩阵的行数相同。

根据这一结论,我们就可以写出矩阵乘法的代码:

struct node{
	ll p[105][105];
};
node X(node a,node b){
	node t;
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			t.p[i][j]=0;
			for(int k=1;k<=n;k++){
				t.p[i][j]+=a.p[i][k]*b.p[k][j];
				t.p[i][j]%=mod;
			}
		}
	}
	return t;
}

矩阵快速幂

了解了矩阵乘法,就可以做矩阵快速幂了。

对于矩阵 \(A\) ,\(A^k\) 就表示 \(k\) 个 \(A\) 矩阵相乘。

先回顾一下快速幂。快速幂的做法是将指数每次折半。这样可以把效率提到 \(log n\)。

例如求 \(a^{2n+1}\) 时,就可以转成 \(a^{n}\times a^{n}\times a\)。求\(a^{2n}\) 时,就可以转成 \(a^{n}\times a^{n}\)。

根据这个思路,就有了快速幂板子:

P1226 【模板】快速幂||取余运算

ll fpow(ll a,ll p){
	ll ans=1;
	while(p){
		if(p&1) ans=ans*a%mod;
		a=a*a%mod;
		p>>=1;
	}	
	return ans;
}

至于光速幂,咱也不懂qwq
光速幂

回归正题,那么矩阵快速幂也同理。只要把普通快速幂中的乘法换成矩阵乘法就可以了。写成代码就是这样:

node fpow(node a,ll k){
	node ans=a,b=a;
	while(k){
		if(k&1) ans=X(b,ans);
		b=X(b,b);
		k>>=1;
	}	
	return ans;
}

需要注意的是,\(k\) 的值在传入时,需要 \(-1\) ,因为矩阵的一次幂就是本身,相当于已经乘过一次了。

结合矩阵乘法,我们就可以得到矩阵快速幂的完整板子:

P3390 【模板】矩阵快速幂

#include<bits/stdc++.h>
#define ll long long
#define mod 1000000007
using namespace std;
ll n,k;
struct node{
	ll p[105][105];
}a;
node X(node a,node b){
	node t;
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			t.p[i][j]=0;
			for(int k=1;k<=n;k++){
				t.p[i][j]+=a.p[i][k]*b.p[k][j];
				t.p[i][j]%=mod;
			}
		}
	}
	return t;
}
node fpow(node a,ll k){
	node ans=a,b=a;
	while(k){
		if(k&1) ans=X(b,ans);
		b=X(b,b);
		k>>=1;
	}	
	return ans;
}
int main(){
    scanf("%lld%lld",&n,&k);
    for(int i=1;i<=n;i++){
    	for(int j=1;j<=n;j++) scanf("%lld",&a.p[i][j]);
	}
	a=fpow(a,k-1);
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++) printf("%lld ",a.p[i][j]);
		printf("\n");
	}
    return 0;
}

一道用到矩阵快速幂的例题:

P1962 斐波那契数列

十分巧妙地用了矩阵快速幂加速递推过程。具体可以看题解,这里不多赘述。

#include<bits/stdc++.h>
#define ll long long
#define mod 1000000007
using namespace std;
ll n;
struct node{
	ll p[15][15];
}a,ans;
node X(node a,node b){
	node t;
	for(int i=1;i<=2;i++){
		for(int j=1;j<=2;j++){
			t.p[i][j]=0;
			for(int k=1;k<=2;k++){
				t.p[i][j]+=a.p[i][k]*b.p[k][j];
				t.p[i][j]%=mod;
			}
		}
	}
	return t;
}
void fpow(ll k){
	while(k){
		if(k&1) ans=X(a,ans);
		a=X(a,a);
		k>>=1;
	}
}
int main(){
    scanf("%lld",&n);
    if(n<=2){
    	printf("1");
    	return 0;
	}
    a.p[1][1]=a.p[1][2]=a.p[2][1]=1;
    ans.p[1][1]=ans.p[1][2]=1;
    fpow(n-1);
    printf("%lld",ans.p[1][1]);
    return 0;
}

高斯消元

高斯消元法用于解线性方程组。

那么什么是线性方程组呢?

线性方程组就是有多个未知数,并且每个未知数的次数均为一次,这样多个未知数组成的方程组为线性方程组。或者我们也可以叫它多元一次方程组。

比如以下这个方程组:

数学老师教我们,多元一次方程组可以用加减消元法和代入消元法求解。

高斯消元法其实就是这样求解的。先进行加减消元,我们可以先求得一个未知数的值,然后可以逐层往回代(代入消元法),依次可以得到第 \(2\) 个、第 \(3\) 个未知数的值,最终解出方程组中各个未知数的值,结束算法。

那如何具体实现?


首先,提出各项系数并转化成一个矩阵。比如上面的方程就可以转化为:

\(\begin{bmatrix} 2&1&-1&|&8\\-3&-1&2&|&11\\-2&1&2&|&-1\end{bmatrix}\)

然后考虑我们做数学时一般解方程思路。我们往往是将系数绝对值最大的方程转移到被减的这一行,方便计算,也可以减小误差。

所以接下来要做的是将 \(x\) 项系数绝对值最大的放到第一行。变化之后矩阵变成这样:

\(\begin{bmatrix} -3&-1&2&|&11\\2&1&-1&|&8\\-2&1&2&|&-1\end{bmatrix}\)

这一步代码实现非常简单,若当前解的是第 \(i\) 个未知数,只要找出 $\max ${ \(fabs(A_{j,i})\) } 即可。其中 \(fabs\) 就是对浮点数取绝对值。设 \(j=r\) 时取到最大,就将第 \(r\) 行与第 \(i\) 行逐项互换就可以了。代码:

r=i;
for(int j=i+1;j<n;j++){
	if(fabs(a[j][i])>fabs(a[r][i])) r=j;
}
if(r!=i) for(int j=0;j<=n;j++) swap(a[i][j],a[r][j]);

之后,对于下面的每一行(记为 \(L_{k}\)),我们可以求出一个比值 \(f=A_{k,i}/A_{i,i}\) 。我们要将 \(L_k\) 的每一项 \(A_{k,j}\) 都减去 $ A_{i,j}\times f $,达到加减消元的目的(可以手算理解)。比如上面这个经处理的矩阵,经过第一次加减消元后会变成:

\(\begin{bmatrix} -3&-1&2&|&11\\0&\frac{1}{3}&\frac{1}{3}&|&\frac{2}{3}\\0&\frac{5}{3}&\frac{2}{3}&|&\frac{13}{3}\end{bmatrix}\)

相当于:

\(L_2\) 变成了 \(L_2-L_1\times(-\frac{2}{3})\)

\(L_3\) 变成了 \(L_3-L_1\times\frac{2}{3}\)


重复上面的步骤:

交换第 \(2,3\) 行(因为\(\frac{5}{3}>\frac{1}{3}\)),得到:

\(\begin{bmatrix} -3&-1&2&|&11\\0&\frac{5}{3}&\frac{2}{3}&|&\frac{13}{3}\\0&\frac{1}{3}&\frac{1}{3}&|&\frac{2}{3}\end{bmatrix}\)

第二次加减消元,得到:

\(\begin{bmatrix} -3&-1&2&|&11\\0&\frac{5}{3}&\frac{2}{3}&|&\frac{13}{3}\\0&0&\frac{1}{5}&|&-\frac{1}{5}\end{bmatrix}\)

相当于:

\(L_3\) 变成了 \(L_3-L_2\times\frac{1}{5}\)


贴上加减消元这整一部分的代码:

for(int i=0;i<n;i++){
	r=i;
	for(int j=i+1;j<n;j++){
		if(fabs(a[j][i])>fabs(a[r][i])) r=j;
	}
	if(r!=i) for(int j=0;j<=n;j++) swap(a[i][j],a[r][j]);
	for(int k=i+1;k<n;k++){
		double f=a[k][i]/a[i][i];
		for(int j=i;j<=n;j++) a[k][j]-=f*a[i][j];
	}
}

下一步就是进行回带了。

先手算:当前矩阵的第三行重新转回数学式子就是 \(\frac{1}{5}z=-\frac{1}{5}\),所以易得 \(z=A_{2,3}/A_{2,2}=-1\)。然后把 \(z=-1\) 带入第二行,可以得到 \(\frac{5}{3}y+\frac{2}{3}\times(-1)=\frac{13}{3}\) 这时我们只要把 \(\frac{13}{3}\) 减掉 \(\frac{2}{3}\times(-1)\),也就是移项,再除以 \(\frac{5}{3}\) 就能得到 \(y\)。以此类推。

所以我们就得到方法:先将等式右边减掉所有等式左边的已知项,再除以未知项系数就好了。具体可以根据代码理解:

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];
	a[i][n]/=a[i][i];
}

最后,我们需要判断无解不唯一解的情况。

仍然从数学上理解,我们知道,对于一次方程,当方程的所有未知数系数都为 \(0\),但常数项不为 \(0\) 时,方程无解。例如 \(0x+0y+0z=114\) ,显然无解。

同样,对于一次方程,当方程的所有未知数系数都为 \(0\),而常数项也为 \(0\) 时,方程有无数解。例如 \(0x+0y+0z=0\) 。

所以,转化到矩阵中,当存在一行 \(L_x\) ,其 \(A_{x,0}\) 到 \(A_{x,n-1}\) 都为 \(0\) 但 \(A_{x,n} \ne 0\) 时,方程组无解。同样,当存在一行 \(L_x\) ,其 \(A_{x,0}\) 到 \(A_{x,n-1}\) 都为 \(0\) 且 \(A_{x,n} = 0\) 时,方程组有无数解。(洛谷 P3389 只问是否有唯一解,所以两种情况可以合并讨论)。

下面是完整代码:

P3389 【模板】高斯消元法

#include<bits/stdc++.h>
using namespace std;
int n,r,fl;
double a[105][105];
void gauss(){
	for(int i=0;i<n;i++){
		r=i;
		for(int j=i+1;j<n;j++){
			if(fabs(a[j][i])>fabs(a[r][i])) r=j;
		}
		if(r!=i) for(int j=0;j<=n;j++) swap(a[i][j],a[r][j]);
		for(int k=i+1;k<n;k++){
			double f=a[k][i]/a[i][i];
			for(int j=i;j<=n;j++) a[k][j]-=f*a[i][j];
		}
	}
	for(int i=0;i<n;i++){
		int kk=0;
		for(int j=0;j<n;j++) if(a[i][j]!=0) kk=1;
		if(!kk){
			fl=1;
			return ;
		}
	}
	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];
		a[i][n]/=a[i][i];
	}
}
int main() {
    scanf("%d",&n);
    for(int i=0;i<n;i++){
    	for(int j=0;j<=n;j++) scanf("%lf",&a[i][j]);
	}
	gauss();
	if(fl){
		printf("No Solution");
		return 0;
	}
	for(int i=0;i<n;i++) printf("%.2lf\n",a[i][n]);
    return 0;
}

线性基

对于一组数 \(A_1\) ~ \(A_n\) ,它们的线性基为 \(P_1\) ~ \(P_n\)。其中 \(P_i\) 表示二进制下出现 \(1\) 的最高位在第 \(i\) 位的数。

所以,线性基可以方便地求最大异或和

先看线性基的构造方法

比如说给定这样一组数构造线性基:\(5,11,17,24\)。

首先将每一个数转化为二进制:

\(5->101\)

\(11->1011\)

\(17->10001\)

\(24->11000\)

然后我们从依次最高位开始找。

如果最高位的 \(1\) 在第 \(i\) 位,并且当前的 \(P_i=0\) 我们就把 \(P_i\) 赋值为这个数。例如处理 \(5\),就把 \(P_2\) 赋值为 \(5\)。

如果当前的 \(P_i \ne 0\) ,那就让当前的数 \(\oplus=P_i\) ,也就是把所有与 \(P_i\) 重复的 \(1\) 抹掉并继续往后搜。
例如处理 \(24\) ,此时 \(P_4\) 已经是 \(17\) 了,就把 \(24 \oplus=17\) 得到 \(9(1001)\)。
又发现此时 \(P_3\) 已经是 \(11\) 了,就把 \(9 \oplus=11\) 得到 \(2(10)\) 。
此时 \(P_1=0\),就把 \(P_1\) 赋值为 \(2\)。

按照手算的过程,我们可以写出代码:

void make(ll a){
	for(int i=52;i>=0;i--){
		if(a>>(ll)i){
			if(!p[i]){
				p[i]=a;
				break;
			}
			else a^=p[i];
		}
	}
	return ;
}

再来看如何寻找答案

简单手算一下,可以发现,求最大异或和是满足贪心的,且对顺序没有要求。

这里可以简单证明:因为异或是不进位运算,所以高位不受低位影响。举几个例子:

  1. \(10001 \oplus 1000 = 11001\) ,显然更优。
  2. \(10111 \oplus 1111 = 11000\) ,虽然低位变小,但高位更大了,所以结果还是更优。
  3. \(11000 \oplus 1111 = 10111\) ,虽然低位变大,但高位更小了,所以结果更劣。

所以,得出:只要将 \(ans\) 对所有 \(P_i\) 从高位到低位异或一遍,就可以得到最优答案。

下面是完整代码:

P3812 【模板】线性基

#include<bits/stdc++.h>
#define ll long long
using namespace std;
ll n,a[60],p[60],ans;
void make(ll a){
	for(int i=52;i>=0;i--){
		if(a>>(ll)i){
			if(!p[i]){
				p[i]=a;
				break;
			}
			else a^=p[i];
		}
	}
	return ;
}
int main() {
    scanf("%d",&n);
    for(int i=1;i<=n;i++){
    	scanf("%lld",&a[i]);
    	make(a[i]);
	}
    for(int i=52;i>=0;i--) ans=max(ans,ans^p[i]);
	printf("%lld",ans);
    return 0;
}

G L & H F

标签:node,frac,int,ll,矩阵,笔记,线性代数,ans
来源: https://www.cnblogs.com/binary1110011/p/16630371.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有