2. 程式人生 > 實用技巧 >矩陣乘法小記

矩陣乘法小記

阿新 • • 發佈：2020-10-22

目錄

移植於原csdn部落格

基礎？遞推

斐波那契，\(F_0=1,F_1=1,F_n=F_{n-2}+F_{n-1}(n\geq 2)\)
所以通過\(O(n)\)解出

矩陣乘法

\[C_{i,j}=\sum A_{i,k}\times B_{k,j} \]

矩陣乘法是具有結合律的，所以通過這個特點可以用到快速冪。

  • 狀態矩陣，長度為\(n\)的一維向量
  • 轉移矩陣，與狀態矩陣相乘的矩陣

洛谷 1962 斐波那契數列

分析

容易可得，狀態矩陣就是\(\left[\begin{matrix}F_{n-1}&F_n\end{matrix}\right]\)
根據斐波那契數列，可以得到

\[\left[\begin{matrix}0&1\\ 1&1\end{matrix}\right] \]

顯然，對於\(F_{n-1}\)在下一時刻會給\(F_{n}\)，而\(F_{n}\)會結合\(F_{n-1}\)在下一時刻變成\(F_{n+1}\)

,所以得到了這個轉移矩陣
那麼

\[F_{n}=F_{0}\times \left[\begin{matrix}0&1\\ 1&1\end{matrix}\right]^n \]

HDU 5015 233 Matrix

題目

給定\(a_{1,0},a_{2,0},\dots,a_{n,0}\)，然後

\[a_{0,0}=0,a_{0,1}=233,a_{0,2}=2333,\dots \]

,規定\(a_{i,j}=a_{i-1,j}+a_{i,j-1}\)
詢問\(a_{n,m} \bmod 10000007\)

分析

首先想要維護字首和，但是233怎麼辦，可以假使\(a_{0,0}=23\)，那麼之後也就是\(a_{0,i}=a_{0,i-1}\times 10+3\)


那麼狀態矩陣是

\[\left[ \begin{matrix} 23&a_{1,0}&a_{2,0}&3 \end{matrix} \right] \]

想到轉移矩陣是

\[\left[ \begin{matrix} 10&0&0&1\\ 10&1&0&1\\ 10&1&1&1\\ 0&0&0&1 \end{matrix} \right] \]

所以

\[F_{m}=F_{0}\times \left[\begin{matrix}10&0&0&1\\10&1&0&1\\10&1&1&1\\0&0&0&1\end{matrix}\right]^m \]

但是如何表示矩陣中的1呢，也就是

\[\left[ \begin{matrix} 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1 \end{matrix} \right] \]

程式碼

#include <cstdio>
#include <cctype>
#include <cstring>
#define rr register
using namespace std;
struct maix{int p[13][13];}A,ANS;
const int mod=10000007; int n;
inline signed iut(){
    rr int ans=0; rr char c=getchar();
    while (!isdigit(c)) c=getchar();
    while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
    return ans;
}
inline maix mul(maix A,maix B){
    rr maix C;
    for (rr int i=0;i<n+2;++i)
    for (rr int j=0;j<n+2;++j) C.p[i][j]=0;
    for (rr int i=0;i<n+2;++i)
    for (rr int j=0;j<n+2;++j)
    for (rr int k=0;k<n+2;++k)
    C.p[i][j]=(1ll*A.p[i][k]*B.p[k][j]+C.p[i][j])%mod;
    return C;
}
signed main(){
    while (scanf("%d",&n)==1){
        for (rr int i=0;i<n+2;++i)
        for (rr int j=0;j<n+2;++j)
            ANS.p[i][j]=i==j;
        for (rr int i=0;i<n+2;++i)
        for (rr int j=0;j<n+2;++j){
            if (!j&&i!=n+1) A.p[i][j]=10;
            else if ((i>=j&&i!=n+1)||j==n+1) A.p[i][j]=1;
                else A.p[i][j]=0;
        }
        for (rr int m=iut();m;m>>=1,A=mul(A,A))
            if (m&1) ANS=mul(ANS,A);
        rr int ans=23ll*ANS.p[n][0];
        for (rr int i=1;i<=n;++i)
            ans=(1ll*iut()*ANS.p[n][i]+ans)%mod;
        ans=(3ll*ANS.p[n][n+1]+ans)%mod;
        printf("%d\n",ans);
    }
    return 0;
}

