[Sdoi2017]序列计数

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;
}