【笔记】矩阵快速幂

矩阵快速幂

矩阵乘法 + 快速幂

矩阵加法：

定义矩阵 \(C=A+B\)。

\(C_{i,j}=A_{i,j}+B_{i,j}\)。

矩阵乘法：

计算两个矩阵的乘法。\(n \times m\) 阶的矩阵 \(A\) 乘以 \(m \times k\) 阶的矩阵 \(B\) 得到的矩阵 \(C\) 是 \(n \times k\) 阶的，且 \(C[i][j]=A[i][0] \times B[0][j]+A[i][1] \times B[1][j]+\) …… \(+A[i][m-1] \times B[m-1][j](C[i][j]\) 表示 \(C\) 矩阵中第 \(i\) 行第 \(j\) 列元素）。

（摘自 Luogu B2105 题面）。

即：第 \(i\) 行与第 \(j\) 列对应乘起来加和放在 \((i,j)\)。

第一个矩阵为 \(n \times m\)，第二个矩阵式 \(m \times k\)，则两个矩阵相乘为：

for(int i=1;i<=n;i++){
	for(int j=1;j<=k;j++){
		for(int d=1;d<=m;d++){
			c[i][j]+=a[i][d]*b[d][j];
		}
	}
}

时间复杂度 \(O(n^3)\)。

Luogu P3390【模板】矩阵快速幂

把快速幂的乘法直接换成矩阵乘法即可。

但是要注意，结果的初始值要定义为对角线全是 \(1\) 的矩阵（单位矩阵）。单位矩阵乘任何矩阵都等于那个矩阵。

#include<iostream>
#include<cstring>
#define int long long 
using namespace std;
const int p=1e9+7;
int n,k;
int a[110][110],s[110][110],c[110][110];
void sta(){
	memset(c,0,sizeof(c));
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			for(int d=1;d<=n;d++){
				c[i][j]=(c[i][j]+(s[i][d]*a[d][j])%p)%p;
			}
		}
	}
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++) s[i][j]=c[i][j];
	}
}
void ata(){
	memset(c,0,sizeof(c));
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			for(int d=1;d<=n;d++){
				c[i][j]=(c[i][j]+(a[i][d]*a[d][j])%p)%p;
			}
		}
	}
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++) a[i][j]=c[i][j];
	}
}
signed main(){
	cin>>n>>k;
	for(int i=1;i<=n;i++) s[i][i]=1;
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			cin>>a[i][j];
		}
	}
	while(k){
		if(k&1) sta();
		ata();
		k>>=1;
	}
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++) cout<<s[i][j]<<" ";
		cout<<endl;
	}
	return 0;
}

运算符重载：把乘号重载成矩阵乘法的形式。

方法：把矩阵都定义成结构体，然后写运算符重载函数。

struct Matrix{
	int mx[110][110];
}a,s,c;

Matrix operator *(const Matrix &a,const Matrix &b){
	memset(c.mx,0,sizeof(c.mx));
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			for(int d=1;d<=n;d++){
				c.mx[i][j]=(c.mx[i][j]+(a.mx[i][d]*b.mx[d][j])%p)%p;
			}
		}
	}
	return c;
}

这样快速幂就能写成：

while(k){
	if(k&1) s=s*a;
	a=a*a;
	k>>=1;
}

同时删去了两个重复的函数及来回赋值的过程，优化了时间常数。

矩阵加速递推

可以用矩阵快速幂来加速递推。

构造矩阵的通用技巧：

先把递推式中出现的要推下去的元素列成一排，构造状态矩阵；
根据递推关系，把状态矩阵转移一次后的矩阵（一次递推）用状态矩阵中的元素表示；
状态矩阵中第 \(i\) 个元素转移一次的式子，是由状态矩阵中的元素构成的多项式，转移第 \(i\) 项的多项式的第 \(j\) 项对应的系数填在转移矩阵的第 \(i\) 列第 \(j\) 行。

Luogu P1962 斐波那契数列

递推式要不断使用 \(f_{n-1}\) 和 \(f_{n-2}\) 推出 \(f_n\)，因此状态矩阵中要维护这两项。转移一次后，\(f_{n-1}\) 和 \(f_{n-2}\) 变成 \(f_n\) 和 \(f_{n-1}\)。

根据算式：

\[f_n=1\times f_{n-1}+1\times f_{n-2} \]

\[f_{n-1}=1\times f_{n-1}+0\times f_{n-2} \]

可以求出每一项的系数，按列填到矩阵里：

\[\begin{bmatrix} f_n & f_{n-1} \end{bmatrix} = \begin{bmatrix} f_{n-1} & f_{n-2} \end{bmatrix} \times \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \]

