[BZOJ4818][Sdoi2017][容斥原理][矩陣優化DP]序列計數

考慮容斥原理
Ans=f滿足和為p的倍數−f滿足和為p的倍數求不含質數

可以DP,f(i,j)表示轉移到第i位，前i位和模P等於j的方案數

那麼顯然f(i,j)=∑f(i−1,k)∗cnt(j−k+p)modp

其中cnti表示1~m中模P為i的個數(計算第二個個數的時候要出去質數)。直接轉移當然不行，進一步可以發現這個DP可以用矩陣乘法來優化，就用一般套路就行。

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <string>
#define N 110
#define P 20170408
#define M 20000010

using namespace std;

int n,m,p;
char vis[M];
int prime[1280000];
int cnt[N];

inline void Add(int &x,int y){
  if((x+=y)>=P) x-=P;
}

struct Mat{
  int a[N][N];
  Mat(){ for(int i=0;i<p;i++)for(int j=0;j<p;j++)a[i][j]=0; }
  int *operator [](int x){ return a[x]; }
  friend Mat operator *(Mat A,Mat B){
    Mat C;
    for(int i=0;i<p;i++)
      for(int j=0;j<p;j++)
    for(int k=0;k<p;k++)
      Add(C[i][j],1ll*A[i][k]*B[k][j]%P);
    return C;
  }
}f1,f2,g;

inline Mat Pow(Mat A,int c){
  Mat ret;
  for(int i=0;i<p;i++) ret[i][i]=1;
  for(;c;c>>=1,A=A*A) if(c&1) ret=ret*A;
  return ret;
}

int main(){
  freopen("1.in","r",stdin);
  freopen("1.out","w",stdout);
  scanf("%d%d%d",&n,&m,&p);
  for(int i=1;i<=m;i++) cnt[i%p]++;
  for(int i=0;i<p;i++)
    for(int j=0;j<p;j++)
      g[i][j]=cnt[(i-j+p)%p];
  f1[0][0]=f2[0][0]=1;
  f1=f1*Pow(g,n);
  vis[1]=1;
  for(int i=2;i<=m;i++){
    if(!vis[i]) prime[++*prime]=i;
    for(int j=1;j<=*prime&&1ll*prime[j]*i<=m;j++)
      if(vis[prime[j]*i]=1,i%prime[j]==0) break;
  }
  memset(cnt,0,sizeof(cnt));
  for(int i=1;i<=m;i++)
    if(vis[i]) cnt[i%p]++;
  for(int i=0;i<p;i++)
    for(int j=0;j<p;j++)
      g[i][j]=cnt[(i-j+p)%p];
  f2=f2*Pow(g,n);
  printf("%d\n",(f1[0][0]-f2[0][0]+P)%P);
  //printf("%d\n",clock());
  return 0;
}