矩阵学习笔记
矩阵学习笔记
持续更新
3Brown1Blue线性代数系列值得你拥有
高斯消元
高斯消元算法的思想是，对于每一个未知量 \(x_i\) ,找到一个 \(x_i\) 的系数非零，但 \(x_1-x_{i-1}\) 的系数为0的方程，然后用初等行变换把其他方程的 \(x_i\) 的系数全部消去
例题
- 球形空间产生器
将 \(n+1\) 个方程彼此相减得到 \(n\) 个线性方程，然后高斯消元板子即可
#include<bits/stdc++.h>
using namespace std;
double eps=1e-10;
double a[20][20],b[20][20],c[20];
int n;
//a是原数组，b为系数矩阵，c为常数
int main(){
scanf("%d",&n);
for(int i=1;i<=n+1;i++)
for(int j=1;j<=n;j++)scanf("%lf",&a[i][j]);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
b[i][j]=2.0*(a[i][j]-a[i+1][j]);
c[i]+=(a[i][j]*a[i][j]-a[i+1][j]*a[i+1][j]);
}
}
for(int i=1;i<=n;i++){
for(int j=i;j<=n;j++){
if(b[j][i]>eps){
for(int k=1;k<=n;k++)swap(b[i][k],b[j][k]);
swap(c[i],c[j]);break;
}
}
for(int j=1;j<=n;j++){
if(i==j)continue;
double bei=b[j][i]/b[i][i];
for(int k=1;k<=n;k++)b[j][k]-=b[i][k]*bei;
c[j]-=c[i]*bei;
}
}
for(int i=1;i<=n;i++){
printf("%.3lf ",c[i]/b[i][i]);
}
printf("\n");
return 0;
}
矩阵加速dp
这个还是直接看例题吧
例题
1. hh去散步
题意：在一张无向图中，有重边，要求从起点出发，且同一条边不能连续走两次，问经过T次边到达终点的路径一共有几条。 \(n<=50,m<=60,T<=2^{30}\)
k短路！！！ 看看标题，让我们用dp来想一想。
假设没有 “有重边、同一条边不能连续走两次” 的条件，设 \(f_{i,j}\) 表示已走 \(i\) 条边，现在在 \(j\) 节点的方案数，\(to(j)\) 表示能够经过一条边到达 \(j\) 的点的集合，则有式子 \(f_{i+1,j}=\sum_{k\in to(j)}f_{i,k}\)
那如果有呢？因为有重边的限制，所以我们可以采用一种很妙的思想，叫做 点边转换 ，将上述过程中的节点变为边。
即设 \(f_{i,j}\) 表示已走 \(i\) 条边，现在走完了第 \(j\) 号边的方案数，\(to(j)\) 表示能够直接到达 \(j\) 号边的集合，则有式子 \(f_{i+1,j}=\sum_{k\in to(j)}f_{i,k}\)
只是小小的改变，我们便解决了重边的问题。可是一看数据范围 \(T<=2^{30}\) ！！！怕得要换银河才能过
再看看 \(n,m\) 的数据范围，这就是经典的矩阵加速啊
不过我们也得想想为什么可以用矩阵。 对于所有 i 确定的 \(f_{i,j}\)，我们可以用一个 2*m 的向量来表示当前状态。而转移方程呢？想想是不是可以用一个矩阵来表示其转移关系呢，然后直接相乘呢？不过在这里我们得要弄清楚矩阵每一位的意义，我们才能写对这一个转移矩阵。
设这个 \(2m*2m\) 的转移矩阵为 \(A\)，考虑 \(f*A\),结果为向量 \(f'\) ,则 \(\forall i\in[1,n],f'_i=\sum_{j=1}^{2m}f_j*A_{j,i}\) 。那么\(A_{j,i}\) 的意义就为已经走到第 j 号边，现在要走到第 i 号边的方案数。再看条件“同一条边不能连续走两次”，我们得要特判一下反向边
这样我们就可以写出向量 \(f\) 的转移矩阵。 \(f_T=f_1*A^{T-1}\)。又因为矩阵具有结合律，所以我们可以用矩阵快速幂算出 \(A^T\)，再与 \(f\) 相乘，最后再统计答案。时间复杂度为 \(O(m^3logT)\)
注意:
-
使用memcpy()一定要确保操作的两个数组空间大小相同的情况下使用，否则会有想不到的错误（修改越界）发生
-
一定要理解矩阵代表的具体含义;矩阵乘法不满足交换律
-
异或(^)的优先级很低，注意要打括号
-
因为定义，我们需要预处理出初始状态
-
感觉写结构体封装矩阵或许更好
代码：
/*
* @Author: fpjo
* @Date: 2020-08-24 19:23:41
* @Last Modified by: fpjo
* @Last Modified time: 2020-08-24 21:09:15
*/
#include<bits/stdc++.h>
using namespace std;
#define int long long
int const MOD=45989;
int n,m,t,a,b,tot=1,ans=0;
int h[200],mat[200][200],ori[200];
struct edge{
int next,to,be;
}e[1000];
void add(int u,int v){
e[++tot].to=v,e[tot].next=h[u],h[u]=tot;
}
void mul(){
int x[200];
memset(x,0,sizeof(x));
for(int j=2;j<=tot;j++){
for(int k=2;k<=tot;k++){
x[j]=(x[j]+(ori[k]*mat[k][j])%MOD)%MOD;
}
}
for(int i=2;i<=tot;i++){
ori[i]=x[i];
}
}
void mulself(){
int x[200][200];memset(x,0,sizeof(x));
for(int i=2;i<=tot;i++){
for(int j=2;j<=tot;j++){
for(int k=2;k<=tot;k++){
x[i][j]=(x[i][j]+(mat[i][k]*mat[k][j])%MOD)%MOD;
}
}
}
for(int i=2;i<=tot;i++){
for(int j=2;j<=tot;j++){
mat[i][j]=x[i][j];
}
}
}
signed main(){
scanf("%lld%lld%lld%lld%lld",&n,&m,&t,&a,&b);
for(int i=1;i<=m;i++){
int u,v;scanf("%lld%lld",&u,&v);
add(u,v);add(v,u);
}
for(int i=2;i<=tot;i++){
int to=e[i].to;
for(int j=h[to];j;j=e[j].next){
if(j!=(i^1)){
mat[i][j]=1;
}
}
}
for(int i=h[a];i;i=e[i].next)ori[i]=1;
t-=1;
for(;t;t>>=1){
if(t&1)mul();
mulself();
}
for(int i=h[b];i;i=e[i].next){
ans=(ans+ori[i^1])%MOD;
}
printf("%lld\n",ans);
return 0;
}
矩阵与数据结构相结合
这个就直接看例题了吧
例题
1. 大魔法师
维护区间上的三亚组 \(A,B,C\) ,满足以下区间操作:
- \(A_i=A_i+B_i\)
- \(B_i=B_i+C_i\)
- \(C_i=C_i+A_i\)
- \(A_i=A_i+v(v为常数)\)
- \(B_i=B_i*v(同上)\)
- \(C_i=v(同上)\)
- 查询 \([l,r]\)
码量对于现在的我来说有点大,100+行.
实际上思路很简单,用线段树维护区间,每个节点用一个向量 \([A,B,C,len]\),然后用矩阵分别表示前6个操作即可
就是细节有亿点点多,写了我一个晚上加一个上午
注意取模可以优化,但我没写,所以贼慢
下面贴代码
/*
* @Author: fpjo
* @Date: 2020-08-27 08:13:18
* @Last Modified by: fpjo
* @Last Modified time: 2020-08-27 09:55:43
*/
#include<bits/stdc++.h>
using namespace std;
int const N=25*1e4+10,MOD=998244353;
int n,m,ansa,ansb,ansc;
int a[N][3];
// tag实际上是一个矩阵!而不是opt
struct Mat{
int mat[4][4];
int l;
Mat() {memset(mat,0,sizeof(mat));l=4;}
const Mat operator *(const Mat &a){
Mat x;
for(int i=0;i<l;i++)
for(int j=0;j<l;j++)
for(int k=0;k<l;k++)
if(mat[i][k] && a.mat[k][j])
x.mat[i][j]=(x.mat[i][j]+(long long)mat[i][k]*a.mat[k][j]%MOD)%MOD;
return x;
}
void operator *= (const Mat &a){
*this=(*this) * a;
}
void operator =(const Mat &a){
for(int i=0;i<l;i++)
for(int j=0;j<l;j++)
mat[i][j]=a.mat[i][j];
}
//调试
void print(){
for(int i=0;i<l;i++){
for(int j=0;j<l;j++){
printf("%d ",mat[i][j]);
}
printf("\n");
}
printf("--------------------\n");
}
}trans,ori;
struct seg{
int l,r;
Mat m,tag;
}seg[N*10];
inline int read(){
int x=0,f=1;char c=getchar();
while(c<'0' || c>'9'){if(c=='-')f=-1;c=getchar();}
while(c>='0' && c<='9'){x=(x<<1)+(x<<3)+c-'0';c=getchar();}
return x*f;
}
void build(int x,int l,int r){
seg[x].tag=ori;
seg[x].l=l,seg[x].r=r;
if(l==r){
seg[x].m.mat[0][0]=a[l][0],seg[x].m.mat[0][1]=a[l][1],seg[x].m.mat[0][2]=a[l][2],seg[x].m.mat[0][3]=1;
return;
}
int mid=(l+r)>>1;
build(x<<1,l,mid);build(x<<1|1,mid+1,r);
for(int i=0;i<4;i++){
seg[x].m.mat[0][i]=(seg[x<<1].m.mat[0][i]+seg[x<<1|1].m.mat[0][i])%MOD;
}
return;
}
void change(int t,int v){
if(t==1)trans.mat[1][0]=1;
else if(t==2)trans.mat[2][1]=1;
else if(t==3)trans.mat[0][2]=1;
else if(t==4)trans.mat[3][0]=v;
else if(t==5)trans.mat[1][1]=v;
else if(t==6)trans.mat[3][2]=v,trans.mat[2][2]=0;
}
void init(){
trans=ori;
}
void pushdown(int x){
seg[x<<1].tag*=seg[x].tag,seg[x<<1|1].tag*=seg[x].tag;
seg[x<<1].m*=seg[x].tag,seg[x<<1|1].m*=seg[x].tag;
seg[x].tag=ori;
}
void solve(int x,int l,int r){
if(l<=seg[x].l&&seg[x].r<=r){
seg[x].tag*=trans;
seg[x].m*=trans;
return;
}
pushdown(x);
int mid=(seg[x].l+seg[x].r)>>1;
if(l<=mid)solve(x<<1,l,r);
if(r>mid)solve(x<<1|1,l,r);
for(int i=0;i<4;i++)seg[x].m.mat[0][i]=((long long)seg[x<<1].m.mat[0][i]+seg[x<<1|1].m.mat[0][i])%MOD;
}
void query(int x,int l,int r){
if(l<=seg[x].l && seg[x].r<=r){
ansa=(ansa+seg[x].m.mat[0][0])%MOD,ansb=(ansb+seg[x].m.mat[0][1])%MOD,ansc=(ansc+seg[x].m.mat[0][2])%MOD;
return;
}
pushdown(x);
int mid=(seg[x].l+seg[x].r)>>1;
if(l<=mid)query(x<<1,l,r);
if(r>mid)query(x<<1|1,l,r);
return;
}
int main(){
ori.mat[0][0]=1,ori.mat[1][1]=1,ori.mat[2][2]=1,ori.mat[3][3]=1;
init();
n=read();
for(int i=1;i<=n;i++)a[i][0]=read(),a[i][1]=read(),a[i][2]=read();
build(1,1,n);
scanf("%d",&m);
for(int i=1;i<=m;i++){
int opt,l,r,v=-1;opt=read(),l=read(),r=read();
if(opt>=4 && opt<=6)v=read();
if(opt<7){
change(opt,v);
solve(1,l,r);
init();
}else{
ansa=0,ansb=0,ansc=0;
query(1,l,r);
printf("%d %d %d\n",ansa,ansb,ansc);
}
}
return 0;
}
尾声
学长说矩阵运用的话基本就这三个方面，所以算完结撒花啦！
当然练习题到时候会补的
