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

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

struct Mat{
int a

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