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
思路
非常容易想到 O(n*m*p)的转移 不再赘述
发现不需要枚举每一个m值 而是可以在一开始就对m分类 (mod p的剩余系) 分类之后的结果就用c数组记录 这样就可以从i向j转移即 s[(i+j)%p]+=s[i]*c[j]; O(n*p*p); 还是会T
那么我们可以构造一个 s[i][j]的(p*p)的一个矩阵 表示从i向j转移的变化 那么可以定义一个初始向量 f 表示第一个位置的选择方案
明显的是 f[i]=c[i]; 即分类之后的选择方案 那么 f*s就是第二个位置选择之后每一个剩余系的方案数 现在我们可以对s做矩阵快速幂(n-1)次 最后再用f去乘s 之后的 f[0][0]就是不考虑质数的答案ans1
然后在把质数从c中抛出去 再做一遍 表示没有质数的方案数ans2 那么ans1-ans2就是答案
代码
#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
typedef long long ll;
const ll MOD=20170408;
const int maxn=200+5;
const int N=20000000+5;
int prime[N],tot;
bool visit[N];
int n,m,p;
ll s[maxn][maxn],c[maxn];
struct matrix
{
int n,m,a[maxn][maxn];//n行 m列
void clear()
{
n=m=0; memset(a,0,sizeof(a));
}
matrix operator * (const matrix &p) const
{
matrix tmp; tmp.clear();
for(int i=0; i<n; i++)
{
for(int j=0; j<p.m; j++)
{
for(int k=0; k<m; k++)
{
tmp.a[i][j]+=(1LL*a[i][k]*p.a[k][j])%MOD;
tmp.a[i][j]%=MOD;
}
}
}
tmp.n=n; tmp.m=p.m;
return tmp;
}
};
void getprime()
{
for(int i=2; i<=m; i++)
{
if(!visit[i])
{
prime[++tot]=i;
}
for(int j=1; j<=tot,prime[j]*i<=m; j++)
{
visit[prime[j]*i]=true;
if(i%prime[j]==0) break;
}
}
}
matrix quick(int n)
{
matrix e; e.clear(); e.clear();
e.n=e.m=p; matrix ans;
ans.n=ans.m=p;
for(int i=0; i<p; i++)
for(int j=0; j<p; j++)
{
e.a[i][j]=s[i][j];
}
for(int i=0; i<p; i++) ans.a[i][i]=1;
while(n)
{
if(n&1) ans=ans*e;
e=e*e;
/* for(int i=0;i<p;i++)
{
for(int j=0;j<p;j++)
{
cout<<e.a[i][j]<<' ';
}
cout<<endl;
}
puts("**********************************");*/
n>>=1;
}
return ans;
}
int main()
{
freopen("test.in","r",stdin);
freopen("test.out","w",stdout);
scanf("%d %d %d",&n,&m,&p);
for(int i=0; i<p; i++)
c[i]+=m/p;
for(int j=1; j<=m%p; j++)
c[j]++;
for(int i=0; i<p; i++)
{
for(int j=0; j<p; j++)
{
s[i][j]=c[(j-i+p)%p]%MOD;
}
}
matrix f1;
f1.n=1,f1.m=p;
for(int i=0; i<p; i++)
{
f1.a[0][i]=c[i];
}
matrix ans1=f1*quick(n-1);
getprime();
for(int i=1; i<=tot; i++)
{
// cout<<prime[i]<<endl;
c[prime[i]%p]--;
}
for(int i=0; i<p; i++)
{
c[i]=((c[i]%MOD)+MOD)%MOD;
}
matrix f2;
f2.n=1,f2.m=p;
for(int i=0; i<p; i++)
{
f2.a[0][i]=c[i];
}
for(int i=0; i<p; i++)
{
for(int j=0; j<p; j++)
{
s[i][j]=c[(j-i+p)%p]%MOD;
// cout<<s[i][j]<<' ';
}
// cout<<endl;
}
matrix ans2=f2*quick(n-1);
cout<<(ans1.a[0][0]-ans2.a[0][0]+MOD)%MOD<<endl;
return 0;
}