526互联

【动态规划】【SDOI2017】序列计数

发布时间 2023-09-03 21:58:54作者: The_Last_Candy

【动态规划】【SDOI2017】序列计数

题目描述

Alice 想要得到一个长度为 \(n\) 的序列，序列中的数都是不超过 \(m\) 的正整数，而且这 \(n\) 个数的和是 \(p\) 的倍数。

Alice 还希望，这 \(n\) 个数中，至少有一个数是质数。

Alice 想知道，有多少个序列满足她的要求。

\(1 \leq n \leq 10^9,1 \leq m \leq 2 \times 10^7,1 \leq p \leq 100\)。

算法概述

先不考虑质数的限制，对第一条限制做 \(dp\) ，容易得到：

设 \(f_{i,j}\) 为前 \(i\) 个数，模 \(p\) 的值为 \(j\) 的方案数。

\[f_{i,j} = \sum_kf_{i - 1,k} \times cntf_{(j - k + p) \% p} \]

其中 \(cntf_i\) 为 \(1 \to m\) 中模 \(p\) 为 \(i\) 的数字个数。

同理，考虑质数限制，设 \(g_{i,j}\) 为前 \(i\) 个数，只取合数，模 \(p\) 的值为 \(j\) 的方案数。

\[g_{i,j} = \sum_kg_{i - 1,k} \times cntg_{(j - k + p) \% p} \]

其中 \(cntg_i\) 为 \(1 \to m\) 的合数中模 \(p\) 为 \(i\) 的数字个数。

发现 \(p\) 比较小，并且 \(j = (k + (j - k + p)) \% p\) ，所以我们考虑矩阵乘法优化这个 \(dp\) 。

考虑从 \(i - 1\) 转移到 \(i\) 的时候，\((x + y) \% p = z\) 的项都会转移到 \(z\) 上，考虑这样设计矩阵：

\[\begin{bmatrix} cntf_0&cntf_1&cntf_2&\dots&cntf_{p - 1}\\ cntf_{p - 1}&cntf_0&cntf_1&\dots&cntf_{p - 2}\\ cntf_{p - 2}&cntf_{p - 1}&cntf_0&\dots&cntf_{p - 3}\\ \dots&\dots&\dots&\dots&\dots\\ cntf_1&cntf_2&cntf_3&\dots&cntf_0 \end{bmatrix} \]

初始矩阵为

\[\begin{bmatrix} cntf_0&cntf_1&cntf_2&\dots&cntf_{p - 1} \end{bmatrix} \]

如果不清楚为什么是这样，建议将两个矩阵手模一下，观察乘法过程。

\(f_{n,0}\) 就是初始矩阵乘以转移矩阵的 \(n - 1\) 次方的第 \(0\) 项。

\(g_{n,0}\) 的求法同理。

根据补集转化，至少有一个质数的答案就是 \(f_{n,0} - g_{n,0}\) 。

Code

#include<bits/stdc++.h>
using namespace std;
const int N = 2e6 + 5,MOD = 20170408;
int prime[N],cnt = 0,p,cntf[105],cntg[105],m,n;
typedef long long ll;
ll f[105],g[105];
bitset <N * 10> vis;
inline void init()
{
	cntf[1]++;cntg[1]++;
	for(int i = 2;i <= m;i++)
	{
		cntf[i % p]++;
		if(!vis[i]) prime[++cnt] = i;
		else cntg[i % p]++;
		for(int j = 1;j <= cnt && 1ll * i * prime[j] <= m;j++)
		{
			vis[i * prime[j]] = 1;
			if(!(i % prime[j])) break;
		}
	}
}
struct Matrix {
	ll a[105][105];
};
Matrix operator *(Matrix x,Matrix y)
{
	Matrix z;
	memset(z.a,0,sizeof(z.a));
	for(int i = 0;i <= p - 1;i++)
		for(int j = 0;j <= p - 1;j++)
			for(int k = 0;k <= p - 1;k++)
				z.a[i][j] = (z.a[i][j] + x.a[i][k] * y.a[k][j] % MOD) % MOD;
	return z;
}
inline Matrix ksm(Matrix base,int pts)
{
	Matrix ret;
	memset(ret.a,0,sizeof(ret.a));
	for(int i = 0;i <= p - 1;i++) ret.a[i][i] = 1;
	for(;pts > 0;pts >>= 1,base = base * base)
		if(pts & 1)
			ret = ret * base;
	return ret;
}
int main()
{
	cin>>n>>m>>p;
	init();
	Matrix kf,kg;
	for(int i = 0;i <= p - 1;i++)
		for(int j = 0;j <= p - 1;j++)
			kf.a[i][j] = cntf[(j - i + p) % p],kg.a[i][j] = cntg[(j - i + p) % p];
	kf = ksm(kf,n - 1);
	kg = ksm(kg,n - 1);
	for(int i = 0;i <= p - 1;i++) 
		for(int j = 0;j <= p - 1;j++)
			f[i] = (f[i] + cntf[j] * kf.a[j][i] % MOD) % MOD,g[i] = (g[i] + cntg[j] * kg.a[j][i] % MOD) % MOD;
	cout<<(f[0] - g[0] + MOD) % MOD;
	return 0;
}