「學習筆記」矩陣乘法與矩陣快速冪
點擊查看目錄
目錄
- 「學習筆記」矩陣乘法與矩陣快速冪
- 矩陣乘
- 演算法
- 代碼
- 矩陣快速冪
- 演算法
- 用處
- 代碼(模板題)
- 練習題
- 斐波那契數列
- 思路
- 代碼
- [SCOI2009] 迷路
- 思路
- 代碼
- 佳佳的 Fibonacci
- 思路
- 代碼
- 選拔隊員(不知道教練從哪里找的)
- 題意
- 思路
- 代碼
- Tr A
- 思路
- 代碼
為什么說文章里多次出現「沖一個矩陣快速冪就行了」:一看斯該喝的 \(Au\) 記就看到了這句話,并對這句話留下了深刻的印象,
釣評論:有人敢幫我@斯該喝嗎,
矩陣乘
演算法
矩陣 \(A\) 規模為 \(n\times m\),矩陣 \(B\) 規模為 \(m\times q\),設 \(C=A\times B\),則:
\[C_{i,j}=\sum_{k=1}^{m}A_{i,k}*B_{k,j} \]代碼
點擊查看代碼
const ll N=110,inf=1ll<<40;
ll n,m,p;
class Matrix{
public:
ll mat[N][N];
inline ll* operator[](ll x){return mat[x];}
inline Matrix operator*(Matrix ma)const{
Matrix mt;
_for(i,1,n){
_for(j,1,p){
mt[i][j]=0;
_for(k,1,m)mt[i][j]+=mat[i][k]*ma[k][j];
}
}
return mt;
}
}a,b;
inline void print(Matrix ma){
_for(i,1,n){
_for(j,1,p)printf("%lld ",ma[i][j]);
puts("");
}
return;
}
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline void In(){
n=rnt(),m=rnt();
_for(i,1,n)_for(j,1,m)a[i][j]=rnt();
p=rnt();
_for(i,1,m)_for(j,1,p)b[i][j]=rnt();
print(a*b);
return;
}
}
矩陣快速冪
演算法
沒啥好說的吧(
多載一下運算子然后沖一個矩陣快速冪就行了,
用處
- 優化遞推式,例:斐波那契數列
代碼(模板題)
點擊查看代碼
const ll N=110,P=1e9+7,inf=1ll<<40;
ll n,k;
class Mat{
public:
ll a[N][N];
inline ll* operator[](ll x){return a[x];}
inline void one(){_for(i,1,n)a[i][i]=1;}
inline Mat operator*(Mat ma){
Mat mt;
_for(i,1,n){
_for(j,1,n){
mt[i][j]=0;
_for(k,1,n)mt[i][j]=(mt[i][j]+a[i][k]*ma[k][j]%P)%P;
}
}
return mt;
}
}a;
inline void printf(Mat ma){
_for(i,1,n){
_for(j,1,n)printf("%lld ",ma[i][j]);
puts("");
}
return;
}
inline Mat fpow(Mat ma,ll b){
Mat ans;ans.one();
while(b){
if(b&1)ans=ans*ma;
ma=ma*ma,b>>=1;
}
return ans;
}
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline void In(){
n=rnt(),k=rnt();
_for(i,1,n)_for(j,1,n)a[i][j]=rnt();
printf(fpow(a,k));
return;
}
}
練習題
斐波那契數列
思路
眾所周知斐波那契數列的遞推式是 \(Fib_n=Fib_{n-1}+Fib_{n-2}\),\(\Theta(n)\) 本可以解決,但本題 \(n<2^{63}\),顯然需要 \(\log_2n\) 演算法,考慮優化,
我們設 \(F(n)\) 表示矩陣 \(\left[Fib_n,Fib_{n-1}\right]\),如果我們要把它變成 \(F(n+1)=\left[Fib_{n+1},Fib_n\right]\),則需要把 \(F(n)_{1,1}\) 挪到 \(F(n+1)_{1,2}\) ,把 \(F(n)_{1,1}+F(n)_{1,2}\) 挪到 \(F(n+1)_{1,1}\),
嘗試用矩陣優化這個東西,
我們可以發現:
\[\begin{aligned} F(n+1)&=\left[\begin{matrix}Fib_n+1&Fib_n\end{matrix}\right]\\ &=\left[\begin{matrix}Fib_n*1+Fib_{n-1}*1&Fib_n*1+Fib_{n-1}*0\end{matrix}\right]\\ &=\left[\begin{matrix}F(n)_{1,1}*1+F(n)_{1,2}*1&F(n)_{1,1}*1+F(n)_{1,2}*0\end{matrix}\right]\\ \end{aligned} \]那么設 \(M=\left[\begin{matrix}1&1\\1&0\end{matrix}\right]\),則:
\[\begin{aligned} F(n+1)&=\left[\begin{matrix}F(n)_{1,1}\times M_{1,1}+F(n)_{1,2}\times M_{2,1}&F(n)_{1,1}\times M_{1,2}+F(n)_{1,2}\times M_{2,2}\end{matrix}\right]\\ &=F(n)\times M \end{aligned} \]然后沖一個矩陣快速冪就行了,
代碼
點擊查看代碼
const ll N=5,inf=1ll<<40;
ll n,m;
class Mat{
public:
ll a[N][N];
inline ll* operator[](ll x){return a[x];}
friend Mat Mul(Mat m1,Mat m2){
Mat an;
_for(i,1,2){
_for(j,1,2){
an[i][j]=0;
_for(k,1,2)an[i][j]=(an[i][j]+m1[i][k]*m2[k][j]%m)%m;
}
}
return an;
}
};
inline Mat fpow(ll b){
Mat ans;ans[1][1]=ans[1][2]=1,ans[2][1]=ans[2][2]=0;
Mat ma;ma[1][1]=ma[1][2]=ma[2][1]=1,ma[2][2]=0;
while(b>0){
if(b&1)ans=Mul(ans,ma);
ma=Mul(ma,ma),b>>=1;
}
return ans;
}
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline void In(){
n=rnt(),m=rnt();
Mat ans=fpow(n-2);
printf("%lld\n",ans[1][1]);
return;
}
}
[SCOI2009] 迷路
思路
設 \(f_{i,j}\) 表示 \(j\) 時刻到點 \(i\) 的方案數,轉移方程:
\[f_{i,j}=\sum_{(i,k)\in\mathbb{E}}f_{k,j-w_{i,k}} \]不是很可做,那我們先考慮邊權只有 \(0\) 和 \(1\)的情況,可以發現轉移方程能直接這樣寫:
\[f_{i,j}=\sum_{1\le j\le n}w_{i,k}\times f_{k,j-1} \]發現又是個矩陣乘法,直接沖一個矩陣快速冪就行了,
但是本題邊權不只有 \(0\) 和 \(1\),不能直接沖,
那么我們對一個點進行拆點,如下圖:

