[SDOI2017]序列计数

Description:

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

Hint:

对\(20\%\)的数据，\(1\leq n,m\leq100\)
对\(50\%\)的数据，\(1\leq m \leq 100\)
对\(80\%\)的数据，\(1\leq m\leq 10^6\)
对\(100\%\)的数据，\(1\leq n \leq 10^9,1\leq m \leq 2\times 10^7,1\leq p\leq 100\)

Solution:

非正解A掉的....

考虑构造生成函数

设多项式 \(x^i\) 项的系数表示总和模p为i的序列的方案数

观察到p非常小,把多项式暴力乘起来,快速幂优化一下就过了??

至于最少有一个素数的约束,不太好算,就转化一下,求所有方案-不含素数的方案

#include <map>
#include <set>
#include <stack>
#include <cmath>
#include <queue>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#define ls p<<1 
#define rs p<<1|1
using namespace std;
typedef long long ll;
const int mxn=2e7+5,mod=20170408;
int n,m,p,tot,pr[8000000],vis[mxn];

inline int read() {
    char c=getchar(); int x=0,f=1;
    while(c>'9'||c<'0') {if(c=='-') f=-1;c=getchar();}
    while(c<='9'&&c>='0') {x=(x<<3)+(x<<1)+(c&15);c=getchar();}
    return x*f;
}
inline int chkmax(int &x,int y) {if(x<y) x=y;}
inline int chkmin(int &x,int y) {if(x>y) x=y;}

struct ed {
    int to,nxt;
}t[mxn<<1];

struct poly {
    int x[205];
    poly() {
        memset(x,0,sizeof(x));
    }
    friend poly operator * (poly a,poly b) {
        poly res;
        for(int i=0;i<p;++i) 
            for(int j=0;j<p;++j) 
                (res.x[(i+j)%p]+=1ll*a.x[i]*b.x[j]%mod)%=mod;
        return res;
    }
}U,C;

void sieve() { //m在线性范围内,可以直接筛
    vis[1]=1; 
    for(int i=2;i<=m;++i) {
        if(!vis[i]) pr[++tot]=i;
        for(int j=1;j<=tot&&i*pr[j]<=m;++j) {
            vis[pr[j]*i]=1;
            if(i%pr[j]==0) break ;
        }
    }
}

void init() 
{
    for(int i=1;i<=m;++i) {
        ++U.x[i%p];
        if(vis[i]) ++C.x[i%p]; //分别预处理
    }
}

poly qpow(poly bs,int b) 
{
    poly res=bs; --b;
    while(b) {
        if(b&1) res=res*bs;
        bs=bs*bs;
        b>>=1;
    }	
    return res;
}

int main()
{
    n=read(); m=read(); p=read(); 
    sieve(); init(); U=qpow(U,n); C=qpow(C,n);
    printf("%d",(U.x[0]-C.x[0]+mod)%mod);
    return 0;
}