[BZOJ4818][Sdoi2017][容斥原理][矩阵优化DP]序列计数

考虑容斥原理
$Ans=f_{满足和为p的倍数}-f_{满足和为p的倍数求不含质数}$
可以DP,$f(i,j)$表示转移到第i位，前i位和模P等于j的方案数

那么显然$f(i,j)=\sum f(i-1,k)*cnt_{(j-k+p)\mod p}$其中$cnt_i$表示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;
}