「luogu3702」[SDOI2017]序列计数

考虑题目中的两个限制：

• 至少有一个数是质数：简单容斥，至少有一个数是质数的方案数=总方案数-没有一个数是质数的方案数。
• 所有数的和是p 的倍数：因为p的范围很小，可以在模p意义下考虑。

然后发现这是一个dp。

dp[i][j] 表示前 i 个数的和模 p 为 j 的方案数。

发现可以用矩阵快速幂优化。

#include<bits/stdc++.h>
using namespace std;
const int N=20000010,P=110,mod=20170408;
int n,m,p,prime[N],tot,cnt[P],all[P];
bool isp[N];
void pre(){
    isp[1]=1,cnt[1%p]++,all[1%p]++;
    for(int i=2;i<=m;i++){
        all[i%p]++;
        if(!isp[i]){prime[++tot]=i;}
        else cnt[i%p]++;
        for(int j=1;j<=tot&&1LL*prime[j]*i<=m;j++){
            isp[i*prime[j]]=1;
            if(!(i%prime[j])) break;
        }
    }
    return;
}
struct Mat{
    int n,m,num[P][P];
    Mat(){
        n=m=0;
        memset(num,0,sizeof(num));
        return;
    }
    void set1(){for(int i=1;i<=n;i++) num[i][i]=1;return;}
    Mat operator*(const Mat& k)const{
        Mat c;
        c.n=n,c.m=k.m;
        for(int i=0;i<n;i++)
            for(int j=0;j<k.m;j++)
                for(int ii=0;ii<m;ii++)
                    c.num[i][j]=((long long)c.num[i][j]+(long long)num[i][ii]*k.num[ii][j])%mod;
        return c;
    }
    void print(){
        for(int i=0;i<n;i++){
            for(int j=0;j<m;j++) cout<<num[i][j]<<" ";
            cout<<endl;
        }
        cout<<endl;
        return;
    }
};
void qpow(Mat x,int y,Mat& res){
    res.n=x.n,res.m=x.m;
    res.set1();
    while(y){
        if(y&1) res=res*x;
        x=x*x,y>>=1;
    }
    return;
}

int main(){
    scanf("%d%d%d",&n,&m,&p);
    pre();
    Mat a,b,c,d;
    a.n=a.m=b.n=b.m=p;
    c.n=d.n=p,c.m=d.m=1;
    for(int i=0;i<p;i++){
        for(int j=0;j<p;j++) a.num[i][j]=cnt[(i-j+p)%p],b.num[i][j]=all[(i-j+p)%p];
    }
    c.num[0][0]=d.num[0][0]=1;
    qpow(a,n-1,a);qpow(b,n-1,b);
    c=a*c,d=b*d;
    int ans=((d.num[0][0]-c.num[0][0])%mod+mod)%mod;
    printf("%d",ans);
    return 0;
}