進階？？？

洛谷 4159 迷路

題目

windy在有向圖中迷路了。 該有向圖有\(N\)個節點，windy從節點1出發，他必須恰好在\(T\)時刻到達節點\(N\)。 現在給出該有向圖，你能告訴windy總共有多少種不同的路徑嗎？ 注意：windy不能在某個節點逗留，且通過某有向邊的時間嚴格為給定的時間。

分析

首先邊權是0或者1的時候，很容易想到矩陣乘法，然而邊權是\(0\sim 9\)，考慮拆點，把一個點拆成9個點，再建立轉移矩陣，對於每個拆的點，連向前一個拆的點，對於連邊，用原來的點連向邊權為\(w\)的拆出的點，嗯，就是這個樣子

程式碼

#include <cstdio>
#define rr register
using namespace std;
struct maix{int p[91][91];}A,ANS;
int n,m,t;
inline maix mul(maix A,maix B){
    rr maix C;
    for (rr int i=1;i<=t;++i)
    for (rr int j=1;j<=t;++j) C.p[i][j]=0;
    for (rr int i=1;i<=t;++i)
    for (rr int j=1;j<=t;++j)
    for (rr int k=1;k<=t;++k)
    C.p[i][j]=(A.p[i][k]*B.p[k][j]+C.p[i][j])%2009;
    return C;
}
signed main(){
    scanf("%d%d",&n,&m); t=n*9;
    for (rr int i=1;i<=n;++i){
        for (rr int j=1;j<9;++j) A.p[i+j*n][i+j*n-n]=1;
        for (rr int j=1,x;j<=n;++j){
            scanf("%1d",&x);
            if (x) A.p[i][j+x*n-n]=1;
        }
    }
    for (rr int i=1;i<=t;++i) ANS.p[i][i]=1;
    for (;m;m>>=1,A=mul(A,A))
        if (m&1) ANS=mul(ANS,A);
    return !printf("%d",ANS.p[1][n]);
}

洛谷 4834 薩塔尼亞的期末考試

題目

有\(n\)個點，假設第1個點被選的概率為\(P\),那麼第\(i\)個點就是\(iP\),釋放的能量為\(Fib(i)\)，問期望釋放的能量

分析

根據概率期望，可以得到答案也就是求

\[\frac{\sum_{i=1}^{n}i\times Fib(i)}{n(n+1)\div2} \]

下面用逆元就可以解決，但是上面嘗試化簡
首先要

\[\sum_{i=1}^n fib(i)=fib(1)+fib(2)+\dots +fib(n-1)+fib(n) \]

\[=1+fib(1)+fib(2)+\dots+fib(n-1)+fib(n)-1 \]

\[=2fib(2)+fib(1)+\dots+fib(n-1)+fib(n)-1=\dots=fib(n+2)-1 \]

那麼

\[\sum_{i=1}^n i\times fib(i)=n\sum_{i=1}^nfib(i)-\sum_{i=1}^{n-1}\sum_{j=1}^{i}fib(j) \]

\[=nfib(n+2)-n-\sum_{i=1}^{n-1}fib(i+2)-1=nfib(n+2)-1-\sum_{i=1}^{n-1}fib(i+2) \]

\[=nfib(n+2)+1-\sum_{i=1}^{n+1}fib(i)=nfib(n+2)+1-(fib(n+3)-1)=nfib(n+2)-fib(n+3)+2 \]

剩下的就是一個斐波那契數列解決的問題了，程式碼就不貼了

洛谷 5107 能量採集

題目

求\(\sum_{i=1}^n\sum_{j=1}^mgcd(i,j)\)

好的並不是這樣
給定一個\(n\)個點\(m\)條邊的有向圖，每個點有初始能量\(a_i\)。

