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
正解：矩阵快速幂/$FFT$+快速幂+中国剩余定理。
昨天晚上考这题，发现是一道$FFT$优化$DP$的裸题。然而数组开小了，所以只有$80$分。。
正解似乎是矩阵快速幂，不过为什么$FFT$跑得快一些。。并且当$p$很大时矩阵快速幂就没用了。。
首先我们可以想到一个$O(n^{3})$的暴力$DP$。设$f[i][j]$表示前$i$个数，模$p$为$j$的方案数，那么$f[i][(j+k) \mod p]+=f[i-1][j]$，其中$k$为枚举选哪个数。
同正解一样，我们求出所有数的情况，然后减去没有质数的情况，最后得到的就是至少有一个质数的情况。
显然，这是一个卷积的形式，那么我们可以考虑用$FFT$来优化$DP$。
首先构造一个模$p$的多项式，$a[i]$表示模$p$为$i$的数有多少个。那么直接用$FFT$和快速幂算出$a^{n}$就行了。
然后线性筛求出所有质数，构造出除去所有质数的多项式$a$，再用一次$FFT$算出$a^{n}$。
第一个多项式的$a[0]$就是模$p$为$0$的情况，第二个多项式的$a[0]$就是没有任何质数且模$p$为$0$的情况，易知两个$a[0]$相减即为答案。
不过这个模数会炸精度，我们把这个模数拆成$8,1091,2311$，分别算出模这$3$个数的答案，最后用中国剩余定理合并就行了。
总复杂度$O(plogplogn)$，所以比矩阵快速幂的$O(p^{3}logn)$要快。不过我没加任何常数优化，所以还是很慢。。
//It is made by wfj_2048~
#include <algorithm>
#include <iostream>
#include <complex>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <cmath>
#include <queue>
#include <stack>
#include <map>
#include <set>
#define rhl (20170408)
#define NN (20000010)
#define pi acos(-1.0)
#define inf (1<<30)
#define il inline
#define RG register
#define ll long long
#define C complex <long double>
#define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout) using namespace std; C a[],b[],c[],ans[]; int prime[NN],rev[],r64[],N,n,m,p,lg,cnt;
ll res[],Ans[],ans1,ans2,aans;
bool vis[NN]; il int gi(){
RG int x=,q=; RG char ch=getchar();
while ((ch<'' || ch>'') && ch!='-') ch=getchar();
if (ch=='-') q=-,ch=getchar();
while (ch>='' && ch<='') x=x*+ch-,ch=getchar();
return q*x;
} il void sieve(){
vis[]=;
for (RG int i=;i<=m;++i){
if (!vis[i]) prime[++cnt]=i;
for (RG int j=,k;j<=cnt;++j){
k=i*prime[j]; if (k>m) break;
vis[k]=; if (i%prime[j]==) break;
}
}
return;
} il void FFT(C *a,RG int n,RG int f){
for (RG int i=;i<n;++i)
if (i<rev[i]) swap(a[i],a[rev[i]]);
for (RG int i=;i<n;i<<=){
C wn(cos(pi/i),sin(f*pi/i)),x,y;
for (RG int j=;j<n;j+=(i<<)){
C w(,);
for (RG int k=;k<i;++k,w*=wn){
x=a[j+k],y=w*a[j+k+i];
a[j+k]=x+y,a[j+k+i]=x-y;
}
}
}
return;
} il void mul(C *a,C *b,RG int pp){
for (RG int i=;i<N;++i) c[i]=b[i]; FFT(a,N,),FFT(c,N,);
for (RG int i=;i<N;++i) a[i]*=c[i]; FFT(a,N,-);
memset(res,,sizeof(res));
for (RG int i=;i<N;++i){
res[i%p]+=(ll)(a[i].real()/N+0.5);
res[i%p]%=pp,a[i]=;
}
for (RG int i=;i<p;++i) a[i]=res[i]; return;
} il void qpow(C *a,RG int b,RG int pp){
for (RG int i=;i<N;++i) ans[i]=a[i]; b--;
while (b){ if (b&) mul(ans,a,pp); mul(a,a,pp),b>>=; }
memset(res,,sizeof(res));
for (RG int i=;i<N;++i)
res[i%p]+=(ll)ans[i].real(),res[i%p]%=pp;
return;
} il void exgcd(RG ll a,RG ll b,RG ll &x,RG ll &y){
if (!b){ x=,y=; return; }
exgcd(b,a%b,y,x); y-=(a/b)*x; return;
} il void work(){
n=gi(),m=gi(),p=gi(); for (N=;N<=(p<<);N<<=) lg++;
for (RG int i=;i<N;++i) rev[i]=rev[i>>]>>|((i&)<<(lg-));
r64[]=,r64[]=,r64[]=,sieve();
for (RG int k=;k<=;++k){
memset(a,,sizeof(a)); for (RG int i=;i<=m;++i) a[i%p].real()++;
qpow(a,n,r64[k]); ans1=res[],memset(a,,sizeof(a));
for (RG int i=;i<=m;++i) if (vis[i]) a[i%p].real()++;
qpow(a,n,r64[k]); ans2=res[],Ans[k]=(ans1-ans2+r64[k])%r64[k];
RG ll u=rhl/r64[k],v=r64[k],x=,y=; exgcd(u,v,x,y);
aans+=u*x%rhl*Ans[k],aans%=rhl;
}
printf("%lld",aans); return;
} int main(){
File("count");
work();
return ;
}