bzoj 4818: [Sdoi2017]序列计数【容斥原理+dp+矩阵乘法】

被空间卡的好惨啊————
参考：https://www.geek-share.com/detail/2703694923.html
容斥，\( ans=ans_{没有限制}-ans{没有质数} \)
动规递推式，\( f[i][j]=\sum_{k=0}^{p-1}f[i-1][k]*cnt[(i-j+p)%p] \)，\( cnt[i] \)表示\( %p==i \)的数，注意计算第二个\( ans \)的时候要用筛子去掉质数
因为\( n\leq 10^9 \)，所以选择矩阵乘法加速递推式。

#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
const long long P=105,N=20000005,mod=20170408;
long long n,m,p,cnt[P],q[1280000];
bool v
;
struct qwe
{
long long a[P][P];
qwe operator * (qwe b)
{
qwe c;
for(long long i=0;i<p;i++)
for(long long j=0;j<p;j++)
{
c.a[i][j]=0;
for(long long k=0;k<p;k++)
c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j]%mod)%mod;
}
return c;
}
}f1,f2,g;
qwe ksm(qwe a,long long b)
{
qwe r;
for(long long i=0;i<p;i++)
r.a[i][i]=1;
while(b)
{
if(b&1)
r=r*a;
a=a*a;
b>>=1;
}
return r;
}
int main()
{
scanf("%lld%lld%lld",&n,&m,&p);
for(long long i=1;i<=m;i++)
cnt[i%p]++;
for(long long i=0;i<p;i++)
for(long long j=0;j<p;j++)
g.a[i][j]=cnt[(i-j+p)%p];
f1.a[0][0]=f2.a[0][0]=1;
f1=f1*ksm(g,n);
v[1]=1;
for(long long i=2;i<=m;i++)
{
if(!v[i])
q[++q[0]]=i;
for(long long j=1;j<=q[0]&&i*q[j]<=m;j++)
{
v[i*q[j]]=1;
if(i%q[j]==0)
break;
}
}
memset(cnt,0,sizeof(cnt));
for(long long i=1;i<=m;i++)
if(v[i])
cnt[i%p]++;//,cout<<i<<endl;;
for(long long i=0;i<p;i++)
for(long long j=0;j<p;j++)
g.a[i][j]=cnt[(i-j+p)%p];
f2=f2*ksm(g,n);//cout<<f1.a[0][0]<<" "<<f2.a[0][0]<<endl;
printf("%lld\n",(f1.a[0][0]-f2.a[0][0]+mod)%mod);
return 0;
}