题解 P4456 【[CQOI2018]交错序列】

非常传统的双矩阵快速幂题目，但是这道题的亮点是卡常数吧…… 在我AC这道题的时候，可能会出现读作TLE写做RE的情况发生，如果发生了这种情况(通常是40pts)，请开始卡常吧…… # 双矩阵快速幂 如果对矩阵快速幂一点都不了解的话，请自行百度自学吧……，这里默认你至少会一般的矩阵快速幂。 所谓双矩阵快速幂呢，就是在同一张矩阵中放置两个(甚至更多的)dp数组，这样可以实现处理两个不同dp数组的相互转移行为，或者出现了$dp_{i,j}$需要由$dp_{i-1,k},dp_{i-2,k'}$转移过来的尴尬情况，此时我们就可以通过一张矩阵上放置两个dp数组的方式来转移(代价是8倍的常数)，最简单的双矩阵快速幂就是求fibonacci数列第n项了，因为你的矩阵里同时保存了$f_{i},f_{i-1}$两个项 那么直接以这道题来讲解一下如何使用双矩阵快速幂以及如何卡常数好了 答案让我们求,有两个元不是十分美观，所以我们呢削去一个元 ## $x^{a}y^{b}$ ## $=(n-y)^ay^b$ 使用二项式定理暴力展开 ## $=\sum_{i=0}^{a}C_{a}^{i}(-1)^{a-i}n^{i}y^{a-i}y^b$ ## $=\sum_{i=0}^{a}C_{a}^{i}(-1)^{a-i}n^{i}y^{a+b-i}$ 然后我们现在要求的是 所有交错序列的上边式子的值，我们发现可以对于每一个交错序列来讲 ## $C_{a}^{i}(-1)^{a-i}n^{i}$ 都是常量，所以可以提出来，最后对所有式子的特征值统一用二项式定理 现在问题变成了，求所有交错序列的$y^{a+b-i}$的和，然后这样我们才可以仅在最后使用一次二项式定理 对于这个问题可以使用dp解决它，令$Dp_{i,j,k}(k∈(0,1))$表示决策到了第i位结尾为0/1，所有方案1的个数的j次方之和 那么转移方程有2个，分别是(0->0)或者(0->1)或者(1->0) 然后转移方程分别是这2个 转移到0的时候因为没有增加1的个数可以直接相加 ## $Dp_{i,j,0}=Dp_{i-1,j,0}+Dp_{i-1,j,1}$ 转移到1的时候增加了1的个数所以可以暴力二项式定理展开，由于不能放两个连续的1，所以只可以由0转移过来 ## $Dp_{i,j,1}=\sum_{k=0}^{j}C_{j}^{k}Dp_{i-1,k,0}$ 因此我们在构造转移矩阵的时候可以一张矩阵上放两个dp数组，前半部分放$Dp_{0}$，后半部分放$Dp_{1}$然后就可以转移啦~ 当然构造转移矩阵时务必注意两点，转移的方向一定是从行到列的，转移从来都是单向转移。 转移矩阵的话，长相大概是左上角是单位矩阵(0->0)，左下角也是单位矩阵(1->0)，右上角是一个旋转90度的杨辉三角(0->1)，右下角全部是0(1->1) 然后矩阵快速幂就好了~ 做完了？ 你兴奋的写完了代码，加上去一看，40pts，只比$2^n$暴力高10分，手测一组数据，2.0s超了时限一倍……，我不服啊…… ## 我要卡常数！ 你仔细一看数据范围，发现m<=1e8而不是1e9 这意味着你可以在矩阵快速幂的时候在c\[i]\[j]通过$O(a+b)$次加法求出值之后再去膜，一次矩阵快速幂仅做了$O((a+b)^2)$次的取膜运算，卡常完毕…… 你说，我不服啊……我要m开到1e9！ 于是你尝试把long long换成了unsigned long long很好，你的程序变成了1.8s 然后你把循环的顺序由i,j,k换成了i,k,j，这会使你连续访问一段内存，大幅度减少cache miss，很好，现在你的程序可以在1.2s完成任务了…… 然后你坐在电脑前卡了半年常数……却发现怎么也卡不进1s……正当你将要弃疗的时候，你把转移矩阵的n次幂打了出来！ 你惊奇的发现，这个转移矩阵的n次幂是由4个上三角矩阵构成的！ 这意味这整个矩阵有1/2的元素是0！，既然是0我们为什么要算它呢？ 所以我们每次矩阵快速幂的时候只计算四个上三角矩阵的乘法，因为那些一定是0的部分我们没必要算……于是我们直接把常数除了一个2，然后现在我们可以在0.6s内通过这道题了…… 上代码~ ```C #include<cstdio> #include<algorithm> #include<ctime> using namespace std;const int M=190;typedef long long ll;int stat;int end; int n;int m;ll mod;ll c[M][M];ll res;ll ans[M];ll siz;ll hsz;int a;int b;ll mi[M]; struct mar { ll mp[M][M]; mar(){for(int i=0;i<siz;i++)for(int j=0;j<siz;j++)mp[i][j]=0;} friend mar operator *(const mar& a,const mar& b)//这里只做了4个上三角矩阵乘法，常数/2 { mar c; for(int i=0;i<hsz;i++)//左上角 for(int k=0;k<siz;k++) for(int j=i;j<hsz;j++) (c.mp[i][j]+=a.mp[i][k]*b.mp[k][j])%=mod; for(int i=0;i<hsz;i++)//右上角 for(int k=0;k<siz;k++) for(int j=hsz+i;j<siz;j++) (c.mp[i][j]+=a.mp[i][k]*b.mp[k][j])%=mod; for(int i=hsz;i<siz;i++)//左下角 for(int k=0;k<siz;k++) for(int j=i-hsz;j<hsz;j++) (c.mp[i][j]+=a.mp[i][k]*b.mp[k][j])%=mod; for(int i=hsz;i<siz;i++)//左下角 for(int k=0;k<siz;k++) for(int j=i;j<siz;j++) (c.mp[i][j]+=a.mp[i][k]*b.mp[k][j])%=mod; return c; } }tr,ret,st; int main() { scanf("%d%d%d%lld",&n,&a,&b,&mod);hsz=a+b+1;siz=2*hsz;mi[0]=1;//因为放了两个dp数组 for(int i=0;i<hsz;i++)c[0][i]=c[i][i]=1;//先打表组合数 for(int i=0;i<hsz;i++) for(int j=1;j<i;j++) c[j][i]=(c[j-1][i-1]+c[j][i-1])%mod; for(int i=0;i<hsz;i++)tr.mp[i][i]=1;for(int i=0;i<hsz;i++)tr.mp[hsz+i][i]=1;//单位矩阵部分 for(int i=0;i<hsz;i++)for(int j=i;j<hsz;j++)tr.mp[i][hsz+j]=c[i][j];//杨辉三角部分 for(int i=0;i<siz;i++)ret.mp[i][i]=1; for(int p=n;p;p>>=1,tr=tr*tr)if(p&1)ret=ret*tr;//矩阵快速幂 for(int i=0;i<hsz;i++)ans[i]=(ret.mp[0][i]+ret.mp[0][hsz+i])%mod;//计算y^(a+b-i)和 for(int i=1;i<hsz;i++){mi[i]=mi[i-1]*n%mod;}//计算幂次 for(int i=0;i<=a;i++)//统一使用二项式定理 { ll del=c[i][a]*mi[i]%mod*ans[a+b-i]%mod; if((a-i)%2){res=(res+mod-del)%mod;}else {res=(res+del)%mod;} }printf("%lld",res);return 0;//拜拜程序~ } ```