题解 P3702 【[SDOI2017]序列计数】

· · 题解

这里介绍一下矩阵快速幂的进阶技巧——双矩阵快速幂

矩阵快速幂优化dp转移是一个非常暴力的手段，通过矩阵快速幂我们可以将原来O(NM)的复杂度偏移至O(M^{3}logN)在M和N差距甚大的情况下真的是再好不过了

前置芝士：矩阵快速幂优化dp

这里为了避免萌新不懂矩阵快速幂所以先讲下优化原理和构造注意事项

(熟练的诸位可以自行跳过)

(在此之前先百度下矩阵和矩阵乘法的定义)

首先有些时候假设我们有一个二维的计数dp，而且可以滚动数组优化

那么我们的dp方程如果恰好长这样，其中Tr_{i,j}表示i状态是否可以转移到j的话。

Dp_{i,j}=\sum_{k=0}^{n}Dp_{i-1,k}Tr_{k,j}

我们会发现这个东西和矩阵乘法(假设a×b=c)的式子非常像

C_{i,j}=\sum_{k=0}^{n}A_{i,k}B_{k,j}

非常像，更准确的说，如果把dp数组的某一行看做方阵的第一行(毕竟没啥人会认真的写非方阵乘法)的话

那么我们可以列出这样的矩阵乘法式子Dp=Dp×Tr就可以完成一次转移

那么如果我们进行了n次转移的话，就相当于Dp=Dp×Tr^{n}

然后根据矩阵乘法的结合律，Tr^{n}可以快速幂计算

这个就是矩阵快速幂的原理

通过上述的推导不难发现，在一次转移中:如果状态i有k种方式转移到状态j的话，我们会发现Tr_{i,j}=k因为tr数组此时替代了一般dp转移过程的判定是不是合法转移的判断语句的作用

当然在计数dp中因为乘法原理的关系，有些时候转移需要乘上一定的系数，而这个作用也被Tr数组包了

刚才说了这么多，就是希望大家意识到，我们的转移矩阵是单向的转移关系，并且每次转移严格的由行状态转移到列状态，构造的时候请务必注意

本题题解

这道题如果是去掉必须有质数的限制的话是非常简单的

直接DP_{i,j}表示长度为i的序列，且序列各元素和%p为j的方案数

然后显然如果状态p1需要转移到p2的话，有(m/p+(p2-p1)是否大于m\%p)种方法此时直接矩阵快速幂硬莽过去就可以了,最后的答案当然是Dp_{n,0}了

但是问题是我们有质数的限制啊……，此时的情况就会变得非常辣手

观察到m<=2*1e7，这是什么？这是卡掉埃式筛放欧拉筛过的数据！(尽管这样的话O(NloglogN)和O(N)也不是特别好卡)

所以先上一个欧拉筛筛出质数再说

然后我们为了方便起见，设zhimo_{i}表示%p为i，且小于等于m的质数有多少个

此时我们发现，其实只需要在刚才的dp中加一维k，(k∈\{0,1\})表示有无质数就可以转移了

怎么转移呢？如果0状态要转移到0状态的话，那么必须放非质数，如果1状态需要转移到1状态的话，当然是随便放了，如果是0状态要转移到1状态的话是只可以放质数的，如果是1状态转移到0状态的话，当然是不可以转移的

所以如果(0，p1)这个状态需要转移到(0,p2)的话，转移系数当然是(%p=p2-p1的数-zhimo_{p2-p1})啦

同理转移(0,p1)到(1,p2)的话，转移系数就是zhimo_{p2-p1}了

转移(1,p1)到(1,p2)的话，转移系数和没有质数限制的时候一样

此时我们发现，我们可以开一个200×200的转移矩阵

前一半是Dp_{0,j}后一半是Dp_{1,j}然后按照上述的算法转移就可以了

对了由于是正整数，所以需要把0减掉

就酱……

上代码~

#include<cstdio>
#include<algorithm>
using namespace std;const int M=2*1e7+10;const int N=210;typedef long long ll;const ll mod=20170408;
int n;int m;int p;int zhi[M];int book[M];int ct;int t[N][N];int siz;
int mo[N];int zmo[N];int nzmo[N];int del[N][N];
struct martix//矩阵类
{
    ll mp[N][N];
    void operator *=(const martix& b)
    {
        for(int i=0;i<siz;i++){for(int j=0;j<siz;j++){t[i][j]=0;}}
        for(int i=0;i<siz;i++)
            for(int j=0;j<siz;j++)
                for(int k=0;k<siz;k++)
                {t[i][j]=(t[i][j]+mp[i][k]*b.mp[k][j])%mod;}
        for(int i=0;i<siz;i++){for(int j=0;j<siz;j++){mp[i][j]=t[i][j];}}
    }
}st,tr,res;
int main()
{
    scanf("%d%d%d",&n,&m,&p);siz=2*p;
    for(int i=2;i<=m;i++)//线性筛素数
    {
        if(!book[i]){zhi[++ct]=i;}
        for(int j=1;j<=ct&&zhi[j]*i<=m;j++)
        {book[zhi[j]*i]=true;if(i%zhi[j]==0){break;}}
    }
    for(int i=1;i<=ct;i++){zmo[zhi[i]%p]++;}//处理出来素数%p为i的方案数
    for(int i=0;i<=m%p;i++){mo[i]=m/p+1;}//然后处理%p为i的方案数
    for(int i=m%p+1;i<p;i++){mo[i]=m/p;}
    mo[0]--;//正整数，去掉0
    for(int i=0;i<p;i++){nzmo[i]=mo[i]-zmo[i];}//%p为i且非质数
    for(int p1=0;p1<p;p1++)
    {for(int p2=0;p2<p;p2++){del[p1][p2]=(p2+p-p1)%p;}}
    for(int p1=0;p1<p;p1++)//0->0
    {
        for(int p2=0;p2<p;p2++)
        {tr.mp[p1][p2]=nzmo[del[p1][p2]];}
    }
    for(int p1=0;p1<p;p1++)//1->1
    {
        for(int p2=0;p2<p;p2++)
        {tr.mp[p+p1][p+p2]=mo[del[p1][p2]];}
    }
    for(int p1=0;p1<p;p1++)//0->1
    {
        for(int p2=0;p2<p;p2++)
        {tr.mp[p1][p+p2]=zmo[del[p1][p2]];}
    }
    for(int i=0;i<siz;i++){res.mp[i][i]=1;}//矩阵快速幂不解释
    for(int i=n;i;i>>=1,tr*=tr){if(i&1){res*=tr;}}
    st.mp[0][0]=1;st*=res;printf("%lld",st.mp[0][p]);return 0;//拜拜程序~
}