#include<iostream>
#include<cstring>
#define int long long 
using namespace std;
const int p=1e9+7;
struct Matrix{
	int mx[3][3];
}a,s,c;
int n;
Matrix operator *(const Matrix &a,const Matrix &b){
	memset(c.mx,0,sizeof(c.mx));
	for(int i=1;i<=2;i++){
		for(int j=1;j<=2;j++){
			for(int d=1;d<=2;d++){
				c.mx[i][j]=(c.mx[i][j]+(a.mx[i][d]*b.mx[d][j])%p)%p;
			}
		}
	}
	return c;
}
signed main(){
	cin>>n;
	n-=1;
	s.mx[1][1]=s.mx[1][2]=s.mx[2][1]=1;s.mx[2][2]=0;
	a.mx[1][1]=a.mx[1][2]=1;
	while(n){
		if(n&1) a=s*a;
		s=s*s;
		n>>=1;
	}
	cout<<a.mx[1][1];
	return 0;
}

Luogu P5175 数列

定义 \(f(n)\) 表示 \(a_n\)，\(s(n)\) 表示 \(\sum_{i=1}^na_i^2\)。

根据完全平方公式写出递推式：

\[f(n)=xf(n-1)+yf(n-2) \]

\[f(n)^2=(xf(n-1)+yf(n-2))=x^2f(n-1)^2+y^2f(n-2)^2+2xyf(n-1)f(n-2) \]

\[f(n)f(n-1)=xf(n-1)^2+yf(n-1)f(n-2) \]

\[s(n)=s(n-1)+f(n)^2 \]

由此构造状态矩阵。发现求 \(s(n)\) 需要用到四个值：\(s(n-1)\)、\(f(n-1)^2\)、\(f(n-2)^2\)、\(f(n-1)f(n-2)\)。

即：

\[\begin{bmatrix} s(n) & f(n)^2 & f(n-1)^2 & f(n)f(n-1) \end{bmatrix} = \begin{bmatrix} s(n-1) & f(n-1)^2 & f(n-2)^2 & f(n-1)f(n-2) \end{bmatrix} \times A \]

根据递推式填上转移矩阵 \(A\)。

\[A=\begin{bmatrix} 1 & 0 & 0 & 0 \\ x^2 & x^2 & 1 & x \\ y^2 & y^2 & 0 & 0 \\ 2xy & 2xy & 0 & y \end{bmatrix} \]

#include<iostream>
#include<cstdio>
#include<cstring>
#define int long long 
using namespace std;
const int p=1e9+7;
struct Matrix{
	int mx[5][5];
}a,s,c;
int T,n,k,f1,f2,x,y;
Matrix operator *(const Matrix &a,const Matrix &b){
	memset(c.mx,0,sizeof(c.mx));
	for(int i=1;i<=4;i++){
		for(int j=1;j<=4;j++){
			for(int d=1;d<=4;d++){
				c.mx[i][j]=(c.mx[i][j]%p+(a.mx[i][d]*b.mx[d][j])%p)%p;
			}
		}
	}
	return c;
}
signed main(){
	cin>>T;
	while(T--){
		memset(s.mx,0,sizeof(s.mx));
		memset(a.mx,0,sizeof(a.mx));
		scanf("%lld%lld%lld%lld%lld",&n,&f1,&f2,&x,&y);
		f1%=p;f2%=p;x%=p;y%=p;
		for(int i=1;i<=4;i++) s.mx[i][i]=1;
		a.mx[1][1]=a.mx[2][3]=1;
		a.mx[2][1]=a.mx[2][2]=x*x%p;
		a.mx[3][1]=a.mx[3][2]=y*y%p;
		a.mx[4][1]=a.mx[4][2]=2*x*y%p;
		a.mx[2][4]=x;a.mx[4][4]=y%p;
		if(n==1) printf("%lld\n",(f1*f1)%p);
		else if(n==2) printf("%lld\n",(f1*f1%p+f2*f2%p)%p);
		else{
			int s2=(f1*f1%p+f2*f2%p)%p;
			int f22=(f2*f2)%p;
			int f11=(f1*f1)%p;
			int f12=(f1*f2)%p;
			n-=2;
			while(n){
				if(n&1) s=s*a;
				a=a*a;
				n>>=1;
			}
			printf("%lld\n",(((s2*s.mx[1][1])%p+(f22*s.mx[2][1]))%p+((f11*s.mx[3][1])%p+(f12*s.mx[4][1])%p)%p)%p);
		}
	}
	return 0;
}

