[SDOI2017]序列计数

题目描述：

雾。

题目分析：

f1[i][j]表示任意选i个1~m内的数，和模p为j的方案数，f2[i][j]表示选i个1~m内的不是质数的数，和模p为j的方案数，f1[n][0]-f2[n][0]就是答案，矩阵乘法加速DP就可以了，复杂度O(m+p^3logn)
。

题目链接：

Ac 代码：

#include<cstdio>
#include<cstring>
#define MN 20000000
#define MP 1270610
#define MK 100
#define MOD 20170408
struct mat
{
    int z[MK][MK];
    mat(){memset(z,0,sizeof(z));}
    mat operator*(const mat&b)
    {
        mat c;int i,j,k;
        for(i=0;i<MK;++i)for(k=0;k<MK;++k)if(z[i][k])for(j=0;j<MK;++j)
            c.z[i][j]=(c.z[i][j]+1LL*z[i][k]*b.z[k][j])%MOD;
        return c;
    }
}a,x;
bool u[MN+5];
int p[MP+5],pn,s[MK];
int main()
{
    int n,m,k,i,j,ans;
    for(u[1]=1,i=2;i<=MN;++i)
    {
        if(!u[i])p[++pn]=i;
        for(j=1;i*p[j]<=MN;++j){u[i*p[j]]=1;if(i%p[j]==0)break;}
    }
    scanf("%d%d%d",&n,&m,&k);
    for(i=1;i<=m;++i)++s[i%k];
    a.z[0][0]=1;
    for(i=0;i<k;++i)for(j=0;j<k;++j)x.z[i][(i+j)%k]=s[j]%MOD;
    for(i=n;i;i>>=1,x=x*x)if(i&1)a=a*x;
    ans=a.z[0][0];
    memset(s,0,sizeof(s));memset(a.z,0,sizeof(a.z));
    for(i=1;i<=m;++i)if(u[i])++s[i%k];
    a.z[0][0]=1;
    for(i=0;i<k;++i)for(j=0;j<k;++j)x.z[i][(i+j)%k]=s[j]%MOD;
    for(i=n;i;i>>=1,x=x*x)if(i&1)a=a*x;
    return printf("%d",(ans-a.z[0][0]+MOD)%MOD)*0;
}