分析
考虑容斥,用总方案减去全是合数的方案数,
可以发现 \(n\) 很大,\(p\) 很小,直接用矩阵乘法转移即可
代码
#include <cstdio>
#include <cctype>
#include <algorithm>
#define rr register
using namespace std;
const int mod=20170408; bool v[mod];
struct maix{int p[100][100];}A0,A1,ANS0,ANS1;
int n,m,k,prime[1300011],Cnt,c0[100],c1[100];
inline signed iut(){
rr int ans=0; rr char c=getchar();
while (!isdigit(c)) c=getchar();
while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
return ans;
}
inline signed mo(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline void mul(maix &A,maix B){
rr maix C;
for (rr int i=0;i<k;++i)
for (rr int j=0;j<k;++j){
C.p[i][j]=0;
for (rr int o=0;o<k;++o)
C.p[i][j]=mo(C.p[i][j],1ll*A.p[i][o]*B.p[o][j]%mod);
}
for (rr int i=0;i<k;++i)
for (rr int j=0;j<k;++j)
A.p[i][j]=C.p[i][j];
}
signed main(){
scanf("%d%d%d",&n,&m,&k),v[1]=1;
for (rr int i=2;i<=m;++i){
if (!v[i]) prime[++Cnt]=i;
for (rr int j=1;j<=Cnt&&prime[j]<=m/i;++j){
v[i*prime[j]]=1;
if (i%prime[j]==0) break;
}
}
for (rr int i=1;i<=m;++i) ++c0[i%k];
for (rr int i=1;i<=m;++i) if (v[i]) ++c1[i%k];
for (rr int i=0;i<k;++i)
for (rr int j=0;j<k;++j)
A0.p[i][(i+j)%k]=c0[j],A1.p[i][(i+j)%k]=c1[j];
ANS0.p[0][0]=ANS1.p[0][0]=1;
for (;n;n>>=1,mul(A0,A0),mul(A1,A1))
if (n&1) mul(ANS0,A0),mul(ANS1,A1);
return !printf("%d",mo(ANS0.p[0][0],mod-ANS1.p[0][0]));
}