Luogu P2233 [HNOI2002] 公交车路线

构造每个站台之间的邻接矩阵，把矩阵平方，就能找到哪些站台到其他站台能换 2 次。把这个矩阵乘 \(n\) 次然后输出 \([1][5]\)。

注意构造矩阵时，E 站台不能有出边。

运用：乘法原理/加法原理。

#include<iostream>
#include<cstring> 
using namespace std;
const int p=1000;
struct Matrix{
	int mx[10][10];
}c,s,a;
Matrix operator *(const Matrix &a,const Matrix &b){
	memset(c.mx,0,sizeof(c.mx));
	for(int i=1;i<=8;i++){
		for(int j=1;j<=8;j++){
			for(int d=1;d<=8;d++){
				c.mx[i][j]=(c.mx[i][j]%p+(a.mx[i][d]*b.mx[d][j])%p)%p;
			}
		}
	}
	return c;
}
int n;
int main(){
	cin>>n;
	for(int i=1;i<=8;i++) s.mx[i][i]=1;
	a.mx[1][2]=a.mx[2][1]=a.mx[2][3]=a.mx[3][2]=a.mx[3][4]=a.mx[4][3]=a.mx[4][5]=a.mx[6][5]=a.mx[6][7]=a.mx[7][6]=a.mx[7][8]=a.mx[8][7]=a.mx[8][1]=a.mx[1][8]=1;
	while(n){
		if(n&1) s=s*a;
		a=a*a;
		n>>=1;
	}
	printf("%d",s.mx[1][5]%p);
	return 0;
}

Luogu P2886 [USACO07NOV] Cow Relays G

Floyd 算法最外层循环相当于加点，每加一个点就在原邻结矩阵上更新所有两点间的最短路。

只多走一条路，也就是每两点间只加一个点，不能在原矩阵上更新，而是要存在初始全为互相不可达的新矩阵中（确保新加点一定优，一定被新加点更新而不是保留边数更少距离更短的值），确保新矩阵里都是正好经过上一次 \(+1\) 条边的答案。

于是形式上就变成矩阵乘法了。每多走一条路，相当于做一次求 \(\min\) 的当前矩阵 \(\times\) 原矩阵的矩阵乘法。

这样的矩阵乘法也可以快速幂优化，因为 \(\min\) 也满足结合律。

最后需要注意，每个点一开始是不能走到自身的。

#include<bits/stdc++.h>
using namespace std;
const int N=1e3+10;
const int INF=1e9;
int n,m,s,t,k;
int tot;
int nd[N],idx;
struct Node{
    int from,to,w;
}e[510];
struct Matrix{
    int mx[110][110];
}a,b,c;
Matrix operator *(const Matrix &a,const Matrix &b){
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++) c.mx[i][j]=INF;
    }
    //for(int i=1;i<=n;i++) c.mx[i][i]=0;
    for(int d=1;d<=n;d++){
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++){
                c.mx[i][j]=min(c.mx[i][j],a.mx[i][d]+b.mx[d][j]);
            }
        }
    }
    return c;
}
void Add(int u,int v,int w){
    tot++;
    e[tot].from=u;
    e[tot].to=v;
    e[tot].w=w;
}
void fstpow(){
    int q=k-1;
    b=a;
    while(q){
        if(q&1) b=b*a;
        a=a*a;
        q>>=1;
    }
}
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0),cout.tie(0);
    cin>>k>>m>>s>>t;
    for(int i=1;i<=m;i++){
        int u,v,w;
        cin>>w>>u>>v;
        Add(u,v,w);
        nd[++idx]=u,nd[++idx]=v;
    }
    sort(nd+1,nd+1+idx);
    n=unique(nd+1,nd+1+idx)-nd-1;
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++) a.mx[i][j]=b.mx[i][j]=INF;
    }
    //for(int i=1;i<=n;i++) a.mx[i][i]=0,b.mx[i][i]=0;
    for(int i=1;i<=m;i++){
        int u=lower_bound(nd+1,nd+1+n,e[i].from)-nd;
        int v=lower_bound(nd+1,nd+1+n,e[i].to)-nd;
        a.mx[u][v]=a.mx[v][u]=e[i].w;
    }
    s=lower_bound(nd+1,nd+1+n,s)-nd;
    t=lower_bound(nd+1,nd+1+n,t)-nd;
    fstpow();
    cout<<b.mx[s][t];
    return 0;
}

posted @ 2025-12-12 23:25 Seqfrel 阅读(38) 评论(0) 收藏举报

刷新页面返回顶部

Loading

Seqfrel