把博客园图标替换成自己的图标
把博客园图标替换成自己的图标end

【BZOJ4818】[SDOI2017] 序列计数(矩乘水题)

点此看题面

大致题意: 问有多少个长度为\(n\)、其中都是不超过\(m\)的正整数的序列,有至少一个质数,且这\(n\)个数的和是\(p\)的倍数。

前言

一道挺水的题目,然而依然能码出三个致命错误的我还是太菜了。。。

容斥

比较显然,有至少一个质数的答案就是总答案减去不含质数的答案。

于是这道题就简单了许多。

动态规划+矩阵优化

考虑设\(f_{i,j}\)表示选择了\(i\)个数,模\(p\)的余数为\(j\)的方案数,则每次只要枚举一个可选的数就能进行转移。

而看到\(n\)这么大,\(p\)却只有\(100\),很自然就想到矩乘。

建出转移矩阵,其中第\(i\)行第\(j\)列表示从余数为\(i\)转移到余数为\(j\)的方案数,显然答案就是转移矩阵的\(n\)次幂第\(0\)行第\(0\)列的值。

在建立总答案的转移矩阵时,可以枚举余数\(i\),计算出\(1\sim m\)中模\(p\)\(i\)的数的个数就是\(\frac{m-i}p+[i\not=0]\)(之所以\(i\not=0\)时才加\(1\),是因为\(0\)不是正整数)。然后将每个\((j,(i+j)\%p)\)加上这一个数。

在建立不含质数的转移矩阵时,可以直接枚举质数,在总答案转移矩阵的基础上减去贡献。

具体实现详见代码。

代码

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define X 20170408
#define Inc(x,y) ((x+=(y))>=X&&(x-=X))
using namespace std;
int n,m,p;
class LinearSieve
{
	private:
		#define SZ 20000000
		int P[SZ+5];
	public:
		int Pt;I int operator [] (CI x) {return P[x];}
		I LinearSieve()//线性筛筛质数
		{
			for(RI i=2,j;i<=SZ;++i) for(!P[i]&&(P[++Pt]=i),
				j=1;j<=Pt&&i*P[j]<=SZ;++j) if(P[i*P[j]]=1,!(i%P[j])) break;
		}
}S;
struct M//矩阵
{
	#define P 100
	int a[P+5][P+5];I int *operator [] (CI x) {return a[x];}
	I M() {for(RI i=0,j;i^p;++i) for(j=0;j^p;++j) a[i][j]=0;}
	I friend M operator * (Con M& x,Con M& y)//矩阵乘法
	{
		M t;for(RI i=0,j,k;i^p;++i) for(j=0;j^p;++j)
			for(k=0;k^p;++k) t[i][j]=(1LL*x.a[i][k]*y.a[k][j]+t[i][j])%X;return t;
	}
	I friend M operator ^ (M x,RI y)//矩阵快速幂
	{
		M t;for(RI i=0;i^p;++i) t[i][i]=1;W(y) y&1&&(t=t*x,0),x=x*x,y>>=1;return t;
	}
}U;
int main()
{
	RI i,j,t,ans;scanf("%d%d%d",&n,&m,&p);
	for(i=0;i^p&&i<=m;++i) for(t=(m-i)/p+(i!=0),j=0;j^p;++j) Inc(U[j][(i+j)%p],t);//计算总答案
	for(ans=(U^n)[0][0],i=1;i<=S.Pt&&S[i]<=m;++i)
		for(j=0;j^p;++j) !U[j][(S[i]+j)%p]--&&(U[j][(S[i]+j)%p]+=X);//计算不含质数答案
	return ans=(ans-(U^n)[0][0]+X)%X,printf("%d\n",ans),0;
}
posted @ 2020-05-26 07:55  TheLostWeak  阅读(558)  评论(0编辑  收藏  举报