[DP 倍增] BZOJ 4818 [Sdoi2017]序列计数

怎么跟这个题这么像啊
容斥一下 用没有限制的 减去 一个质数都不选的
复杂度$O(p^2\log n)$

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
#define cl(x) memset(x,0,sizeof(x))
using namespace std;
typedef long long ll;

#define read(x) scanf("%d",&(x))

const int M=2e7+5;
int prime[M/10],num; bool vst[M];

inline void Pre(int n){
  for (int i=2;i<=n;i++){
    if (!vst[i]) prime[++num]=i;
    for (int j=1;j<=num && (ll)i*prime[j]<=n;j++){
      vst[i*prime[j]]=1;
      if (i%prime[j]==0) break;
    }
  }
}

const int N=105;
const int P=20170408;

int n,m,p;
int a[N],b[N],t[N];

inline void mul(int *a,int *b){
  for (int i=0;i<p;i++)
    for (int j=0;j<p;j++)
      (t[(i+j)%p]+=(ll)a[i]*b[j]%P)%=P;
  for (int i=0;i<p;i++) a[i]=t[i],t[i]=0; 
}

int main(){
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  read(n); read(m); read(p); Pre(m);
  if (m>=p){
    for (int i=0;i<p;i++) a[i]=(m-i)/p+1; a[0]--;
  }else{
    for (int i=1;i<=m;i++) a[i]=1;
  }
  b[0]=1;
  for (int y=n;y;y>>=1,mul(a,a))
    if (y&1)
      mul(b,a);
  int ans=b[0];
  cl(a); cl(b); b[0]=1;
  if (m>=p){
    for (int i=0;i<p;i++) a[i]=(m-i)/p+1; a[0]--;
  }else{
    for (int i=1;i<=m;i++) a[i]=1;
  }
  for (int i=1;i<=num;i++) a[prime[i]%p]--;
  for (int y=n;y;y>>=1,mul(a,a))
    if (y&1)
      mul(b,a);
  ans=(ans+P-b[0])%P;
  printf("%d\n",ans);
  return 0;
}