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

被空间卡的好惨啊————
参考：http://blog.csdn.net/coldef/article/details/70305596
容斥，\( 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[N];
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;
}