矩阵乘法优化dp
矩阵乘法优化dp
前几天学姐说要给我们考矩乘dp和期望dp...其实我们都（也可能只有我）不会。
那只能现学了，然而学了一天突然通知不考了qwq
矩阵乘法
A矩阵为m*k，B矩阵为k*n，两矩阵相乘则得C矩阵为m*n;
for (int i=1;i<=M;++i) for (int j=1;j<=N;++j) for (int k=1;k<=K;++k) c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%mod;
时间复杂度为$O(N^3)$，数学一本通上提到了一种奇妙的矩阵乘法$Strassen：O(N^{2.81})$，但是非常非常复杂非常非常难写，于是就不写了（而且我也不会），最近又看到一种神奇的算法叫做$Coppersmith-Winograd$,$O(N^{2.36})$但是我觉得目前没必要学这个。
矩阵快速幂
其实和数字快速幂的思想是一样的。有一个小技巧:如果一个矩阵的主对角线全为1，其他部分为0，对于另一个矩阵来说就相当于单位元，相乘不变；
矩阵快速幂：https://www.luogu.org/problemnew/show/P3390
# include <cstdio> # include <iostream> using namespace std; const int P=1e9+7; struct s { long long a[101][101]; }A,ans; long long k; int n; s mul(s a,s b) { s c; for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) c.a[i][j]=0; for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) for (int k=1;k<=n;k++) c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%P)%P; return c; } void qui (long long k) { while (k) { if(k&1) ans=mul(ans,A); A=mul(A,A); k=k>>1; } return ; } int main() { scanf("%d%lld",&n,&k); for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) scanf("%lld",&A.a[i][j]); for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) if(i==j) ans.a[i][j]=1; qui(k); for (int i=1;i<=n;i++) { for (int j=1;j<=n;j++) printf("%lld ",ans.a[i][j]); printf("\n"); } return 0; }
矩乘优化数列递推
矩阵乘法可以用来加速线性的递推式(也许也可以加速别的递推吧)，思路非常棒;
首先构造一列数列，保存与递推$f[n]$有关的项，以斐波那契数列为例：
$f[n]=f[n-1]+f[n-2]$ 所以构造的矩阵$A$为
$\begin{pmatrix}f_{n-1} \\f_{n-2}\end{pmatrix}$
接下来考虑怎么转移，首先可以看出，需要转移到状态$B$
$\begin{pmatrix}f_{n} \\f_{n-1}\end{pmatrix}$
（因为我们知道这是矩阵乘法的题目），所以肯定要用矩阵来做，譬如设矩阵$X$，使得$AX=B$。
因为$f[n]=1*f[n-1]+1*f[n-2]$,$f[n-1]=1*f[n-1]+0*f[n-2]$
又有
$\begin{pmatrix}a&b\\c&d\end{pmatrix} \times \begin{pmatrix}e\\f\end{pmatrix}=\begin{pmatrix}a \times e+b \times f\\ c \times e+d \times f\end{pmatrix}$
于是就可以构造出转移矩阵$X$了：
$\begin{pmatrix} 1&1\\1&0\end{pmatrix}$
将原始的矩阵
$\begin{pmatrix} 1\\1\end{pmatrix}$
与X乘上(n-2)次后，[1][1]即为答案，因为矩阵乘法满足结合律，所以可以先将X矩阵进行n-2次的矩阵快速幂，时间复杂度$O(2^3*log(n))$
斐波那契数列：https://www.luogu.org/problemnew/show/P1962
# include <cstdio> # include <iostream> # include <cstring> using namespace std; const int wzx=1000000007; long long n; struct mat { long long a[3][3]; }a,ans; void init() { a.a[1][1]=a.a[1][2]=a.a[2][1]=1; } mat mul(mat a,mat b) { mat c; memset(c.a,0,sizeof(c.a)); for (int i=1;i<=2;i++) for (int j=1;j<=2;j++) for (int k=1;k<=2;k++) c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%wzx; return c; } void qui(mat a,long long x) { while (x) { if(x&1) ans=mul(ans,a); a=mul(a,a); x=x>>1; } } int main() { scanf("%lld",&n); memset(a.a,0,sizeof(a.a)); init(); ans.a[1][1]=ans.a[2][2]=1; qui(a,n-2); printf("%lld\n",(ans.a[1][1]+ans.a[1][2])%wzx); return 0; }
[模板]矩阵加速（数列）：https://www.luogu.org/problemnew/show/P1939
题意概述：$a[1]=a[2]=a[3]=1,a[x]=a[x-3]+a[x-1] (x>3)$，求$a[n]$
和刚刚的题目并没有什么差别qwq，构造矩阵：
$\begin{pmatrix} 1&0&1\\1&0&0\\0&1&0\end{pmatrix}$
# include <cstdio> # include <iostream> # include <cstring> using namespace std; const int wzx=1000000007; int T; long long n; struct mat { long long a[4][4]; }a,ans; void init() { a.a[1][1]=a.a[2][1]=a.a[1][3]=a.a[3][2]=1; ans.a[1][1]=ans.a[2][2]=ans.a[3][3]=1; } mat mul(mat a,mat b) { mat c; memset(c.a,0,sizeof(c.a)); for (int i=1;i<=3;i++) for (int j=1;j<=3;j++) for (int k=1;k<=3;k++) c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%wzx; return c; } void qui(mat a,long long x) { while (x) { if(x&1) ans=mul(ans,a); a=mul(a,a); x=x>>1; } } int main() { scanf("%d",&T); while (T--) { scanf("%lld",&n); if(n<4) { printf("1\n"); continue; } memset(a.a,0,sizeof(a.a)); memset(ans.a,0,sizeof(ans.a)); init(); qui(a,n-3); printf("%lld\n",(ans.a[1][1]+ans.a[1][2]+ans.a[1][3])%wzx); } return 0; }
Pell数列：from openjudge；
题意概述:$a[i]=2*a[i-1]+a[i-2]$，求$a[i]$。
因为并没有太大差别，所以没有写，只是推了一下矩阵：
$\begin{pmatrix} 2&1\\1&0\end{pmatrix}$
[NOI2012]随机数生成器：https://www.luogu.org/problemnew/show/P2044
题意概述：$X[n+1]=(aX[n]+c)%m$ 输出$X[n]%g$。
题目本身并没有什么难度：原始矩阵
$\begin{pmatrix} X_{0}\\ 1 \end{pmatrix}$
转移矩阵：
$\begin{pmatrix} a&c\\ 1&0 \end{pmatrix}$
注意:n是longlong范围的数，往快速幂里传参的时候别忘了！因为这个wa了3次；
还学了一个有趣的东西叫快速(慢速)乘，和快速幂很像，主要用于longlong相乘后取模，担心爆longlong的时候用：
LL mu(LL a,LL b) { LL res=0; while(b) { if(b&1) res=(res+a)%m; a=(a<<1)%m; b>>=1; } return res; }
其实是比直接乘慢的，如果非要说为什么叫快速乘，我猜是因为和快速幂长得像。
# include <cstdio> # include <iostream> # include <cstring> # define LL long long using namespace std; LL n,m,a,c,X,g; struct mat { LL a[3][3]; }A,ans; LL mu(LL a,LL b) { LL res=0; while(b) { if(b&1) res=(res+a)%m; a=(a<<1)%m; b>>=1; } return res; } mat mul(mat a,mat b) { mat anss; memset(anss.a,0,sizeof(anss.a)); for (int i=1;i<=2;++i) for (int j=1;j<=2;++j) for (int k=1;k<=2;++k) { anss.a[i][j]=(anss.a[i][j]+mu(a.a[i][k],b.a[k][j]))%m; if(anss.a[i][j]<0) anss.a[i][j]+=m; } return anss; } void qui(LL n) { while (n) { if(n&1) ans=mul(ans,A); A=mul(A,A); n=n>>1; } } int main() { scanf("%lld%lld%lld%lld%lld%lld",&m,&a,&c,&X,&n,&g); A.a[1][1]=a; A.a[1][2]=c; A.a[2][2]=1; ans.a[1][1]=ans.a[2][2]=1; qui(n); printf("%lld",((mu(ans.a[1][1],X)+ans.a[1][2])%m)%g); return 0; }
斐波那契公约数:https://www.luogu.org/problemnew/show/P1306
题意概述：给定$a,b$,求$gcd(F[a],F[b])$ $F[a]$为斐波那契数列的第a项。
首先这里有一个结论：$gcd(f[a],f[b])=f[gcd(a,b)]$。
然后就可以做了。
# include <cstdio> # include <iostream> # include <cstring> # define mod 100000000 using namespace std; int g_c_d(int a,int b) { return b?g_c_d(b,a%b):a; } struct mat { long long a[3][3]; }ans,a; int n,m,x; void init() { a.a[1][1]=a.a[1][2]=a.a[2][1]=1; ans.a[1][1]=ans.a[2][2]=1; } mat mul(mat a,mat b) { mat c; memset(c.a,0,sizeof(c.a)); for (int i=1;i<=2;++i) for (int j=1;j<=2;++j) for (int z=1;z<=2;++z) c.a[i][j]=(c.a[i][j]+a.a[i][z]*b.a[z][j])%mod; return c; } void qui(int x) { while (x) { if(x&1) ans=mul(ans,a); a=mul(a,a); x=x>>1; } } int main() { scanf("%d%d",&n,&m); x=g_c_d(n,m); if(x>2) { init(); qui(x-2); printf("%lld",(ans.a[1][1]+ans.a[1][2])%mod); return 0; } else printf("1"); return 0; }
emmm很明显我跳过了这道题真正的考点，于是我们现在可以来证明一下了。
假设$f[n]=a,f[n+1]=b$
所以$f[n+2]=a+b,f[n+3]=a+2b$
所以$f[m]=f[m-n+1]*a+f[m-n]*b$
$f[m]=f[m-n-1]*f[n]+f[m-n]*f[n+1]$
$gcd(f[n],f[m])=gcd(f[n],f[m-n-1]*f[n]+f[m-n]*f[n+1])$
很显然$f[n] | f[m-n-1]*f[n]$
$gcd(f[n],f[m-n]*f[n+1])$ (更相减损术)
现在把这个结论先放一放，去证明一个有趣的引理。
引理1：$gcd(f[n],f[n+1])=gcd(f[n],f[n]+f[n-1])$
$ =gcd(f[n-1],f[n]) $
$ ...... $
$ =gcd(f[1],f[2])=1$
于是就证明了一个结论，斐波那契数列相邻的两项互质。(让我们回到之前的结论),不，我们还有第二个引理。
引理2：当$gcd(a,b)=1$时，$gcd(a,b*c)=gcd(a,c)$
最近我发现关于gcd的题目，一个是递归证明，还有一个很棒的方法就是唯一分解定理。
$gcd(a,b)=p_{1}^{min(a_{1},b_{1})}*p_{2}^{min(a_{2},b_{2})}……$
因为$gcd(a,b)=1$
所以对于任意的$i$,$min(a_{i},b_{i})=0$
$gcd(a,b*c)=p_{1}^{min(a_{1},b_{1}+c_{1})}*p_{2}^{min(a_{2},b_{2}+c_{2})}……$
$gcd(a,c)=p_{1}^{min(a_{1},c_{1})}*p_{2}^{min(a_{2},c_{2})}……$
分情况讨论，如果$a_{i}=0$,那么两个$min$都等于0；
如果$b_{i}=0$,那么$b_{i}+c_{i}=c_{i}$
这样就证明了引理2。
终于可以回到原先的结论了。
$gcd(f[n],f[m])=gcd(f[n],f[m-n]*f[n+1])$
由引理1，$gcd(f[n],f[n+1])=1$,运用引理2，$gcd(f[n],f[m])=gcd(f[n],f[m-n])$
不停地减下去，就成了$gcd(f[n],f[m\%n])$
所以等于$gcd(f[0],f[gcd(n,m)])$
于是终于可以证明这个奇妙的定理啦。证明的过程还是很好玩的，只是两个引理之前不知道，在上面花了很长的时间。
[NHOI2011]数学作业：https://www.lydsy.com/JudgeOnline/problem.php?id=2326
题意概述：把1-n首尾相接连成一长串数后对m取模。举例：1~11=1234567891011;$n<=10^{18}$
bz上随便找道题就好神啊！怎么会有这么神的题目呢？
cnblog的公式真的难用啊，常常匹配不上又不能动态显示，决定使用印象笔记的Markdown编辑器。
首先把这个问题简化一下，假设所有的数字都一样长，怎么做呢？构造一个矩阵作为转移，可以做到$logN$的复杂度，见下：
事实上长度是不固定的，那怎么做呢？其实也很简单，分段进行矩阵乘法就可以啦。什么是分段乘啊？自己感觉一下吧。
1 # include <cstdio> 2 # include <iostream> 3 # include <cstring> 4 # define LL long long 5 6 LL n; 7 int m; 8 LL po[25],lans=0,zs,sum=0; 9 10 struct mat 11 { 12 LL a[3][3]; 13 }a,ans,c; 14 15 void init (int i) 16 { 17 ans.a[0][0]=ans.a[1][1]=ans.a[2][2]=1LL; 18 ans.a[0][1]=ans.a[0][2]=ans.a[1][0]=ans.a[1][2]=ans.a[2][0]=ans.a[2][1]=0LL; 19 a.a[0][0]=po[i]; 20 a.a[0][1]=a.a[0][2]=a.a[1][1]=a.a[1][2]=a.a[2][2]=1LL; 21 a.a[1][0]=a.a[2][0]=a.a[2][1]=0LL; 22 } 23 24 LL mu (LL a,LL b) 25 { 26 LL ans=0; 27 while (b) 28 { 29 if(b&1LL) ans=(ans+a)%m; 30 a=(a+a)%m; 31 b=b>>1LL; 32 } 33 return ans; 34 } 35 36 mat mul (mat a,mat b) 37 { 38 std::memset(c.a,0,sizeof(c.a)); 39 for (int i=0;i<=2;++i) 40 for (int j=0;j<=2;++j) 41 for (int k=0;k<=2;++k) 42 c.a[i][j]=(c.a[i][j]+mu(a.a[i][k],b.a[k][j]))%m; 43 return c; 44 } 45 46 void qui (mat a,long long n) 47 { 48 while (n) 49 { 50 if(n&1LL) ans=mul(ans,a); 51 a=mul(a,a); 52 n=n>>1LL; 53 } 54 } 55 56 int main() 57 { 58 scanf("%lld%d",&n,&m); 59 po[0]=1; 60 if(m==1) 61 { 62 printf("0"); 63 return 0; 64 } 65 for (int i=1;i<=18;++i) 66 po[i]=po[i-1]*10LL; 67 for (int i=1;i<=18;++i) 68 { 69 if(po[i]-1<=n) 70 zs=po[i]-po[i-1]; 71 else 72 zs=n-po[i-1]+1; 73 init(i); 74 qui(a,zs); 75 lans=(mu(lans,ans.a[0][0])+mu(ans.a[0][1],po[i-1]-1)+ans.a[0][2])%m; 76 if(po[i]-1>=n) break; 77 } 78 printf("%lld",lans); 79 return 0; 80 }
[NOI2013]矩阵游戏：https://www.lydsy.com/JudgeOnline/problem.php?id=3240
题意概述：很简约的三个式子，答案却很麻烦。求f[n][m]%1,000,000,007.$n,m<=10^{1000000}$
哇这个数据范围！第一反应：这题的n，m是不是和答案没关系啊...那是不可能的。
首先不要想n,m的数据范围，想一想这道题怎么做(模拟除外)。
如果只有一行很明显可以矩阵加速，但是两行之间的转移关系就不很简单了，经过一番思索发现行与行之间的矩阵虽然不一样，但是依然是乘法，所以可以乘在一起啊。比如说A矩阵是横行的转移矩阵，再乘一个B矩阵就可以转移到下一行，每行要乘(m-1)次A，一共要向下转移(n-1)行，最后再乘上(m-1)个A以及初始矩阵就是答案啦。
注意这个地方千万不要为了计算方便瞎化简，因为矩阵乘法不遵循交换律！平时做矩阵快速幂的时候不大容易注意到是因为每次乘的都是相同的矩阵，但是这里就不一样了。
这样就做完了吗？别忘了m,n的真实数据范围，这时候就要请出费马小定理了，因为模数是个质数，所以采用费马小定理减小m,n。然后WA。后来上洛谷交了一次发现如果a或c等于1需要特判一下。这样就愉快的AC了。先贴代码再分析特判的问题吧。
1 # include <cstdio> 2 # include <iostream> 3 # include <cstring> 4 # include <string> 5 # define wzx 1000000007 6 7 long long n=0,m=0,A,B,C,D; 8 char c1[1000010],c2[1000010]; 9 struct mat 10 { 11 long long a[2][2]; 12 }a,b,c; 13 int l1,l2; 14 15 long long mu (long long a,long long b) 16 { 17 long long ans=0; 18 while (b) 19 { 20 if(b&1LL) ans=(ans+a)%wzx; 21 a=(a+a)%wzx; 22 b=b>>1LL; 23 } 24 return ans; 25 } 26 27 mat mul (mat a,mat b) 28 { 29 mat c; 30 c.a[0][0]=c.a[1][0]=c.a[0][1]=c.a[1][1]=0LL; 31 for (int i=0;i<=1;++i) 32 for (int j=0;j<=1;++j) 33 for (int k=0;k<=1;++k) 34 c.a[i][j]=(c.a[i][j]+mu(a.a[i][k],b.a[k][j]))%wzx; 35 return c; 36 } 37 38 mat qui (mat a,long long x) 39 { 40 mat ans; 41 ans.a[0][0]=ans.a[1][1]=1LL; 42 ans.a[0][1]=ans.a[1][0]=0LL; 43 while (x) 44 { 45 if(x&1LL) ans=mul(ans,a); 46 a=mul(a,a); 47 x=x>>1LL; 48 } 49 return ans; 50 } 51 52 int main() 53 { 54 long long modd; 55 scanf("%s%s",c1+1,c2+1); 56 scanf("%lld%lld%lld%lld",&A,&B,&C,&D); 57 l1=strlen(c1+1); 58 l2=strlen(c2+1); 59 if(A==1) modd=wzx; else modd=wzx-1; 60 for (int i=1;i<=l1;++i) 61 n=(n*10LL+c1[i]-'0')%modd; 62 n=(n-1+modd)%modd; 63 if(C==1) modd=wzx; else modd=wzx-1; 64 for (int i=1;i<=l2;++i) 65 m=(m*10LL+c2[i]-'0')%modd; 66 m=(m-1+modd)%modd; 67 a.a[0][0]=A; 68 a.a[0][1]=B; 69 a.a[1][0]=0LL; 70 a.a[1][1]=1LL; 71 a=qui(a,m); 72 c=a; 73 b.a[0][0]=C; 74 b.a[0][1]=D; 75 b.a[1][0]=0LL; 76 b.a[1][1]=1LL; 77 b=mul(a,b); 78 b=qui(b,n); 79 b=mul(b,c); 80 printf("%lld",(b.a[0][0]+b.a[0][1])%wzx); 81 }
关于特判的解释：咕咕咕，别急，在路上了。
佳佳的$Fibonacci$:https://loj.ac/problem/10222
题意概述:$f_i$表示斐波那契第$i$项,$T_i=\sum_{i=1}^n f_i \times i$,求$T_n$,$n<=2^{31}-1$
一本通上给出的解法比较数学化,类似于等比数列求和的做法,只需要一个$4 \times 4$的矩阵，我的做法比较暴力，需要$5 \times 5$矩阵,其实也差不了多少嘛.
$X \times \begin{pmatrix}T_{n-1} \\ (n-1) \times F_{n-1} \\ (n-2) \times F_{n-2} \\ F_{n-1}\\F_{n-2} \end{pmatrix} =\begin{pmatrix}T_{n} \\ n\times F_{n} \\ (n-1) \times F_{n-1} \\ F_{n}\\F_{n-1} \end{pmatrix}$
$X=\begin{pmatrix} 1& 1 &1 & 1 &2\\0 & 1 & 1 &1 &2\\ 0&1&0&0&0\\ 0&0&0&1&1\\ 0&0&0&1&0\end{pmatrix}$
这公式竟然没挂!震惊.
1 # include <cstdio> 2 # include <iostream> 3 # include <cstring> 4 # define R register int 5 # define ll long long 6 7 using namespace std; 8 9 int n,m; 10 struct mat 11 { 12 int n,m; 13 ll a[5][5]; 14 void init() 15 { 16 memset(a,0,sizeof(a)); 17 for (R i=0;i<5;++i) 18 a[i][i]=1; 19 } 20 void cler() { memset(a,0,sizeof(a)); } 21 mat operator * (const mat &a) const 22 { 23 mat c; 24 c.cler(); 25 for (R i=0;i<this->n;++i) 26 for (R j=0;j<this->m;++j) 27 for (R k=0;k<a.m;++k) 28 c.a[i][j]=(c.a[i][j]+this->a[i][k]*a.a[k][j]%m)%m; 29 c.n=n; 30 c.m=a.m; 31 return c; 32 } 33 }; 34 35 mat qui (mat a,int p) 36 { 37 mat ans; 38 ans.init(); 39 ans.n=ans.m=5; 40 while(p) 41 { 42 if(p&1) ans=ans*a; 43 a=a*a; 44 p>>=1; 45 } 46 return ans; 47 } 48 49 int main() 50 { 51 freopen("data.in","r",stdin); 52 53 scanf("%d%d",&n,&m); 54 if(n==1) 55 { 56 printf("%d",1%m); 57 return 0; 58 } 59 mat a,ans; 60 a.cler(),ans.cler(); 61 a.n=5,a.m=1; 62 a.a[0][0]=a.a[1][0]=a.a[3][0]=1; 63 ans.n=5,ans.m=5; 64 ans.a[0][0]=1,ans.a[0][1]=1,ans.a[0][2]=1,ans.a[0][3]=1,ans.a[0][4]=2; 65 ans.a[1][0]=0,ans.a[1][1]=1,ans.a[1][2]=1,ans.a[1][3]=1,ans.a[1][4]=2; 66 ans.a[2][0]=0,ans.a[2][1]=1,ans.a[2][2]=0,ans.a[2][3]=0,ans.a[2][4]=0; 67 ans.a[3][0]=0,ans.a[3][1]=0,ans.a[3][2]=0,ans.a[3][3]=1,ans.a[3][4]=1; 68 ans.a[4][0]=0,ans.a[4][1]=0,ans.a[4][2]=0,ans.a[4][3]=0,ans.a[4][4]=1; 69 ans=qui(ans,n-1); 70 ans=ans*a; 71 printf("%lld",ans.a[0][0]); 72 return 0; 73 }
Fib数列2：https://www.lydsy.com/JudgeOnline/problem.php?id=5118
题意概述：求$f[2^n]%1125899839733759,n<=10^15$。
数据范围极大，但是$2^n$似乎给了我一些启发。暴力分解质因数发现模数是个质数。“质数，幂，模”，几乎可以断定是费马小定理。但是$f[i]\%p=f[i\%p]$这个式子显然不成立，所以考虑一下别的方法。矩阵加速递推，这样就是矩阵的$2^n$次方了，感觉可做。平常写的矩阵加速通常是计算矩阵的$n-2$次方，但是对于这道题这样会很麻烦，所以对矩阵进行一些微调，最后使得正好乘$2^n$次。对于指数运用费马小定理降幂即可.
1 # include <cstdio> 2 # include <iostream> 3 # include <cstring> 4 # define R register int 5 # define LL long long 6 # define wzx 1125899839733759 7 8 int T; 9 LL n; 10 struct mat 11 { 12 LL a[2][2]; 13 }a,ans; 14 15 void init() 16 { 17 a.a[0][1]=a.a[1][1]=a.a[1][0]=1LL; 18 a.a[0][0]=0LL; 19 ans.a[0][0]=ans.a[1][1]=1LL; 20 ans.a[1][0]=ans.a[0][1]=0LL; 21 } 22 23 LL mu (LL a,LL b,LL p) 24 { 25 LL ans=0; 26 while (b) 27 { 28 if(b&1LL) ans=(ans+a)%p; 29 a=(a+a)%p; 30 b=b>>1LL; 31 } 32 return ans; 33 } 34 35 mat mul (mat a,mat b) 36 { 37 mat ans; 38 memset(ans.a,0,sizeof(ans.a)); 39 for (R i=0;i<=1;++i) 40 for (R j=0;j<=1;++j) 41 for (R k=0;k<=1;++k) 42 ans.a[i][j]=(ans.a[i][j]+mu(a.a[i][k],b.a[k][j],wzx))%wzx; 43 return ans; 44 } 45 46 void qui (LL x) 47 { 48 while (x) 49 { 50 if(x&1LL) ans=mul(ans,a); 51 a=mul(a,a); 52 x=x>>1LL; 53 } 54 } 55 56 LL ksm (LL a,LL n) 57 { 58 LL ans=1; 59 while (n) 60 { 61 if(n&1LL) ans=mu(ans,a,wzx-1); 62 a=mu(a,a,wzx-1); 63 n=n>>1LL; 64 } 65 return ans; 66 } 67 68 int main() 69 { 70 scanf("%d",&T); 71 while (T--) 72 { 73 scanf("%lld",&n); 74 init(); 75 n=ksm(2,n); 76 qui(n); 77 printf("%lld\n",ans.a[0][1]); 78 } 79 return 0; 80 }
[SNOI2017]礼物:https://www.lydsy.com/JudgeOnline/problem.php?id=5015
题意概述:设$f_1=1,s_n=\sum_{i=1}^nf_i,f_n=s_{n-1}+n^k$,求$f_n\%mod$,$n<=10^{18},k<=10$
$n$这么大的范围显然是矩阵乘法了,问题是怎样维护$n^k$.因为$k$很小,可以用二项式定理递推实现:$(x+1)^k=\sum_{i=0}^k \binom{k}{i}x^i$
把$n^{0...k}$全部压到状态里维护一下即可.
$\begin {pmatrix}s_{n} &n^k & n^{k-1} & ...& 1 \end {pmatrix} \times \begin {pmatrix} 1&C_k^k &C_{k}^{k-1} &...\\0&0 &C_{k-1}^{k-1}& ... \\ 0 & 0&0 \end {pmatrix}=\begin {pmatrix}s_{n+1} & (n+1)^k & (n+1)^{k-1} & ...& 1 \end {pmatrix}$
1 # include <cstdio> 2 # include <iostream> 3 # include <cstring> 4 # include <string> 5 # define ll long long 6 # define mod 1000000007 7 # define R register int 8 9 using namespace std; 10 11 const int maxn=12; 12 int K,c[maxn][maxn]; 13 ll n; 14 struct mat 15 { 16 int a[maxn][maxn]; 17 }ans,a; 18 19 mat mul (mat a,mat b) 20 { 21 mat c; 22 memset(c.a,0,sizeof(c.a)); 23 for (R k=0;k<=K+1;++k) 24 for (R i=0;i<=K+1;++i) 25 for (R j=0;j<=K+1;++j) 26 c.a[i][j]=(c.a[i][j]+1LL*a.a[i][k]*b.a[k][j])%mod; 27 return c; 28 } 29 30 void qui (mat a,ll x) 31 { 32 while(x) 33 { 34 if(x&1) ans=mul(ans,a); 35 a=mul(a,a); 36 x>>=1; 37 } 38 } 39 40 void init() 41 { 42 memset(ans.a,0,sizeof(ans.a)); 43 memset(a.a,0,sizeof(a.a)); 44 for (R i=0;i<=K;++i) 45 ans.a[i][i]=1; 46 a.a[0][0]=2; 47 for (R j=1;j<=K+1;++j) 48 a.a[0][j]=c[K][K-j+1]; 49 for (R i=1;i<=K+1;++i) 50 for (R j=i;j<=K+1;++j) 51 a.a[i][j]=c[K-i+1][K+1-j]; 52 } 53 54 ll getans (ll n) 55 { 56 init(); 57 qui(a,n-1); 58 ll fans=ans.a[0][0]; 59 for (R i=1;i<=K;++i) 60 fans=(fans+ans.a[0][i])%mod; 61 return fans; 62 } 63 64 int main() 65 { 66 scanf("%lld%d",&n,&K); 67 c[0][0]=1; 68 for (R i=1;i<=K;++i) 69 { 70 c[i][0]=1; 71 for (R j=1;j<=K;++j) 72 c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod; 73 } 74 ll fans=getans(n); 75 printf("%lld",fans); 76 return 0; 77 } 78 //1000000000000000000 10
矩阵优化图论
之前刚学矩乘的时候就觉得不止能用来做数列优化，今天果然见到了奇妙的应用。
图 -> 存图 -> 邻接矩阵...有没有想到矩阵乘法呢qwq
在一个有向图中，从A走到B，k步的方案数？将存图的邻接矩阵进行矩阵乘法k次，g[a][b]即为答案。为什么呢?因为矩阵乘法的本质啊...$c[i][j]+=a[i][k]*b[k][j]$ 很像floyd，即为用k中转一下从i到达j的方案数。
[HNOI2002]公交车路线：https://www.luogu.org/problemnew/show/P2233
题意概述：在一个环形路上乘公交车，求走n步到达4号点的路线数，注意：路上只要到达过四号点就不会再继续前进了，也就是说，只有在第n步才能去4号点。
这个题看起来很像是矩阵图论的模板，但是不能路过4号这个条件看起来有一点奇怪，不过回到矩阵优化图论的本质上去思考，$k$点其实就是用来中转的，而4号点不能用来中转，所以在$k=4$时$continue$一下就可以啦。
# include <cstdio> # include <iostream> # include <cstring> # define R register int # define mod 1000 using namespace std; int n; struct mat { int a[10][10]; }a,ans,c; void init() { for (R i=1;i<=8;++i) ans.a[i][i]=1; for (R i=1;i<=7;++i) a.a[i+1][i]=a.a[i][i+1]=1; a.a[1][8]=a.a[8][1]=1; } mat mul(mat a,mat b) { for (R i=1;i<=8;++i) for (R j=1;j<=8;++j) c.a[i][j]=0; for (R i=1;i<=8;++i) for (R j=1;j<=8;++j) for (R k=1;k<=8;++k) { if(k==5) continue; c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%mod; } return c; } void qui(int x) { while (x) { if(x&1) ans=mul(ans,a); a=mul(a,a); x=x>>1; } } int main() { cin>>n; init(); qui(n); printf("%d",ans.a[1][5]%mod); return 0; }
[SCOI2009]迷路:https://www.lydsy.com/JudgeOnline/problem.php?id=1297
题意概述：一个有向图中，边权都在$[1,9]$内,求$T$时间内恰好走到$n$的方案数.$n<=10,T<=10^9$
看看这几个题的年份,我觉得矩阵乘法可能真的过气了...
如果边权都是$1$,这就是一个矩乘快速幂的板子题了,可是这里出现了一个边权,难道是倍增$floyd$吗?并不是,因为要求的不是最短路.
其实这是一道网络流的题目.我的意思是它考察的是建图的能力.
发现边权非常小,那么可以将一个点拆成$9$ 个,表示从这个点出发已经走了$k$的长度但是还没走到下一个点的方案数,那么每次读入长度为$k$的边,就将起点的第$k$个虚点和终点连起来就好了.
1 # include <cstdio> 2 # include <iostream> 3 # include <cstring> 4 # include <string> 5 # define R register int 6 7 using namespace std; 8 9 const int maxn=91; 10 int id[11][10],cnt; 11 int n,t,mod=2009,len; 12 char s[maxn]; 13 14 struct mat 15 { 16 int n,m; 17 int a[maxn][maxn]; 18 void cler() 19 { 20 memset(a,0,sizeof(a)); 21 n=m=0; 22 } 23 void init() 24 { 25 cler(); 26 for (R i=0;i<=90;++i) 27 a[i][i]=1; 28 } 29 mat operator * (const mat &a) 30 { 31 mat c; 32 c.cler(); 33 c.n=this->n; 34 c.m=a.m; 35 for (R i=0;i<this->n;++i) 36 for (R j=0;j<a.m;++j) 37 for (R k=0;k<this->m;++k) 38 c.a[i][j]=(c.a[i][j]+this->a[i][k]*a.a[k][j]%mod)%mod; 39 return c; 40 } 41 void print() 42 { 43 for (R i=0;i<n;++i) 44 { 45 for (R j=0;j<m;++j) 46 printf("%3d",a[i][j]); 47 if(i!=n-1) printf("\n"); 48 } 49 } 50 }; 51 52 mat qui (mat a,int b) 53 { 54 mat c; 55 c.init(),c.n=a.n,c.m=a.m; 56 while(b) 57 { 58 if(b&1) c=c*a; 59 a=a*a; 60 b>>=1; 61 } 62 return c; 63 } 64 65 mat a; 66 67 int main() 68 { 69 scanf("%d%d",&n,&t); 70 for (R i=0;i<n;++i) 71 for (R j=1;j<=9;++j) 72 id[i][j]=++cnt; 73 a.n=a.m=cnt+1; 74 for (R i=0;i<n;++i) 75 for (R j=1;j<9;++j) 76 a.a[ id[i][j] ][ id[i][j+1] ]=1; 77 for (R i=0;i<n;++i) 78 { 79 scanf("%s",s); 80 for (R j=0;j<n;++j) 81 { 82 len=s[j]-'0'; 83 if(len) a.a[ id[i][len] ][ id[j][1] ]=1; 84 } 85 } 86 a=qui(a,t); 87 printf("%d",a.a[ id[0][1] ][ id[n-1][1] ]); 88 return 0; 89 }
---shzr
