矩阵加速线性递推

前置：【模板】矩阵快速幂

斐波那契数列

定义 $Fi{b}_i$ 为斐波那契数列第 $i$ 项，斐波那契数列定义如下

$$Fi{b}_n=\{\begin{array}{r}1(n\le 2)\\ Fi{b}_{n-1}+Fi{b}_{n-2}(n\ge 3)\end{array}$$

给定一个 $n$，求出 $Fi{b}_{n}mod{10}^9+7$ 的值。其中 $n\le {10}^9$。

此题的一个经典做法是递推，时间复杂度为 $O(n)$，当 $n\le {10}^9$ 时无法快速得到答案。而下面则会介绍一种能够在 $O(\log n)$ 的时间复杂度内求解的办法。

定义一个矩阵 $F_n=\left[\begin{array}{cc}Fi{b}_n& Fi{b}_{n+1}\end{array}\right]$，
定义一个矩阵 $A=\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right]$。
我们就会发现

$$\begin{array}{rl}{F}_n\times A& =\left[\begin{array}{cc}Fi{b}_n\times 0+Fi{b}_{n+1}\times 1& Fi{b}_n\times 1+Fi{b}_{n+1}\times 1\end{array}\right]\\ & =\left[\begin{array}{cc}Fi{b}_{n+1}& Fi{b}_n+Fi{b}_{n+1}\end{array}\right]\\ & =\left[\begin{array}{cc}Fi{b}_{n+1}& Fi{b}_{n+2}\end{array}\right]\\ & =F_{n+1}\end{array}$$

也就是说，$F_n\times A$ 相当于往前递推一次为 $F_{n+1}$，根据定义，我们让 $F_0=\left[\begin{array}{cc}0& 1\end{array}\right]$，只要计算出 $F_0\times A^n$，就能得到 $Fi{b}_n$ 的值，而 $A^n$ 则可以通过快速幂来求解。这个问题就可以在 $O(\log n)$ 内求解。

关于矩阵乘法：
计算两个矩阵的乘法。$n\times m$ 阶的矩阵 $A$ 乘以 $m\times k$ 阶的矩阵 $B$ 得到的矩阵 $C$ 是 $n\times k$ 阶的，且 $C_{i,j}=\sum _{k=1}^m{A}_{i,k}\times B_{k,j}$

题目在这里：斐波那契数列（矩阵加速）

例题 $1$

P1939 矩阵加速（数列）

做法还是矩阵加速递推，定义 $F_n=\left[\begin{array}{ccc}{a}_n& a_{n+1}& a_{n+2}\end{array}\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{array}{ccc}0& 1& 0\end{array}\right]$ .同理，第二维是 $\left[\begin{array}{ccc}0& 0& 1\end{array}\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{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 1& 0& 1\end{array}\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{array}{cc}0& 1\\ q& p\end{array}\right],F_0=\left[\begin{array}{cc}{a}_1& a_2\end{array}\right]$

很多递推的矩阵 $A$ 要视情况而定。

例题 $3$

[NOI2012] 随机数生成器

$F_n=\left[\begin{array}{cc}{X}_n& 1\end{array}\right],A=\left[\begin{array}{cc}a& c\\ 0& 1\end{array}\right],F_0=\left[\begin{array}{cc}{X}_0& 1\end{array}\right]$

对于这种在递推中要加一个常数项的情况，在矩阵 $F$ 中多加一维为 $1$，然后 $A$ 中加上一个系数，表示 $+（1\times c）$。

例题 $4$

[SDOI2008] 递归数列

本题的目标是 $\sum _{i=m}^n{a}_i\mathrm{mod}p$，我们可以分别求出 $\sum _{i=1}^{m-1}{a}_i,\sum _{i=1}^n{a}_i$ 的值，相减后再 $\mathrm{mod}p$。问题就在于如何求出$\sum _{i=1}^n{a}_i$ 的值

我们另 $S_n=\sum _{i=1}^n{a}_i$，构建一个 $1\times (k+1)$ 的矩阵 $F$，其中 $F_n=\left[\begin{array}{cccc}{S}_n& S_{n+1}& \dots & S_{n+k}\end{array}\right],F_0=\left[\begin{array}{ccccc}0& b_1& b_1+b_2& \dots & \sum _{i=1}^k{b}_i\end{array}\right]$

矩阵 $A$ 的前 $k$ 行很好构建，$A_{i,i+1}=1(i<=k)$，其他均为 $0$，但是第 $k+1$ 行却很难表示。

我们不妨思考一下，$S_{n+k+1}$如何表示

$$\begin{array}{rl}{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{array}$$

上面式子太抽象，举个例子：

当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{array}{ccccc}0& 1& 0& \dots & 0\\ 0& 0& 1& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & 1\\ -c_k& c_k-c_{k-1}& c_{k-1}-c_{k-2}& \dots & c_1+1\end{array}\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] 甲苯先生的字符串

矩阵加速优化动态规划。

在本题中，我们令 $d{p}_{i,j}$ 表示已经选定前 $i-1$ 个字符，第 $i$ 个字符选 $j$ 的总方案数。最终的答案就是 $\sum d{p}_{n,j}$。预处理出在 $s_1$ 中字符 $j$ 前面出现的字符集 $C_j$。那么则有一个很明显的转移方程

$$d{p}_{i,j}=\sum _{k\notin C_j}d{p}_{i,k}(i>1)$$

在 $i=1$ 的初始情况下，$d{p}_{1,j}=1$。

注意到 $n\le {10}^{15}$，字符集的大小只有 $26$。考虑矩阵乘法。另 $F_n=\left[\begin{array}{c}d{p}_{n,1},d{p}_{n,2},\dots ,d{p}_{n,26}\end{array}\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;
}

小结

矩阵乘法加速递推可以说是我第一个学明白的数论算法了，虽然感觉在实战中的使用率不大，但是还是感觉是一个很有用的算法