每過一秒，每個點的能量便會等量地流向所有出邊，另外，會有一份流向自己（你可以當做有一個自環）。

現在dkw有\(q\)次詢問，每次詢問會給你一個時間\(t\) ，dkw想知道\(t\)秒時每個點的能量。

不保證圖中沒有重邊和自環，答案對998244353取模。

分析

可以發現這就是一個矩陣乘法，按照題目建就可以了，時間複雜度\(O(qn^3logt)\)，

然而會TLE，怎麼辦，考慮二進位制拆分，也就是預處理\(2^x\)的答案，

然後拼湊答案，那麼時間複雜度是\(O((n^3+qn^2)logt)\)，雖然這個時間複雜度好像很難辦，但至少可以過。

程式碼

50分警告還是不可以過？？發現矩陣乘法對於詢問時顯然是浪費了，於是就有了

還是T掉了，考慮取模優化

還是T掉了，考慮把long long改成int

#include <cstdio>
#include <cctype>
#define rr register
using namespace std;
const int mod=998244353; int n,m,q;
struct maix{int p[51][51];}A[31],T;
inline signed iut(){
    rr int ans=0; rr char c=getchar();
    while (!isdigit(c)) c=getchar();
    while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
    return ans;
}
inline signed ksm(int x,int y){
    rr int ans=1;
    for (;y;y>>=1,x=1ll*x*x%mod)
        if (y&1) ans=1ll*ans*x%mod;
    return ans; 
}
inline signed mo(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline maix mul(maix A,maix B,int t){
    rr maix C;
    for (rr int i=0;i<t;++i)
    for (rr int j=0;j<n;++j) C.p[i][j]=0;
    for (rr int i=0;i<t;++i)
    for (rr int j=0;j<n;++j)
    for (rr int k=0;k<n;++k)
    C.p[i][j]=mo(C.p[i][j],1ll*A.p[i][k]*B.p[k][j]%mod);
    return C;
}
signed main(){
    n=iut(); m=iut(); q=iut();
    for (rr int i=0;i<n;++i) T.p[0][i]=iut(),A[0].p[i][i]=1;
    while (m--) ++A[0].p[iut()-1][iut()-1];
    for (rr int i=0;i<n;++i){
        rr int sum=0;
        for (rr int j=0;j<n;++j) sum+=A[0].p[i][j];
        sum=ksm(sum,mod-2);
        for (rr int j=0;j<n;++j) A[0].p[i][j]=1ll*A[0].p[i][j]*sum%mod;
    }
    for (rr int i=1;i<30;++i) A[i]=mul(A[i-1],A[i-1],n);
    while (q--){
        rr int t=iut(); rr maix ANS=T;
        for (rr int i=0;i<30;++i)
            if ((t>>i)&1) ANS=mul(ANS,A[i],1);
        rr int ans=0;
        for (rr int i=0;i<n;++i) ans^=ANS.p[0][i];
        printf("%d\n",ans>=mod?ans-mod:ans);
    }
    return 0;
}

洛谷 3216 數學作業

分析

根據dp，可以知道

\[F_n=F_{n-1}+10^{\lfloor \log n\rfloor+1}+n \]

那麼轉移矩陣即為

\[\left[\begin{matrix}10^k&0&0\\1&1&0\\1&1&1\end{matrix}\right] \]

然而\(10^k\)比較難搞，考慮分塊，列舉\(k\)把它分成\(k\)個部分依次相乘

程式碼

#include <cstdio>
#include <cstring>
#define rr register
using namespace std;
typedef long long ll;
struct maix{ll p[3][3];}A,ANS;
ll n; int mod;
inline maix mul(maix A,maix B){
    rr maix C; memset(C.p,0,sizeof(C.p));
    for (rr int i=0;i<3;++i)
    for (rr int j=0;j<3;++j)
    for (rr int k=0;k<3;++k)
    C.p[i][j]=(A.p[i][k]*B.p[k][j]+C.p[i][j])%mod;
    return C;
}
signed main(){
    scanf("%lld%d",&n,&mod);
    ANS.p[0][0]=ANS.p[1][1]=ANS.p[2][2]=1;
    for (rr ll t=10,k;;t=t*10){
        if (t>n) k=n-t/10+1; else k=t-t/10;
        memset(A.p,0,sizeof(A.p));
        A.p[0][0]=t%mod; A.p[1][0]=1; A.p[2][0]=1; A.p[1][1]=1; A.p[2][1]=1; A.p[2][2]=1;
        for (;k;k>>=1,A=mul(A,A))
            if (k&1) ANS=mul(ANS,A);
        if (t>n) break;
    }
    return !printf("%lld\n",ANS.p[2][0]);
}

洛谷 5110 塊速遞推

題目

\(a_n=233*a_{n-1}+666*a_{n-2},a_0=0,a_1=1\),多組資料求\(a_n\)

分析

特徵方程做法
可以輕鬆列出轉移矩陣是

\[\left[\begin{matrix}233&666\\1&0\end{matrix}\right] \]

但是顯然是TLE的\(O(Tn^3logm)\)
二進位制優化後是\(O(n^3logm+Tn^2logm)\)
但是\(T\)太大了，可能還是會T，至少我沒打這種方法
因為\(a^b=a^{b\bmod \varphi(p)}(\bmod p)[gcd(a,p)==1]\),所以可以把時間複雜度優化到\(O(n^3logmod+Tn^2logmod)\)
PS：這題是因為出題人模數設定的好，不然不能這麼亂用尤拉定理
然而這樣還是會TLE，嘗試用矩陣分塊，設\(k=\lceil\sqrt{mod}\rceil\)
那麼其實\(a^b=(a^k)^{\lfloor\frac{b}{k}\rfloor}+a^{b\bmod k}\)，兩個東西都可以預處理，從而把時間複雜度優化到\(O(n^3\sqrt k+Tn^3)\)，然後迴圈展開就可以草率的理解為\(O(\sqrt k+T)\)

程式碼

#include <cstdio>
#include <cstring>
#define rr register
using namespace std;
const int N=31623,mod=1000000007;
typedef unsigned long long ull;
struct maix{int p[2][2];}A[N+1],B[N+1];
inline signed mo(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline maix mul(maix A,maix B){
    rr maix C;
    C.p[0][0]=mo(1ll*A.p[0][0]*B.p[0][0]%mod,1ll*A.p[0][1]*B.p[1][0]%mod);
    C.p[0][1]=mo(1ll*A.p[0][0]*B.p[0][1]%mod,1ll*A.p[0][1]*B.p[1][1]%mod);
    C.p[1][0]=mo(1ll*A.p[1][0]*B.p[0][0]%mod,1ll*A.p[1][1]*B.p[1][0]%mod);
    C.p[1][1]=mo(1ll*A.p[1][0]*B.p[0][1]%mod,1ll*A.p[1][1]*B.p[1][1]%mod);
    return C;
}
signed main(){
    A[1].p[0][0]=233,A[1].p[0][1]=666,A[0].p[0][0]=A[0].p[1][1]=A[1].p[1][0]=1,B[0]=A[0];
    for (rr int i=2;i<=N;++i) A[i]=mul(A[i-1],A[1]);
    for (rr int i=1;i<=N;++i) B[i]=mul(B[i-1],A[N]);
    rr int t,n,ans=0; rr ull sa,sb,sc,temp; scanf("%d",&t);
    for (scanf("%llu %llu %llu",&sa,&sb,&sc);t;--t){
        sa^=sa<<32,sa^=sa>>13,sa^=sa<<1,temp=sa,sa=sb,sb=sc,sc^=temp^sa;
        n=(sc-1)%(mod-1); rr maix C=mul(B[n/N],A[n%N]); ans^=C.p[0][0];
    }
    return !printf("%d",ans);
}

POJ 3233 Matrix Power Series

題目

給定一個矩陣\(A\)，求\(\sum_{i=1}^kA_k\)

分析

這道題按照藍書的版本應該是分治，但是貌似有點慢，考慮找到轉移矩陣，根據dalao的解法

\[\left[\begin{matrix}A&E\\0&E\end{matrix}\right] \]

即為轉移矩陣，雖然解釋可能會有點問題，但是驚人的發現

\[\left[\begin{matrix}A&E\\0&E\end{matrix}\right]^k=\left[\begin{matrix}A&E+\sum_{i=1}^kA_k\\0&E\end{matrix}\right] \]

所以最後還要減掉一個單位矩陣

程式碼

#include <cstdio>
#include <cctype>
#include <cstring>
#define rr register
using namespace std;
struct maix{int p[61][61];}A,ANS; int n,mod,k;
inline signed iut(){
    rr int ans=0; rr char c=getchar();
    while (!isdigit(c)) c=getchar();
    while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
    return ans;
}
inline maix mul(maix A,maix B){
    rr maix C; memset(C.p,0,sizeof(C.p));
    for (rr int i=1;i<=n;++i)
    for (rr int j=1;j<=n;++j)
    for (rr int k=1;k<=n;++k)
    C.p[i][j]=(A.p[i][k]*B.p[k][j]+C.p[i][j])%mod;
    return C;
}
signed main(){
	n=iut(); k=iut()+1; mod=iut();
	for (rr int i=1;i<=n;++i)
	for (rr int j=1;j<=n;++j) A.p[i][j]=iut();
	for (rr int i=1;i<=n;++i) A.p[i+n][i+n]=A.p[i][i+n]=1;
	n<<=1; for (rr int i=1;i<=n;++i) ANS.p[i][i]=1;
	for (;k;k>>=1,A=mul(A,A)) if (k&1) ANS=mul(ANS,A);
	n>>=1; for (rr int i=1;i<=n;++i) --ANS.p[i][i+n];
	for (rr int i=1;i<=n;++i)
	for (rr int j=1;j<=n;++j)
	    printf("%d%c",ANS.p[i][j+n]<0?ANS.p[i][j+n]+mod:ANS.p[i][j+n],j==n?10:32);
	return 0;
}

後續

如何優化卡常
減少取模次數，少用long long，減少不必要的迴圈，減少矩陣大小,迴圈展開
對於只有單個數據，時間複雜度顯然是\(O(n^3logm)\)，空間複雜度\(O(n^2)\)
對於多組資料，時間複雜度顯然是\(O(Tn^3logm)\)，那麼需要通過預處理減少詢問時的時間
可以二進位制拆分，也就是預處理\(2^k\)次冪的結果再拼湊，時間複雜度為\(O(n^3logm+Tn^2logm)\)，空間複雜度\(O(n^2logm)\)
當然還有矩陣分塊，大致就是\(a^b=a^{b \bmod r}\times (a^r)^{\lfloor\frac{b}{r}\rfloor}\)，所以時間複雜度應該是\(O(n^3r+Tn^3)\)。這個\(r\)就非常關鍵了
然而空間相應會耗損較大,為\(O(n^2r)\)，所以要權衡利弊，二進位制拆分會比較平衡，但也非常容易被卡，貌似還有\(k\)進位制拆分這種操作，然而我啥也不會
完結，撒花

相關推薦

矩陣乘法小記

目錄基礎？遞推矩陣乘法洛谷 1962 斐波那契數列分析HDU 5015 233 Matrix題目分析程式碼洛谷 4159 迷路題目分析程式碼洛谷 4834 薩塔尼亞的期末考試題目分析洛谷 5107 能量採集題目分析程式碼洛谷 3216 數學作業分析程

[HNOI2008] GT考試 - KMP,矩陣乘法,dp

Description 准考證號為 \\(N\\) 位數 \\(X_1,X_2…X_n(0\\le X_i\\le9)\\)，他不希望准考證號上出現不吉利的數字。他的不吉利數 \\(A_1,A_2…A_m(0\\le A_i\\le 9)\\) 有 \\(M\\) 位，不出現是指 \\(X_1,X_2…X_n\\

【學習筆記】矩陣乘法題目選講

這是校內講課課件。其實這些程式碼我覺得並沒有必要放上，不過作為許多題解的合集，海蝕放上了。

CF1380F.Strange Addition（線段樹維護矩陣乘法）

題意： 定義一種特殊的加法，這種加法沒有進位，其他與普通加法一樣，比如16+16=212。

快速冪+龜速乘+費馬小定理+逆元+矩陣乘法

我是這個機房最菜的 我今天覆習的是： 王者吃雞CF，上分小隊等你來 扯遠了，接下來才是乾貨

LOJ#3058. 「HNOI2019」白兔之舞 單位根反演+矩陣乘法+MTT

毒瘤多項式. 找原根+矩陣乘法+拆係數FFT大雜燴. code: #include <vector> #include <cmath>

矩陣乘法(HJ69)

一：解題思路 二：完整程式碼示例 (C++版和Java版) C++程式碼： #include <iostream>

#矩陣乘法，斐波那契#洛谷 2544 [AHOI2004] 數字迷陣

題目 分析 oeis找規律得到第一列和第二列的通項公式，然後矩陣乘法 程式碼 #include <cstdio>

牛客-拆分（矩陣乘法）

題目大意 有\\(n\\)個數，現在把它們劃分成若干段，使存在一段的數的個數不為給出的\\(k\\)個數中的任意一個。問有多少種劃分方式

【題解】CF645E. Intellectual Inquiry / sequence【20201020 CSP 模擬賽】【貪心 矩陣乘法】

題目連結 題目連結（略有不同） 題意 有一長為 \\(n\\) 的序列 \\(a\\)，\\(a_i\\in [1,k]\\cap \\mathbf{N}\\)。你要在後面添 \\(m\\) 個數（\\(a_i\\in [1,k]\\cap \\mathbf{N}\\)），使得新序列的本質不同子序列個

矩陣乘法優化DP複習

矩陣乘法優化DP複習 前言 最近做毒瘤做多了……聯賽難度的東西也該複習複習了。

#矩陣乘法#洛谷 5343 【XR-1】分塊

題目 分析 考慮dp，\\(dp[i]=\\sum dp[i-j]\\) 既然\\(j\\)很小，那麼這顯然可以用矩陣乘法優化

洛谷1527 [國家集訓隊]矩陣乘法（整體二分）

題意：給你一個的矩陣，每次詢問一個子矩形的第k小數 輸入格式：第一行有兩個整數，分別表示矩陣大小n和詢問組數q

【LOJ NOI Round#2 Day1 T1】單槍匹馬（矩陣乘法）

題目傳送門 操作二要求的東西是一個迴圈迭代的東西，手推相鄰兩項找下規律，發現相鄰兩項的分子分母間含有線性關係，考慮用矩陣乘法求解。對於 \\([1,n]\\)的詢問，從後往前倒推， \\(x_{n-1}=a_{n-1} \\times x_{n}

學習筆記--矩陣乘法

定義 數有數的乘法，矩陣乘上矩陣就叫做矩陣乘法 在數學上，我們這樣定義矩陣乘法

torch中向量、矩陣乘法大總結

技術標籤：# pytorchpython人工智慧ptorch LawsonAbs的認知與思考，還請各位讀者批判閱讀。

python numpy for迴圈矩陣乘法

技術標籤：# 矩陣python線性代數numpy 把向量放到一個list中，用for i in list進行遍歷

python numpy中array和mat型別矩陣乘法的三種方式比較（np.dot()、*、np.multiply()）

技術標籤：# 矩陣pythonnumpy線性代數 import numpy as np # np.array() # 1. np.dot() 矩陣乘法

例5：求斐波那契數列變形前n項的和【矩陣乘法】【快速冪】【斐波那契數列】

技術標籤：題解矩陣乘法快速冪斐波那契數列 題目描述 數列 f [ n ] = f [ n − 1 ] + f [ n − 2 ]

Matrix Power Series【矩陣乘法】

技術標籤：矩陣乘法快速冪 Time Limit: 3000MS Memory Limit: 131072K Total Submissions: 35006 Accepted: 14029