表示个数,模为,有无质数的方案
倍增就完了啊
为什么要的矩乘啊
#include<bits/stdc++.h>
using namespace std;
const int RLEN=(1<<20)|5;
inline char gc(){
static char ibuf[RLEN],*ib,*ob;
(ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
return (ib==ob)?EOF:*ib++;
}
#define gc getchar
inline int read(){
char ch=gc();
int res=0,f=1;
while(!isdigit(ch))f^=ch=='-',ch=gc();
while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
return f?res:-res;
}
#define pb push_back
#define re register
#define fi first
#define se second
#define pii pair<int,int>
#define cs const
#define bg begin
#define ll long long
#define poly vector<int>
template<class tp>inline void chemx(tp &a,tp b){a<b?a=b:0;}
template<class tp>inline void chemn(tp &a,tp b){a>b?a=b:0;}
cs int mod=20170408;
inline int add(int a,int b){a+=b-mod;return a+(a>>31&mod);}
inline void Add(int &a,int b){a+=b-mod;a+=a>>31&mod;}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline void Dec(int &a,int b){a-=b;a+=a>>31&mod;}
inline int mul(int a,int b){static ll r;r=1ll*a*b;return r>=mod?r%mod:r;}
inline void Mul(int &a,int b){static ll r;r=1ll*a*b;a=r>=mod?r%mod:r;}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,Mul(a,a))(b&1)&&(Mul(res,a),1);return res;}
inline int Inv(int x){return ksm(x,mod-2);}
cs int N=20000005;
char xxx;
bitset<N> vis;
int pr[2000005],tot;
int n,m,p;
struct node{
int a0,a1;
node(int _a0=0,int _a1=0):a0(_a0),a1(_a1){}
inline void operator +=(cs node &x){
Add(a0,x.a0),Add(a1,x.a1);
}
friend inline node operator *(cs node &a,cs node &b){
return node(mul(a.a0,b.a0),(1ll*a.a0*b.a1+1ll*b.a0*a.a1+1ll*a.a1*b.a1)%mod);
}
}f[105],res[105],g[105];
inline void Mult(node *a,node *b,node *tmp){
for(int i=0;i<p;i++)g[i]=node();
for(int i=0;i<p;i++)
for(int j=0;j<p;j++)
g[(i+j)>=p?(i+j-p):(i+j)]+=a[i]*b[j];
for(int i=0;i<p;i++)tmp[i]=g[i];
}
char yyy;
inline void init(int l){
cs int len=l;
vis[1]=1;
for(int i=2;i<=len;i++){
if(!vis[i])pr[++tot]=i;
for(int j=1;j<=tot&&i*pr[j]<=len;j++){
vis[i*pr[j]]=1;
if(i-i/pr[j]*pr[j]==0)break;
}
}
}
int main(){
n=read(),m=read(),p=read();
init(m);
for(int i=1;i<=m;i++)if(vis[i])f[i-i/p*p].a0++;else f[i-i/p*p].a1++;
res[0].a0=1;
for(;n;n>>=1,Mult(f,f,f))if(n&1)Mult(res,f,res);
cout<<res[0].a1<<'
';
}