[BZOJ4818][SDOI2017]序列计数（矩阵加速dp）

题目：

我是超链接

题解：

首先设计状态$f[i][j]$表示选i个数%p余数为j的方案数
据说还至少要有一个质数？那就再设$g[i][j]$表示选i个数%p余数为j且不含一个质数的方案数
那么最后的答案就是$f[i][0]-g[i][0]$
然后看看这个数据范围啊，一看就要想到用矩阵快速幂加速转移。
然后我们发现每次转移的时候都要有一个转移系数，即1~m之内能给这个余数做出多少贡献和个数，形成的矩阵就可以进行加速转移了

代码：

#include <cstdio>
#include <cstring>
#define LL long long
using namespace std;
const int mod=20170408;
const int N=20000005;
int n,m,p,pri[N],tot,cnt[105][2];bool ss[N];
struct mat
{
    LL sq[105][105];
    void clear(){memset(sq,0,sizeof(sq));}
}f,g,F,G;
mat cf(mat a,mat b)
{
    mat c;c.clear();
    for (int i=0;i<p;i++)
      for (int j=0;j<p;j++)
        for (int k=0;k<p;k++)
          c.sq[i][j]=(c.sq[i][j]+a.sq[i][k]*b.sq[k][j]%mod)%mod;
    return c;
}
mat ksm(mat a,int k)
{
    mat ans=a;k--;
    for (;k;k>>=1,a=cf(a,a))
      if (k&1) ans=cf(ans,a);
    return ans;
}
void init(int n)
{
    ss[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!ss[i]) pri[++tot]=i;
        for (int j=1;j<=tot&&pri[j]*i<=n;j++)
        {
            ss[pri[j]*i]=1;
            if (i%pri[j]==0) break;
        }
    }
    for (int i=1;i<=n;i++)
    {
        if (ss[i])//不是质数 
        {
            cnt[i%p][0]++;
            if (cnt[i%p][0]==mod) cnt[i%p][0]=0;
        }
        cnt[i%p][1]++;
        if (cnt[i%p][1]==mod) cnt[i%p][1]=0;
    }
}
int main()
{
    scanf("%d%d%d",&n,&m,&p);
    init(m);
    for (int i=0;i<p;i++)
      f.sq[0][i]=cnt[i][1],
      g.sq[0][i]=cnt[i][0];
    for (int i=0;i<p;i++)
      for (int j=0;j<p;j++)
        F.sq[i][(i+j)%p]=cnt[j][1],G.sq[i][(i+j)%p]=cnt[j][0];
    F=ksm(F,n-1);G=ksm(G,n-1);
    f=cf(f,F); g=cf(g,G);
    LL ansf=f.sq[0][0],ansg=g.sq[0][0];
    printf("%lld",(ansf-ansg+mod)%mod);
}