拆成

看起來有點復雜,但懂了之后好理解的,
拆完之后直接沖一個矩陣快速冪就行了,
開做發現沖一個矩陣快速冪就行了,
然而有一個邊權不一定為一的限制,所以暴力把點拆進一個矩陣,就可以只用矩陣快速冪來做這道題了,
代碼
點擊查看代碼
const ll N=110,P=2009,inf=1ll<<40;
ll n,m,t,g[N][N];
class Mat{
public:
ll a[N][N];
inline ll* operator[](ll x){return a[x];}
inline void one(){_for(i,1,m)a[i][i]=1;return;}
inline Mat operator*(Mat ma){
Mat ans;
_for(i,1,m){
_for(j,1,m){
ans[i][j]=0;
_for(k,1,m)ans[i][j]=(ans[i][j]+a[i][k]*ma[k][j]%P)%P;
}
}
return ans;
}
}tu;
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline ll rsnt(){
char c=getchar();
while(!isdigit(c))c=getchar();
return (c^48);
}
inline Mat fpow(Mat ma,ll b){
Mat ans;ans.one();
while(b){
if(b&1)ans=ans*ma;
ma=ma*ma,b>>=1;
}
return ans;
}
inline void Zhuan(){
_for(i,1,n){
_for(j,1,9){
ll k=10*(i-1)+j;
tu[k][k+1]=1;
}
_for(j,1,n){
ll k1=10*(i-1);
ll k2=10*(j-1);
if(g[i][j])tu[k1+g[i][j]][k2+1]=1;
}
}
return;
}
inline void In(){
n=rnt(),t=rnt();
_for(i,1,n)_for(j,1,n)g[i][j]=rsnt();
m=10*n,Zhuan();
Mat ans=fpow(tu,t);
printf("%lld\n",ans[1][m-9]);
return;
}
}
佳佳的 Fibonacci
思路
題目背景有時候不是白給的,
這道題中,聯系題目背景可以發現原式可以化簡為:
\[\begin{aligned} T(n)&=(F_1+F_2+\cdots+F_n)+(F_2+F_3\cdots+F_n)+(F_3+F_4+\cdots+F_n)+\cdots+(F_{n-1}+F_n)+F_n\\ &=(S(n)-S(0))+(S(n)-S(1))+(S(n)-S(3))+\cdots+(S(n)-S(n-2))+(S(n)-S(n-1))\\ &=\sum_{i=1}^{n}S(n)-S(i-1) \end{aligned} \]通過題目背景可知:
可這(指 \(S(n)\))對佳佳來說還是小菜一碟,
那么我們就去算 \(S(n)\),好像有點斷章取義,
同時,\(S(n)=S(n-1)+F_n\),所以 \(S(n)=F_{n+2}-1\),
為啥 $F_{-1}=1$ 啊?
\[\begin{aligned} \because &F_n=F_{n-1}+F_{n-2}\\ \therefore &F_1=F_0+F_{-1}\\ &1=0+F_{-1}\\ \therefore &F_{-1}=1\\ &同理可得:\\ &F_{-2}=-1,F_{-3}=2,F_{-4}=-3,F_{-5}=5,\cdots,F_{-n}=(-1)^{n+1}F_n \end{aligned} \]然后往原式子里帶:
\[\begin{aligned} T(n)&=\sum_{i=1}^{n}S(n)-S(i-1)\\ &=\sum_{i=1}^{n}(F_{n+2}-1)-(F_{i+1}-1)\\ &=nF_{n+2}-\sum_{i=2}^{n+1}F_{i}\\ &=nF_{n+2}-(S(n+1)-1)\\ &=nF_{n+2}-F_{n+3}+2 \end{aligned} \]然后沖一個矩陣快速冪就行了,
代碼
點擊查看代碼
const ll N=110,inf=1ll<<40;
ll n,m;
class Mat{
public:
ll a[N][N];
inline ll* operator[](ll x){return a[x];}
inline void one(){a[1][1]=a[1][2]=1;}
inline Mat operator*(Mat ma){
Mat ans;
_for(i,1,2){
_for(j,1,2){
ans[i][j]=0;
_for(k,1,2)ans[i][j]=(ans[i][j]+a[i][k]*ma[k][j]%m)%m;
}
}
return ans;
}
};
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline Mat fpow(Mat ma,ll b){
Mat ans;ans.one();
while(b){
if(b&1)ans=ans*ma;
ma=ma*ma,b>>=1;
}
return ans;
}
inline void In(){
n=rnt(),m=rnt();
Mat mat;mat[1][1]=mat[1][2]=mat[2][1]=1;
Mat ans=fpow(mat,n+1);
printf("%lld\n",(n*ans[1][2]%m-ans[1][1]+m+2)%m);
return;
}
}
選拔隊員(不知道教練從哪里找的)
題意
選出若干個男生和若干多個女生(即男女生的數目隨便定)安排到機房內的 \(N\) 個位置上去,要求任意兩位女生不能相鄰(即任意兩個女生之間必須有至少一個男生),求方案數 \(\bmod{M}\),
思路
直接用排列推是不行的,嘗試寫個 \(\text{dp}\),
設 \(f_{i,0}\) 表示有 \(n\) 個人坐了上去,最后一個人是男; \(f_{i,1}\) 表示有 \(n\) 個人坐了上去,最后一個人是女,
初始狀態為 \(f_{i,0}=f_{i,1}=1\),顯然有遞推式:
\[\begin{aligned} f_{i,0}&=f_{i-1,0}+f_{i-1,1}\\ f_{i,1}&=f_{i-1,0} \end{aligned} \]用矩陣優化:
\[\left[\begin{matrix} f_{n,0}&f_{n,1} \end{matrix}\right] \times \left[\begin{matrix} 1&1\\ 1&0 \end{matrix}\right] = \left[\begin{matrix} f_{n+1,0}&f_{n+1,1} \end{matrix}\right] \]誒這玩意兒不就是斐波那契數列嗎?!
然后沖一個矩陣快速冪就行了,
代碼
點擊查看代碼
const ll N=110,inf=1ll<<40;
ll T,m,n;
class Mat{
public:
ll a[5][5];
inline ll* operator[](ll x){return a[x];}
inline Mat operator*(Mat ma){
Mat ans;
_for(i,1,2){
_for(j,1,2){
ans[i][j]=0;
_for(k,1,2)ans[i][j]=(ans[i][j]+a[i][k]*ma[k][j]%m)%m;
}
}
return ans;
}
};
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline Mat fpow(ll b){
Mat ans;ans[1][1]=ans[1][2]=1,ans[2][1]=ans[2][2]=0;
Mat ma;ma[1][1]=ma[1][2]=ma[2][1]=1,ma[2][2]=0;
while(b){
if(b&1)ans=ans*ma;
ma=ma*ma,b>>=1;
}
return ans;
}
inline void In(){
n=rnt();
printf("%lld\n",fpow(n)[1][1]%m);
return;
}
}
Tr A
思路
不能理解這道題出的為什么這么靠后,
沖一個矩陣快速冪就行了,
代碼
點擊查看代碼
const ll N=20,P=9973,inf=1ll<<40;
ll T,n,k;
class Mat{
public:
ll a[N][N];
inline ll* operator[](ll x){return a[x];}
inline void one(){_for(i,1,n)a[i][i]=1;return;}
inline void zero(){memset(a,0,sizeof(a));return;}
inline Mat operator*(Mat ma){
Mat ans;ans.zero();
_for(i,1,n){
_for(j,1,n){
ans[i][j]=0;
_for(k,1,n)ans[i][j]=(ans[i][j]+a[i][k]*ma[k][j]%P)%P;
}
}
return ans;
}
}a;
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline Mat fpow(Mat ma,ll b){
Mat ans;ans.zero();ans.one();
while(b){
if(b&1)ans=ans*ma;
ma=ma*ma,b>>=1;
}
return ans;
}
inline ll GetAnswer(Mat ma){
ll num=0;
_for(i,1,n)num=(num+ma[i][i]);
return num%P;
}
inline void In(){
n=rnt(),k=rnt();
a.zero();
_for(i,1,n)_for(j,1,n)a[i][j]=rnt();
Mat ans=fpow(a,k);
printf("%lld\n",GetAnswer(ans));
return;
}
}
本文來自博客園,作者:Keven-He,轉載請注明原文鏈接:https://www.cnblogs.com/Keven-He/p/Matrix.html
轉載請註明出處,本文鏈接:https://www.uj5u.com/qita/501088.html
標籤:其他
上一篇:leetcode 240. Search a 2D Matrix II 搜索二維矩陣 II(中等)
下一篇:API 網關簡述
