这里提供一种复杂度 $O(p^2\log n)$ 的做法。
首先设 $f[i][j]$ 表示填了前 $i$ 个数,总和 $\bmod p=j$ 的方案数
然后有 $f[i][k]=f[i-1][i]*cnt[j]$,其中 $i+j \equiv k(\bmod p)$,$cnt[j]$ 表示在 $1$ 到 $m$ 中 $\equiv j(\bmod p)$ 的数量
然后考虑还有一层限制就是必须得有一个质数,那么正难则反我们可以容斥,即用:总方案数-没有质数的方案数=ans;
质数筛然后更改一下 $cnt[]$ 即可。
然后注意到这是一个线性递推,可以写成矩阵的方式。
然后就有经典的 $O(p^3 \log n)$ 的做法,不再赘述。
优化
我们观察到这个转移矩阵,他是个循环矩阵。具体来说:$ A_{i,j} = A_{0,(j-i) \bmod n}$
可以证明两个循环矩阵乘起来还是一个循环矩阵
$$ C_{i,j}=\sum\limits_{k=0}^{n}A_{i,k} \times B_{k,j}$$
对于 $C$ 的第一行有:
$$C_{0,j}=\sum\limits_{k=0}^{n} A_{0,k} \times B_{k,j}$$
$$C_{0,j}=\sum\limits_{k=0}^{n} A_{0,k} \times B_{0,(j-k)\bmod n}$$
对于剩余的部分:
$$C_{i,j}=\sum\limits_{k=0}^{n}A_{0,(k-i)\bmod n} \times B_{0,(j-k) \bmod n}$$
$$C_{i,j}=C_{0,(j-i) \bmod n}$$
现在我们知道其实只有第一行才是有效的,那么我们就可以 $O(n^2)$ 矩阵乘法然后再 $O(n^2)$ 填数字。
具体的有简化形式:
$$C_{0,k}=\sum\limits_{i+j \equiv k (\bmod n)}A_{0,i} \times B_{0,j}$$
对于这道题,转移矩阵是个循环矩阵,因此我们做 $O(n^2)$ 的矩阵乘法下的快速幂,然后再填满矩阵,获得最终的转移矩阵。
最后再乘上我们的初始向量即可。
代码:
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=1e2+5,M=2e7+10,mod=20170408;
int n,m,p;
int prs[M],idx;
bool st[M];
int cnt[2][N];
struct Matrix
{
int n;
int a[N][N];
inline int* operator[](const int i)
{
return a[i];
}
inline void resize(int _n)
{
n=_n;
memset(a,0,sizeof a);
}
inline void make_I()
{
for(int i=0;i<n;i++) a[i][i]=1;
}
inline Matrix operator*(Matrix B)const
{
Matrix res;
res.resize(n);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(int k=0;k<n;k++)
(res[i][j]+=a[i][k]*B[k][j]%mod)%=mod;
return res;
}
}F,B[2];
Matrix mul(Matrix A,Matrix B)
{
Matrix res;res.resize(A.n);
for(int i=0;i<A.n;i++)
for(int j=0;j<A.n;j++)
(res[0][(i+j)%p]+=A[0][i]*B[0][j]%mod)%=mod;
return res;
}
Matrix qmi(Matrix A,int k)
{
Matrix res;
res.resize(A.n),res[0][0]=1;
while(k)
{
if(k&1) res=mul(res,A);
A=mul(A,A);
k>>=1;
}
return res;
}
void get_prs(int n)
{
st[1]=true;
for(int i=2;i<=n;i++)
{
if(!st[i]) prs[++idx]=i;
for(int j=1;prs[j]<=n/i;j++)
{
st[prs[j]*i]=true;
if(i%prs[j]==0) break;
}
}
}
void get_cnt()
{
for(int i=1;i<=m;i++)
{
cnt[0][i%p]++;
if(st[i]) cnt[1][i%p]++;
}
}
void init_M()
{
for(int t=0;t<2;t++) B[t].resize(p);
for(int t=0;t<2;t++)
for(int i=0;i<p;i++)
B[t][0][i]=cnt[t][i];
for(int t=0;t<2;t++) B[t]=qmi(B[t],n);
for(int t=0;t<2;t++)
for(int i=1;i<p;i++)
for(int j=0;j<p;j++)
B[t][i][j]=B[t][0][(j-i+p)%p];
F.resize(p);
F[0][0]=1;
}
signed main()
{
ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
cin>>n>>m>>p;
get_prs(m);
get_cnt();
init_M();
Matrix res[2];
for(int t=0;t<2;t++) res[t]=F*B[t];
int ans=(res[0].a[0][0]-res[1].a[0][0]+mod)%mod;
cout<<ans<<"\n";
return 0;
}