d
p
[
i
]
[
j
]
[
0
/
1
]
dp[i][j][0/1]
dp[i][j][0/1]表示前
i
i
i个数的和模
p
p
p余数为
j
j
j,没有质数/有质数。
c
n
t
[
i
]
cnt[i]
cnt[i]表示
1
1
1~
m
m
m中模
p
p
p余数为
i
i
i的数的个数。
p
c
n
t
[
i
]
pcnt[i]
pcnt[i]表示
1
1
1~
m
m
m中模
p
p
p余数为
i
i
i的质数的个数。
下面的数组下标减法为模
p
p
p意义下的。
那么
d
p
[
i
]
[
j
]
[
0
]
dp[i][j][0]
dp[i][j][0]必须用非质数转移:
d
p
[
i
]
[
j
]
[
0
]
=
∑
t
=
0
p
−
1
d
p
[
i
−
1
]
[
t
]
[
0
]
∗
(
c
n
t
[
j
−
t
]
−
p
c
n
t
[
j
−
t
]
)
dp[i][j][0]=\sum_{t=0}^{p-1} dp[i-1][t][0]*(cnt[j-t]-pcnt[j-t])
dp[i][j][0]=t=0∑p−1dp[i−1][t][0]∗(cnt[j−t]−pcnt[j−t])
d
p
[
i
]
[
j
]
[
1
]
dp[i][j][1]
dp[i][j][1]由前面的两个状态转移:
d
p
[
i
]
[
j
]
[
1
]
=
∑
t
=
0
p
−
1
d
p
[
i
−
1
]
[
t
]
[
0
]
∗
p
c
n
t
[
j
−
t
]
+
d
p
[
i
−
1
]
[
t
]
[
1
]
∗
c
n
t
[
j
−
t
]
dp[i][j][1]=\sum_{t=0}^{p-1} dp[i-1][t][0]*pcnt[j-t]+dp[i-1][t][1]*cnt[j-t]
dp[i][j][1]=t=0∑p−1dp[i−1][t][0]∗pcnt[j−t]+dp[i−1][t][1]∗cnt[j−t]
构造矩阵转移即可(矩阵分块)。复杂度约为 O ( p 3 ∗ l o g n ) O(p^3*logn) O(p3∗logn),还要加上一个 O ( m ) O(m) O(m)的线筛。
#include<bits/stdc++.h>
#define re register
#define cs const
cs int M=2e7+10,P=205,mod=20170408;
int n,m,p,lp,cnt[P],pcnt[P];
inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int mul(int x,int y){return 1ll*x*y%mod;}
inline void Add(int &x,int y){x=(x+y>=mod)?x+y-mod:x+y;}
inline int pdec(int x,int y){return x-y<0?x-y+p:x-y;}
struct matrix{
int a[P][P];
matrix(int t=0){
memset(a,0,sizeof a);
for(int re i=0;i<P;++i) a[i][i]=t;
}
friend inline matrix operator*(cs matrix &a,cs matrix &b){
matrix c;
for(int re i=0;i<lp;++i)
for(int re k=0;k<lp;++k)
for(int re j=0;j<lp;++j)
Add(c.a[i][j],mul(a.a[i][k],b.a[k][j]));
return c;
}
friend inline matrix operator^(matrix a,int b){
matrix ret(1);
for(;b;b>>=1,a=a*a) if(b&1) ret=ret*a;
return ret;
}
}trans,A;
int mark[M],Pri[M],tot=0;
inline void linear_sieves(){
mark[0]=mark[1]=1;
for(int re i=2;i<=m;++i){
if(!mark[i]) Pri[++tot]=i,++pcnt[i%p];
for(int re j=1;j<=tot&&i*Pri[j]<=m;++j){
mark[i*Pri[j]]=1;
if(i%Pri[j]==0) break;
}
}
}
inline void init(){
for(int re i=0,k;i<p;++i){
for(int re j=0;j<p;++j){
k=pdec(i,j);
trans.a[i][j]=cnt[k]-pcnt[k];
trans.a[i][j+p]=pcnt[k];
trans.a[i+p][j+p]=cnt[k];
}
}A.a[0][0]=1;
}
int main(){
// freopen("2302.in","r",stdin);
scanf("%d%d%d",&n,&m,&p);
lp=p<<1,linear_sieves();
for(int re i=0;i<p;++i) cnt[i]=m/p;
for(int re i=1;i<=(m%p);++i) ++cnt[i];
init(),printf("%d\n",(A*(trans^n)).a[0][p]);
}