Description
Alice想要得到一个长度为n的序列，序列中的数都是不超过m的正整数，而且这n个数的和是p的倍数。
Alice还希望，这n个数中，至少有一个数是质数。
Alice想知道，有多少个序列满足她的要求。
Input
一行三个数，n,m,p。
1<=n<=10^9,1<=m<=2×10^7,1<=p<=100
Output
一行一个数，满足Alice的要求的序列数量，答案对20170408取模。
Sample Input
3 5 3
Sample Output
33
题解Here!
据说正解是DP+矩阵快速幂。。。
然而我用了生成函数+快速幂，FFT都没用到。。。
首先这题显然补集转化，就是用全部方案减去不含任何质数的方案。
考虑m比较小，我们能大力把<=m的质数全都筛出来。
发现n很大，要么倍增要么快速幂。。。
发现p相当小，所以我们能在mod p的同余系下做啊。
一看到同余系下求方案数立刻想到卷积和生成函数。。。
假设我们有一个多项式f(x)，其中 xi 的系数为 a个数的序列之和 mod p为 i 的方案数(a为我们引入的变量)。
同时我们有另一个多项式g(x),其中 xi 的系数为b个数的序列之和 mod p为 i 的方案数(b为我们引入的变量)。
那么，我们如果让f(x)和g(x)做卷积的话，新的多项式 xi 的系数就是(a+b)个数的序列之和 mod p为 i 的方案数。
这就是生成函数了。
回到这个题，我们先初始化多项式f(x)，令 xi 的系数为 1个数的序列之和 mod p为 i 的方案数。
然后我们求出这个多项式的n次方，就是我们需要的答案了。
发现这道题的p很小，我们连FFT都不用，直接用一个多项式类暴力快速幂就行了。复杂度O(m+p2*log2n)，跑的飞快。。。
附代码：
#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#define MAXN 110
#define MAXM 2000010
#define MOD 20170408
using namespace std;
int n,m,p,k=0,prime[MAXM];
bool np[MAXM*10];
struct node{
long long num[MAXN];
node(){memset(num,0,sizeof(num));}
}one,two;
inline int read(){
int date=0,w=1;char c=0;
while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}
while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}
return date*w;
}
void make(){
np[0]=np[1]=true;
for(int i=2;i<=m;i++){
if(!np[i])prime[++k]=i;
for(int j=1;j<=k&&(long long)prime[j]*i<=m;j++){
np[prime[j]*i]=true;
if(i%prime[j]==0)break;
}
}
for(int i=1;i<=m;i++){
one.num[i%p]++;
if(np[i])two.num[i%p]++;
}
}
node operator *(const node &x,const node &y){
node ret;
for(int i=0;i<p;i++)
for(int j=0;j<p;j++)
ret.num[(i+j)%p]=(ret.num[(i+j)%p]%MOD+x.num[i]%MOD*y.num[j]%MOD)%MOD;
return ret;
}
node mexp(node base,int k){
node ans=base;
k--;
while(k){
if(k&1)ans=ans*base;
base=base*base;
k>>=1;
}
return ans;
}
int main(){
n=read();m=read();p=read();
make();
one=mexp(one,n);
two=mexp(two,n);
long long ans=(one.num[0]-two.num[0]+MOD)%MOD;
printf("%lld
",ans);
return 0;
}