斐波那契数列
定义 \(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]\)。
我们就会发现
也就是说，\(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\)
做法还是矩阵加速递推，定义 \(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\)
\(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\)
本题的目标是 \(\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\) 行
- \(A_{k+1,k+1}=c_1+1\)
- \(A_{k+1,1}=-c_k\)
- 对于剩下的 \(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\)
矩阵加速优化动态规划。
在本题中，我们令 \(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;
}
小结
矩阵乘法加速递推可以说是我第一个学明白的数论算法了，虽然感觉在实战中的使用率不大，但是还是感觉是一个很有用的算法
标签：begin,right,end,matrix,矩阵,times,线性,递推 From： https://www.cnblogs.com/zuoqingyuan/p/18367078