矩阵加速递推浅谈
前言:
以下是一位矩阵蒟蒻学习矩阵加速的悲哀历史
(我才不会说是我不认真学)
前置芝士——矩阵快速幂:
矩阵加速,顾名思义,用到矩阵进行加速,而这个算法的加速核心就是矩阵快速幂 。(我只会背板子,逃~
在此对矩阵快速幂略作解释 :
快速幂大家应该都知道,矩阵快速幂就是把普通乘法换为矩阵乘法。
而矩阵乘法大家也应都有了解,大致方式如下:
\[\begin{bmatrix}a&c\\b&d\end{bmatrix}
\times
\begin{bmatrix}e&g\\f&h\end{bmatrix}=
\begin{bmatrix}ae\!+\!cf&ag\!+\!ch\\be\!+\!df&bg\!+\!dh\end{bmatrix}
\quad
\]
因此矩阵平方:
\[\begin{bmatrix}a_1&a_3\\a_2&a_4\end{bmatrix}^2=
\begin{bmatrix}a_1\!a_1\!+\!a_3\!a_2&a_1\!a_3\!+\!a_3\!a_4 \\a_2\!a_1\!+\!a_4\!a_2&a_2\!a_3\!+\!a_4\!a_4\end{bmatrix}
\]
各位应该也能发现一个矩阵乘法的性质:如果两个矩阵能相乘,则第一个矩阵的行数与第二个矩阵的列数必须相等。
所以我们接下来构造矩阵,一般构造成 \(n\!\times\!n\)的形式,以便计算。
切入正题——矩阵加速:
DP or 递推是我们在经常用到的做题方法,一般时间复杂度为\(O(n)\)或\(O(n^2)\),但遇到一些比较恶心的出题人,给你来一个\(n\le10^{16}\),你就可以和这个可爱的世界说再见了。而这个时候,便可能要用到矩阵加速。
例题1 :(P1962)斐波那契数列
这道题我们已经很熟了,但大家也应该清楚,随着数据量的增加,橙->绿->紫,紫题今天尚且不讲 (其实是我不会),我们来谈谈\(n<2^{63}\)。
本蒟蒻语文不好,没法子引入,所以也不拐弯抹角了,直接上矩阵。
众所周知,\(f[n]=f[n-1]+f[n-2]\),所以我们可以看成$ f[n+2]=f[n+1]* 1+f[n]* 1 $ 因此矩阵就构成了:
\[\begin{bmatrix}f[n]&f[n+1]\end{bmatrix}
\times
\begin{bmatrix}?&1\\?&1\end{bmatrix}=
\begin{bmatrix}?&f[n+2]\end{bmatrix}
\]
为了保证矩阵不断相乘的连续性,即因为第一个矩阵左边\(f\)的下标比右边小一,所以得出矩阵也应如此,所以:
\[\begin{bmatrix}f[n]&f[n+1]\end{bmatrix}\times
\begin{bmatrix}?&1\\?&1\end{bmatrix}=
\begin{bmatrix}f[n+1]&f[n+2]\end{bmatrix}
\]
接下来中间矩阵大家就应该知道长啥样了,所以整个柿子为:
\[\begin{bmatrix}f[n]&f[n+1]\end{bmatrix}\times
\begin{bmatrix}0&1\\1&1\end{bmatrix}=
\begin{bmatrix}f[n+1]&f[n+2]\end{bmatrix}
\]
又因为:
\[\begin{bmatrix}f[n+1]&f[n+2]\end{bmatrix}\times
\begin{bmatrix}0&1\\1&1\end{bmatrix}=
\begin{bmatrix}f[n+2]&f[n+3]\end{bmatrix}
\]
所以以此类推:
\[\begin{bmatrix}f[1]&f[2]\end{bmatrix}\times
\begin{bmatrix}0&1\\1&1\end{bmatrix}^{n-1}=
\begin{bmatrix}f[n]&f[n+1]\end{bmatrix}
\]
柿子推完了,接下来就是矩阵快速幂的事了,代码网上满天飞,我就不贴了。(就是我懒)
另外,
这题的柿子还可以写成如下形式:
\[\begin{bmatrix}f[1]\\f[2]\end{bmatrix}\times
\begin{bmatrix}0&1\\1&1\end{bmatrix}^{n-1}=
\begin{bmatrix}f[n]\\f[n+1]\end{bmatrix}
\]
虽然说相乘方式变了,但为什么就不能说第一个矩阵是第二个,第二个矩阵是第一个?大家可以选择喜欢的方式推柿子。
(本人推荐第二种,原因在当你做了一道duliu矩阵加速题时便会知道)
例题2 :(P1707)刷题比赛
本题比较复杂 (恶心出题人),一大堆常数一下让本蒟蒻不知所措,但按照套路列两个矩阵,在慢慢算中间矩阵,还是可以的,就是比较烦。。。
但如此傲娇的我,怎能在这种事情上停住脚步?
矩阵如下:
\[\begin{bmatrix}a[n]\\a[n+1]\\b[n]\\b[n+1]\\c[n]\\c[n+1]\\k^2\\k\\1\\w^{k}\\z^{k}\end{bmatrix}\times
\begin{bmatrix}0&1&0&0&0&0&0&0&0&0&0\\q&p&0&1&0&1&r&t&1&0&0\\0&0&0&1&0&0&0&0&0&0&0\\0&1&v&u&0&1&0&0&0&1&0\\0&0&0&0&0&1&0&0&0&0&0\\0&1&0&1&y&x&0&1&2&0&1\\0&0&0&0&0&0&1&2&1&0&0\\0&0&0&0&0&0&0&1&1&0&0\\0&0&0&0&0&0&0&0&1&0&0\\0&0&0&0&0&0&0&0&0&w&0\\0&0&0&0&0&0&0&0&0&0&z\end{bmatrix}=
\begin{bmatrix}a[n+1]\\a[n+2]\\b[n+1]\\b[n+2]\\c[n+1]\\c[n+2]\\k^2+2k+1\\k+1\\1\\w^{k+1}\\z^{k+1}\end{bmatrix}
\]
终于打完了,累呀~~~
但我相信大家以后遇到这种duliu题应该都能有思路了。
终于到贴代码的时候了:
#include<bits/stdc++.h>
#define re register
using namespace std;
long long p,q,r,t,u,v,w,x,y,z;
long long n,mod,f[12],mp[12][12];
inline long long Mul(long long a,long long b)
{
re long long ans=0;
for(;b;b>>=1)
{
if(b&1)ans=(ans+a)%mod;
a=(a<<1)%mod;
}
return ans;
}
inline void work()
{
re long long mid[12];
memset(mid,0,sizeof(mid));
for(re int i=1;i<=11;i++)
for(re int j=1;j<=11;j++)
mid[i]=(mid[i]+Mul(f[j],mp[i][j]))%mod;
for(re int i=1;i<=11;i++)f[i]=mid[i];
}
inline void Pow()
{
re long long mid[12][12];
memset(mid,0,sizeof(mid));
for(re int i=1;i<=11;i++)
for(re int j=1;j<=11;j++)
for(re int k=1;k<=11;k++)
mid[i][j]=(mid[i][j]+Mul(mp[i][k],mp[k][j]))%mod;
for(re int i=1;i<=11;i++)
for(re int j=1;j<=11;j++)
mp[i][j]=mid[i][j];
}
int main()
{
scanf("%lld%lld%lld%lld%lld%lld%lld%lld%lld%lld%lld%lld",&n,&mod,&p,&q,&r,&t,&u,&v,&w,&x,&y,&z);
f[1]=f[3]=f[5]=f[7]=f[8]=f[9]=1;
f[2]=f[4]=f[6]=3,f[10]=w,f[11]=z;
mp[1][2]=mp[2][4]=mp[2][6]=mp[2][9]=mp[3][4]=mp[4][2]=mp[4][6]=mp[4][10]=mp[5][6]=mp[6][2]=mp[6][4]=mp[6][8]=mp[6][11]=mp[7][7]=mp[7][9]=mp[8][8]=mp[8][9]=mp[9][9]=1;
mp[2][1]=q,mp[2][2]=p,mp[2][7]=r,mp[2][8]=t,mp[4][3]=v,mp[4][4]=u,mp[6][5]=y,mp[6][6]=x,mp[10][10]=w,mp[11][11]=z,mp[6][9]=mp[7][8]=2;
re long long res=n-1ll;
for(;res;Pow(),res>>=1)
if(res&1)work();
printf("nodgd %lld\nCiocio %lld\nNicole %lld\n",f[1],f[3],f[5]);
return 0;
}
/*
nodgd 74
Ciocio 80
Nicole 59
*/
提醒:本题需要用到龟速乘(会爆long long)~
代码如下:
const long long mod=998244353;
inline long long Mul(long long a,long long b)
{
re long long ans=0;
for(;b;b>>=1)
{
if(b&1)ans=(ans+a)%mod;
a=(a<<1)%mod;
}
return ans;
}