矩陣加速線性遞推

zuoqingyuan111發表於2024-08-19
矩陣

前置：【模板】矩陣快速冪

斐波那契數列

定義 \(Fib_i\) 為斐波那契數列第 \(i\) 項，斐波那契數列定義如下

\[Fib_n = \left\{\begin{aligned} 1 \space (n \le 2) \\ Fib_{n-1}+Fib_{n-2} \space (n\ge 3) \end{aligned}\right. \]

給定一個 \(n\)，求出 \(Fib_n \mod 10^9+7\) 的值。其中 \(n \le 10^9\)。

此題的一個經典做法是遞推，時間複雜度為 \(O(n)\)，當 \(n \le 10^9\) 時無法快速得到答案。而下面則會介紹一種能夠在 \(O(\log n)\) 的時間複雜度內求解的辦法。

定義一個矩陣 \(F_n=\left[\begin{matrix} Fib_n&Fib_{n+1}\\ \end{matrix} \right]\)，
定義一個矩陣 \(A=\left[\begin{matrix} 0&1\\ 1&1 \end{matrix} \right]\)。
我們就會發現

\[\begin{aligned} F_n\times A&=\left[\begin{matrix} Fib_n\times 0+Fib_{n+1}\times1&Fib_n\times 1+Fib_{n+1}\times 1\\ \end{matrix} \right]\\ &=\left[\begin{matrix} Fib_{n+1}&Fib_n+Fib_{n+1}\\ \end{matrix} \right]\\ &=\left[\begin{matrix} Fib_{n+1}&Fib_{n+2}\\ \end{matrix} \right]\\ &=F_{n+1} \end{aligned}\]

也就是說，\(F_n\times A\) 相當於往前遞推一次為 \(F_{n+1}\)，根據定義，我們讓 \(F_0=\left[\begin{matrix} 0&1\\ \end{matrix} \right]\)，只要計算出 \(F_0\times A^n\)，就能得到 \(Fib_n\) 的值，而 \(A^n\) 則可以透過快速冪來求解。這個問題就可以在 \(O(\log n)\) 內求解。

關於矩陣乘法：
計算兩個矩陣的乘法。\(n \times m\) 階的矩陣 \(A\) 乘以 \(m \times k\) 階的矩陣 \(B\) 得到的矩陣 \(C\) 是 \(n \times k\) 階的，且 \(C_{i,j}=\sum\limits_{k=1}^m A_{i,k}\times B_{k,j}\)

題目在這裡：斐波那契數列（矩陣加速）

例題 \(1\)

P1939 矩陣加速（數列）

做法還是矩陣加速遞推，定義 \(F_n=\left[\begin{matrix} a_n&a_{n+1}&a_{n+2}\\ \end{matrix} \right]\)。重點在與如何構造矩陣 \(A\)，使得 \(F_n \times A=F_{n+1}\)。

\(F_n,F_{n+1}\) 的矩陣都是 \(1\times 3\) 大小，根據 \(n \times m\) 階的矩陣 \(A\) 乘以 \(m \times k\) 階的矩陣 \(B\) 得到的矩陣 \(C\) 是 \(n \times k\) 階的。所以矩陣 \(A\) 是一個 \(3\times 3\) 階的矩陣。

\(F_{n+1}\) 矩陣的第一項為 \(a_{n+1}\)，也就是 \(a_n\times 0+a_{n+1}\times 1+a_{n+2}\times 0\)，所以 \(A\) 矩陣第一維是 \(\left[\begin{matrix} 0&1&0\\ \end{matrix} \right]\) .同理，第二維是 \(\left[\begin{matrix} 0&0&1\\ \end{matrix} \right]\)。

\(F_n\) 的第三個數字是 \(a_{n+3}=a_{n+2}+a_{n}=a_n\times 1+a_{n+1}\times 0+a_{n+2}\times 1\)。那麼矩陣 \(A=\left[\begin{matrix} 0&1&0\\ 0&0&1\\ 1&0&1\\ \end{matrix} \right]\)。

另外千萬記得在寫矩陣乘法是一定要記得初始化。

點選檢視程式碼 
#include <iostream>
#include <cstring>
#define int long long
using namespace std;
const long long mod=1e9+7;
int f[3],a[3][3],n,q;
void cheng1(){//矩陣F*矩陣A
    int c[3];
    memset(c,0,sizeof(c));
    for(int i=0;i<=2;i++){
        for(int j=0;j<=2;j++){
            c[i]=(long long)(c[i]+f[j]*a[i][j])%mod;
        }
    }
    for(int i=0;i<=2;i++)f[i]=c[i];
    return;
}
void cheng2(){//矩陣A*矩陣A
    int b[3][3];
    memset(b,0,sizeof(b));
    for(int i=0;i<=2;i++){
        for(int j=0;j<=2;j++){
            for(int k=0;k<=2;k++){
                b[i][j]=(long long)(b[i][j]+a[i][k]*a[k][j])%mod;
            }
        }
    }
    for(int i=0;i<=2;i++)for(int j=0;j<=2;j++)a[i][j]=b[i][j];
    return;
}
void ksm(int b){//ksm
    while(b){
        if(b&1)cheng1();
        cheng2();
        b>>=1;
    }
    return;
}
signed main(){
    cin>>q;
    while(q--){
        cin>>n;
        a[0][0]=0,a[0][1]=1,a[0][2]=0,a[1][0]=0,a[1][1]=0,a[1][2]=1,a[2][0]=1,a[2][1]=0,a[2][2]=1,f[1]=1,f[2]=1,f[0]=0;
        ksm(n);
        cout<<f[0]<<endl;
    }
    return 0;
}

可以看出矩陣乘法最佳化線性遞推的特點是，遞推的每一項只會與前面若干項取值相關。矩陣乘法的遞推難點主要在於構建矩陣 \(A\)。且大部分題目中矩陣 \(A\) 的值都要視情況而定。

另外矩陣快速冪的時間複雜度不是標準的 \(O(\log n)\)，而是附帶一個常數，就是矩陣 \(F\) 第二維的三次方，如斐波那契數列問題中，\(F_n\) 的第二維長度為 \(2\)，則時間複雜度是 \(O(2^3\log n)\)；數列加速問題中，\(F_n\) 的第二維長度為 \(3\)，則時間複雜度是 \(O(3^3\log n)\)。因此，如果第二維長度過大，則程式也會超時。

例題 \(2\)

廣義斐波那契數列

矩陣 \(F\) 同普通斐波那契，矩陣 \(A=\left[\begin{matrix} 0&1\\ q&p \end{matrix} \right],F_0=\left[\begin{matrix} a_1&a_2\\ \end{matrix} \right]\)

很多遞推的矩陣 \(A\) 要視情況而定。

例題 \(3\)

[NOI2012] 隨機數生成器

\(F_n=\left[\begin{matrix} X_n&1 \end{matrix} \right],A=\left[\begin{matrix} a&c\\ 0&1\\ \end{matrix} \right],F_0=\left[\begin{matrix} X_0&1 \end{matrix} \right]\)

對於這種在遞推中要加一個常數項的情況，在矩陣 \(F\) 中多加一維為 \(1\)，然後 \(A\) 中加上一個係數，表示 \(+（1\times c）\)。

例題 \(4\)

[SDOI2008] 遞迴數列

本題的目標是 \(\sum\limits_{i=m}^{n}{a_{i}} \mod p\)，我們可以分別求出 \(\sum\limits_{i=1}^{m-1}a_i,\sum\limits_{i=1}^{n}a_i\) 的值，相減後再 \(\mod p\)。問題就在於如何求出\(\sum\limits_{i=1}^{n}a_i\) 的值

我們另 \(S_n=\sum\limits_{i=1}^{n}a_i\)，構建一個 \(1\times (k+1)\) 的矩陣 \(F\)，其中 \(F_n=\left[\begin{matrix} S_n&S_{n+1}&\dots&S_{n+k} \end{matrix} \right],F_0=\left[\begin{matrix} 0&b_1&b_1+b_2&\dots&\sum\limits_{i=1}^k b_i \end{matrix} \right]\)

矩陣 \(A\) 的前 \(k\) 行很好構建，\(A_{i,i+1}=1(i<=k)\)，其他均為 \(0\)，但是第 \(k+1\) 行卻很難表示。

我們不妨思考一下，\(S_{n+k+1}\)如何表示

\[\begin{aligned}S_{n+k+1}&=S_{n+k}+a_{n+k+1}\\ &=a_{n+k+1}=\sum_{j=1}^{k}{c_{j} \times a_{n+k+1-j}}\\ &=a_{n+k+1}=\sum_{j=1}^{k}{c_{j} \times (S_{n+k+1-j}-S_{n+k-j})}\\ &=S_{n+k+1}=S_{n+k}+\sum_{j=1}^{k}{c_{j}\times S_{n+k+1-j}-c_j\times S_{n+k-j}}\\ &=S_{n+k+1}=S_{n+k}\times(c_k+1)+\sum_{j=1}^{k-1}S_{n+k-j}\times(c_{k-j+1}-c_{k-j})-S_n\times c_k \end{aligned} \]

上面式子太抽象，舉個例子：

當k=2時
F[n]=[s[n],s[n+1],s[n+2]]
s[n+3]=s[n+2]+a[n+3]
a[n+3]=c[1]*a[n+2]+c[2]*a[n+1]
a[n+3]=c[1]*(s[n+2]-s[n+1])+c[2]*(s[n+1]-s[n])
a[n+3]=c[1]*s[n+2]-c[1]*s[n+1]+c[2]*s[n+1]-c[2]*s[n]
s[n+3]=s[n+2]*(c[1]+1)+s[n+1]*(c[2]-c[1])+s[n]*(-c[2])

把他們拆開，計算每一個 \(S\) 的係數，購造出矩陣 \(A\) 的第 \(k+1\) 行

  1. \(A_{k+1,k+1}=c_1+1\)
  2. \(A_{k+1,1}=-c_k\)
  3. 對於剩下的 \(A_{k+1,i}\)，\(A_{k+1,i}=c_{k-i+2}-c_{k-i+1}\)

按照上述思路，就可以寫出矩陣 \(A\)。

\[\left[ \begin{matrix} 0&1&0&\dots&0\\ 0&0&1&\dots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\dots&1\\ -c_k&c_k-c_{k-1}&c_{k-1}-c_{k-2}&\dots&c_1+1 \end{matrix} \right]\]

我們以樣例為例

0  1 0
0  0 1
-1 0 2

這就是樣例的矩陣 \(A\)。

最後的程式碼在過程中不斷取模就可以了

點選檢視程式碼 
#include <iostream>
#include <cstdio>
#include <cstring> 
using namespace std;
typedef long long ll;
ll k,n,m,b[25],c[25],mod,f[25],a[25][25];
inline ll read(){
    ll r=0;char ch=getchar();
    while(ch<'0'||ch>'9')ch=getchar();
    while(ch>='0'&&ch<='9')r=(r<<1)+(r<<3)+(ch^48),ch=getchar();
    return r;
}
void write(ll r){
    if(r>9)write(r/10);
    putchar(r%10+'0');
    return;
}
void init(){
	memset(a,0,sizeof(a));
	memset(f,0,sizeof(f));
	for(int i=2;i<=k+1;i++)f[i]=b[i-1]+f[i-1];
	for(int i=1;i<=k;i++)a[i][i+1]=1;
	for(int i=1;i<=k+1;i++){
		if(i==1)a[k+1][i]=-c[k];
		else if(i==k+1)a[k+1][i]=c[1]+1;
		else a[k+1][i]=c[k-i+2]-c[k-i+1];
	} 
	return;
}
//省略一部分
int main(){
	k=read();
	for(int i=1;i<=k;i++)b[i]=read();
	for(int i=1;i<=k;i++)c[i]=read();
	m=read(),n=read(),mod=read();
	write(((ksm(n)-ksm(m-1))%mod+mod)%mod);
	return 0;
}

蒟蒻第一個不看題解過掉的紫題

例題 \(5\)

[TJOI2019] 甲苯先生的字串

矩陣加速最佳化動態規劃。

在本題中，我們令 \(dp_{i,j}\) 表示已經選定前 \(i-1\) 個字元，第 \(i\) 個字元選 \(j\) 的總方案數。最終的答案就是 \(\sum dp_{n,j}\)。預處理出在 \(s_1\) 中字元 \(j\) 前面出現的字符集 \(C_j\)。那麼則有一個很明顯的轉移方程

\[dp_{i,j}=\sum\limits_{k\notin C_j}dp_{i,k}(i>1) \]

在 \(i=1\) 的初始情況下，\(dp_{1,j}=1\)。

注意到 \(n\le 10^{15}\)，字符集的大小隻有 \(26\)。考慮矩陣乘法。另 \(F_n=\left[\begin{matrix} dp_{n,1},dp_{n,2},\dots,dp_{n,26} \end{matrix}\right]\)。對於矩陣 \(A\)。如果 \(j\in C_i\)，則 \(A_{i,j}=0\)，否則為 \(A_{i,j}=1\)。

點選檢視程式碼 
#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
typedef long long ll;
const int mod=1e9+7,N=1e5+10;
ll f[30],a[30][30],n;
char s[N];
int ksm(ll b){
    ll cnt=0;
    while(b){
        if(b&1)cheng1();
        cheng2();
        b>>=1;
    }
    for(int i=1;i<=26;i++)cnt=(cnt+f[i])%mod;
    return cnt;
}
//省略一部分
int main(){
    scanf("%lld",&n);
    scanf("%s",s+1);
    int len=strlen(s+1);
    for(int i=1;i<=26;i++)for(int j=1;j<=26;j++)a[i][j]=1;
    for(int i=2;i<=len;i++){
        a[s[i]-'a'+1][s[i-1]-'a'+1]=0;
    }
    for(int i=1;i<=26;i++)f[i]=1;
    cout<<ksm(--n)<<endl;
    return 0;
}

小結

矩陣乘法加速遞推可以說是我第一個學明白的數論演算法了，雖然感覺在實戰中的使用率不大，但是還是感覺是一個很有用的演算法

相關文章