BZOJ-4818: [Sdoi2017]序列计数

Description

Alice想要得到一个长度为n的序列，序列中的数都是不超过m的正整数，而且这n个数的和是p的倍数。Alice还希望，这n个数中，至少有一个数是质数。Alice想知道，有多少个序列满足她的要求。

设dp[i]表示%p=i的方案数，矩乘优化求方案数，总方案数减去没有质数的方案数即可。

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
const ll mod=20170408;
ll ans;
int n,m,p;
int pri[1501000],cnt;
bool v[20001000];
void initialize()
{
    int i,j; v[1]=1;
    for(i=2;i<=m;++i)
    {
        if(!v[i]) pri[++cnt]=i;
        for(j=1;i*pri[j]<=m&&j<=cnt;++j)
        {
            v[i*pri[j]]=1;
            if(i%pri[j]==0) break;
        }
    }
}
struct matrix1
{
    ll a[110][110];
    matrix1(){memset(a,0,sizeof a);}
    matrix1 operator *(const matrix1 z)
    {
        matrix1 re; int i,j,k;
        for(i=0;i<p;++i)
            for(j=0;j<p;++j)
                for(k=0;k<p;++k)
                    re.a[i][j]=(re.a[i][j]+a[i][k]*z.a[k][j])%mod;
        return re;
    }
};
struct matrix2
{
    ll a[110];
    matrix2(){memset(a,0,sizeof a);}
    matrix2 operator *(const matrix1 z)
    {
        matrix2 re; int i,j;
        for(i=0;i<p;++i)
            for(j=0;j<p;++j)
                re.a[i]=(re.a[i]+a[j]*z.a[j][i])%mod;
        return re;
    }
};
int main()
{
    scanf("%d%d%d",&n,&m,&p);
    initialize();
    int i,j;
    matrix2 x,y; matrix1 I1,I2;
    for(i=1;i<=m;++i)
    {
        if(v[i]) x.a[i%p]++;
        y.a[i%p]++;
    }
    for(i=0;i<p;++i)
        for(j=0;j<p;++j)
        {
            I1.a[i][(i+j)%p]=(I1.a[i][(i+j)%p]+x.a[j])%mod;
            I2.a[i][(i+j)%p]=(I2.a[i][(i+j)%p]+y.a[j])%mod;
        }
    n--;
    while(n)
    {
        if(n&1) x=x*I1,y=y*I2;
        I1=I1*I1; I2=I2*I2; n>>=1;
    }
    printf("%lld",(y.a[0]-x.a[0]+mod)%mod);
    return 0;
}