[SDOI2017]序列计数

也许更好的阅读体验

(mathcal{Description})

给你一个长度为(n)的序列，序列中的每个数都是不超过(m)的正整数，求满足以下两个条件的序列数量：

1.序列中至少有一个质数

2.序列中(n)个数之和为(p)的倍数

​​
数据范围：(1leq nleq 10^9,1leq mleq 2∗10^7,1leq pleq 100)

(mathcal{Solution})

有两个条件，第一个条件可以先算出所有的序列数量然后减去没有质数的序列数量
第二个条件则是一个比较套路的(DP)了，通常(p)都会很小

设(f_{i,j})表示长度为(i)的序列的每个数之和在模(p)意义下为(j)的序列数量
考虑把长度(i)的序列分成两部分,然后把这两部分加起来模(p)后要得到(j)
可以得到(f_{i,j}=f_{i-k,left(j-h+p ight)\% p}*f_{k,h})，这个(k,h)任取
于是我们取(k=frac{i}{2})，这样就可以类似快速幂倍增地计算了

而对于(f)的初值，有(f_{1,i}=m/p+(m\%p>=i) - (i==0))

(mathcal{Code})

/*******************************
Author:Morning_Glory
LANG:C++
Created Time:2019年10月06日 星期日 07时54分22秒
*******************************/
#include <cstdio>
#include <fstream>
#include <cstring>
using namespace std;
const int maxn = 20000007;
const int maxp = 105;
const int mod = 20170408;
//{{{cin
struct IO{
	template<typename T>
	IO & operator>>(T&res){
		res=0;
		bool flag=false;
		char ch;
		while((ch=getchar())>'9'||ch<'0')	flag|=ch=='-';
		while(ch>='0'&&ch<='9')	res=(res<<1)+(res<<3)+(ch^'0'),ch=getchar();
		if (flag)	res=~res+1;
		return *this;
	}
}cin;
//}}}
int n,m,p,cnt;
int prime[maxn],f[maxn],g[maxn];
bool vis[maxn];
//{{{sieve
void sieve (int n)
{
	for (int i=2;i<=n;++i){
		if (!vis[i])	vis[prime[++cnt]=i]=true;
		for (int j=1;j<=cnt&&prime[j]*i<=n;++j){
			vis[prime[j]*i]=true;
			if (i%prime[j]==0)	break;
		}
	}
}
//}}}
//{{{mul
void mul (int *f,int *g)
{
	int t[maxp];
	for (int i=0;i<p;++i)	t[i]=0;
	for (int i=0;i<p;++i)
		for (int j=0;j<p;++j)
			t[i]=(t[i]+1ll*f[j]*g[(i-j+p)%p]%mod)%mod;
	for (int i=0;i<p;++i)	f[i]=t[i];
}
//}}}
//{{{ksm
void ksm (int *f,int b)
{
	int s[maxp];
	for (int i=0;i<p;++i)	s[i]=0;
	s[0]=1;
	for (;b;b>>=1,mul(f,f))
		if (b&1)	mul(s,f);
	for (int i=0;i<p;++i)	f[i]=s[i];
}
//}}}
int main()
{
	cin>>n>>m>>p;
	sieve(m);
	for (int i=0;i<p;++i)	f[i]=g[i]=m/p+(m%p>=i);
	--f[0],--g[0];
	for (int i=1;i<=cnt;++i)	--g[prime[i]%p];
	ksm(f,n),ksm(g,n);
	printf("%d
",(f[0]-g[0]+mod)%mod);
	return 0;
}

如有哪里讲得不是很明白或是有错误，欢迎指正
如您喜欢的话不妨点个赞收藏一下吧