BZOJ4818: [Sdoi2017]序列计数

4818: [Sdoi2017]序列计数

Time Limit: 30 Sec  Memory Limit: 128 MB
Submit: 754  Solved: 461
[Submit][Status][Discuss]

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

 思路{

　　水题一道....首先容斥一下,$Ans=$所有的方案数-不含质数的方案数.

　　二维DP转移是很好想的 $ F [ i ] [ j ] $ 表示前$ i $个数和 $ %p $为$ j $的方案数.不含质数的同理

　　 $ F [ i ] [ (j+x) ] + = F [ i - 1] [ j]  $

　　这个暴力DP不对,看到n那么大,不含$Max,Min$等操作,直接矩阵快速幂了.

　　构造转移矩阵时可以第一行枚举m,求出来.

　　然后发现第二行以后都是上一行往左移动一位.

　　那么不用$ O( m*p) $的构造,直接$ O(m*m) $了.

　　在BZOJ上跑得贼慢......

}

#include<bits/stdc++.h>
#define il inline
#define RG register
#define ll long long
#define db double
#define N 110
#define mod 20170408
using namespace std;
bool vis[20000000];int P[2000000];
void pre(){
  vis[1]=true;
  for(int i=2;i<=20000000;++i){
    if(!vis[i])P[++P[0]]=i;
    for(int j=1;j<=P[0]&&P[j]*i<=20000000;++j){
      vis[i*P[j]]=true;
      if(!(i%P[j]))break;
    }
  }
}
int ss[N];
struct matrix{
  int n,m;
  ll ma[N][N];
  matrix operator *(const matrix & a)const{
    matrix c;memset(c.ma,0,sizeof(c.ma));
    if(m!=a.n)return c;c.n=n,c.m=a.m;
    for(RG int i=1;i<=c.n;++i)
      for(RG int j=1;j<=c.m;++j){
	for(RG int h=1;h<=a.n;++h){
	    c.ma[i][j]+=(ma[i][h]*a.ma[h][j])%mod;
	    if(c.ma[i][j]>=mod)c.ma[i][j]-=mod;
	  }
      }
    return c;
  }
}ans,temp;
matrix qp(matrix a,ll b){
  matrix Ans;Ans.n=Ans.m=a.n;
  for(int i=1;i<=a.n;++i)
    for(int j=1;j<=a.n;++j)
      Ans.ma[i][j]=(i==j);
  for(;b;b>>=1,a=a*a)
    if(b&1)Ans=Ans*a;
  return Ans;
}
int n,m,p;
int main(){
  scanf("%d%d%d",&n,&m,&p);pre();
  ans.n=p,ans.m=1;
  for(int i=1;i<=m;++i){
    ans.ma[(i%p)+1][1]++;
    if(ans.ma[(i%p)+1][1]>=mod)ans.ma[(i%p)+1][1]-=mod;
  }
  for(int i=1;i<=m;++i){
    int tmp=i%p;tmp=p-tmp;tmp++;
    if(tmp>p)tmp=1;
    temp.ma[1][tmp]++;
    if(temp.ma[1][tmp]>=mod)temp.ma[1][tmp]-=mod;
  }
  for(int i=2;i<=p;++i){
    for(int j=2;j<=p;++j)
      temp.ma[i][j]=temp.ma[i-1][j-1];
    temp.ma[i][1]=temp.ma[i-1][p];
  }temp.n=temp.m=p;
  temp=qp(temp,n-1);
  ans=temp*ans;
  ll Ans1=ans.ma[1][1];
  memset(ans.ma,0,sizeof(ans.ma));
  for(int i=1;i<=m;++i)
    if(vis[i]){
      ans.ma[(i%p)+1][1]++;
      if(ans.ma[(i%p)+1][1]>=mod)ans.ma[(i%p)+1][1]-=mod;
    }
  memset(temp.ma,0,sizeof(temp.ma));
  for(int i=1;i<=m;++i){
    if(!vis[i])continue;
    int tmp=i%p;tmp=p-tmp;tmp++;
    if(tmp>p)tmp=1;
    temp.ma[1][tmp]++;
    if(temp.ma[1][tmp]>=mod)temp.ma[1][tmp]-=mod;
  }
  for(int i=2;i<=p;++i){
    for(int j=2;j<=p;++j)
      temp.ma[i][j]=temp.ma[i-1][j-1];
    temp.ma[i][1]=temp.ma[i-1][p];
  }temp.n=temp.m=p;
  ans.n=p,ans.m=1;
  temp=qp(temp,n-1);
  ans=temp*ans;
  ll Ans2=ans.ma[1][1];
  cout<<(Ans1-Ans2+mod)%mod;
  return 0;
}