BZOJ 4818 [Sdoi2017]序列计数

题目大意

求长度为$n$,每个数都为不超过$m$的正整数，和为$p$的倍数，至少有一个数为质数的序列个数。

正难则反，想到可以用所有和为$p$的倍数的序列个数减去没有质数的和为$p$的倍数的序列个数。

设$f_{i,j}$为前$i$个数和$\mod p = j$的总方案数，用$cnt_i$表示$1 \sim m$中$\mod p = i$的数的个数，于是有$f_{i,j} = \sum\limits_{k=0}^{p-1} f_{i-1,k} + cnt_{(j-k) \text{mod} p}$.这个东西发现可以矩乘优化($n \le 10^9$)所以搞一下就好啦。

同理用$cnt_i$再表示一次$1 \sim m$中$\mod p = i$的合数个数，once again，减掉就可以了。

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
#define mo 20170408
int n,m,p,cnt[101];
ll ans;
struct matrix
{
	ll v[101][101];
	ll *operator[](int a){return v[a];}
	matrix(){memset(v,0,sizeof v);}
	matrix operator*(matrix x)
	{
		int i,j,k;
		matrix ret;
		for(k=0;k<p;++k)
			for(i=0;i<p;++i)
				for(j=0;j<p;++j)
					ret[i][j]=(ret[i][j]+v[i][k]*x[k][j])%mo;
		return ret;
	}
	matrix operator^(int x)
	{
		matrix ret,tmp=*this;
		for(int i=0;i<p;++i) ret.v[i][i]=1;
		while(x)
		{
			if(x&1) ret=ret*tmp;
			tmp=tmp*tmp,x>>=1;
		}
		return ret;
	}
}res,lim;
int prm[20000010],tot;
bool v[20000010];
void init()
{
	int i,j;
	for(i=2;i<=m;++i)
	{
		if(!v[i]) prm[++tot]=i;
		for(j=1;j<=tot&&i*prm[j]<=m;++j)
		{
			v[i*prm[j]]=true;
			if(i%prm[j]==0) break;
		}
	}
	for(i=1;i<=tot;++i) cnt[prm[i]%p]--;
}
int main()
{
	scanf("%d%d%d",&n,&m,&p);
	int i,j,t;
	for(i=0,t=m/p;i<p;++i) cnt[i]=t;
	for(i=1,t=m%p;i<=t;++i) cnt[i]++;
	for(i=0;i<p;++i) for(j=0;j<p;++j) lim[i][j]=cnt[(j-i+p)%p];
	memset(res.v,0,sizeof res.v),res[0][0]=1,res=res*(lim^n),ans=res[0][0];
	init();
	for(i=0;i<p;++i) for(j=0;j<p;++j) lim[i][j]=cnt[(j-i+p)%p];
	memset(res.v,0,sizeof res.v),res[0][0]=1,res=res*(lim^n),ans=(ans-res[0][0]+mo)%mo;
	printf("%lld\n",ans);
	return 0;
}