自动搬运

查看原文

来自洛谷，原作者为

	sky_of_war **

搬运于2025-08-24 21:51:46，当前版本为作者最后更新于2019-04-19 21:16:47，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

我觉得这篇题解下很多矩阵解法没有写清楚矩阵怎么构造...所以我来详细地写一下

同步更新于我的博客：「LUOGU P3702」[SDOI2017]序列计数

首先这个形式显然非常不好计算，容斥一波：$\text{ans}=\text{满足和为}p\text{的倍数的方案数}-\text{满足和为}p\text{的倍数且不含质数的方案数}$

我们很容易得到$20\rm pts$的优秀做法：令$f_{i,j}$表示$i$个数$\mod p$等于$j$的方案数，则转移为$f_{i,j} = \sum_{k} f_{i-1,k} \times \text{cnt}_{(j - k + p)\bmod \; p}$，其中$\text{cnt}_i$表示$1\sim m$中$\bmod P$为$i$的个数。令$g_{i,j}$表示$i$个合数$\mod p$等于$j$的方案数，则转移为$g_{i,j} = \sum_{k} g_{i-1,k} \times \text{compo}_{(j - k + p)\bmod \; p}$，其中$\text{compo}_i$表示$1\sim m$中的合数$\bmod P$为$i$的个数。 考虑用矩阵乘法优化。对于$f$，矩阵如下所示：

$$\mathbf P = \begin{bmatrix} \text{cnt}_0 & \text{cnt}_{p-1} & \text{cnt}_{p-2} & \text{cnt}_{p-3} & \text{cnt}_{p-4} & ... & \text{cnt}_1 \\ \text{cnt}_1 & \text{cnt}_0 & \text{cnt}_{p-1} & \text{cnt}_{p-2} & \text{cnt}_{p-3} & ... & \text{cnt}_2 \\ \text{cnt}_2 & \text{cnt}_1 & \text{cnt}_0 & \text{cnt}_{p-1} & \text{cnt}_{p-2} & ... & \text{cnt}_3 \\ \text{cnt}_3 & \text{cnt}_2 & \text{cnt}_1 & \text{cnt}_0 & \text{cnt}_{p-1} & ... & \text{cnt}_4 \\ ... & ... & ... & ... & ... & ... & ...\\ \text{cnt}_{p-1} & \text{cnt}_{p-2} & \text{cnt}_{p-3} & \text{cnt}_{p-4} & \text{cnt}_{p-5} & ... & \text{cnt}_0 \end{bmatrix} $$

向量长这样：

$$\mathbf V = \begin{bmatrix} \text{cnt}_0 \\ \text{cnt}_1 \\ \text{cnt}_2 \\ ... \\ \text{cnt}_{p-1} \end{bmatrix} $$

最后答案矩阵为$\mathbf F = \mathbf P ^ {n-1}\times \mathbf V$。 对于$g$，矩阵如下所示：

$$\mathbf Q = \begin{bmatrix} \text{compo}_0 & \text{compo}_{p-1} & \text{compo}_{p-2} & \text{compo}_{p-3} & \text{compo}_{p-4} & ... & \text{compo}_1 \\ \text{compo}_1 & \text{compo}_0 & \text{compo}_{p-1} & \text{compo}_{p-2} & \text{compo}_{p-3} & ... & \text{compo}_2 \\ \text{compo}_2 & \text{compo}_1 & \text{compo}_0 & \text{compo}_{p-1} & \text{compo}_{p-2} & ... & \text{compo}_3 \\ \text{compo}_3 & \text{compo}_2 & \text{compo}_1 & \text{compo}_0 & \text{compo}_{p-1} & ... & \text{compo}_4 \\ ... & ... & ... & ... & ... & ... & ...\\ \text{compo}_{p-1} & \text{compo}_{p-2} & \text{compo}_{p-3} & \text{compo}_{p-4} & \text{compo}_{p-5} & ... & \text{compo}_0 \end{bmatrix} $$

向量长这样：

$$\mathbf W = \begin{bmatrix} \text{compo}_0 \\ \text{compo}_1 \\ \text{compo}_2 \\ ... \\ \text{compo}_{p-1} \end{bmatrix} $$

最后答案矩阵为$\mathbf G = \mathbf Q ^ {n-1}\times \mathbf W$。 最终答案为$\mathbf F_{1,1} - \mathbf G_{1,1}$

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1e2 + 10, mo = 20170408, MAXM = 2e7 + 10;
struct mat
{
	int m, n, ma[MAXN][MAXN];
	mat() {}
	mat(int _m, int _n) :
		m(_m), n(_n)
	{
		memset(ma, 0, sizeof(ma));
	}
	friend inline mat operator * (mat a, mat b)
	{
		mat res = mat(a.m, b.n);
		for(int i = 1; i <= res.m; i++)
			for(int j = 1; j <= res.n; j++)
				for(int k = 1; k <= a.n; k++)
					res.ma[i][j] = (res.ma[i][j] + 1ll * a.ma[i][k]
									* b.ma[k][j] % mo) % mo;
		return res;
	}
	friend inline mat operator ^ (mat a, int b)
	{
		mat c = a, res = mat(a.m, a.n);
		for(int i = 1; i <= res.m; i++)
			res.ma[i][i] = 1;
		while(b)
		{
			if(b & 1) res = res * c;
			c = c * c;
			b >>= 1;
		}
		return res;
	}
} P, Q, V, W;
int n, m, p, cnt[MAXN], compo[MAXN];
bool prime[MAXM];
int pr[MAXM];
template <class T>
inline void _read(T &x)
{
	x = 0;
	char t = getchar();
	while(!isdigit(t) && t != '-') t = getchar();
	if(t == '-')
	{
		_read(x);
		x *= -1;
		return ;
	}
	while(isdigit(t))
	{
		x = x * 10 + t - '0';
		t = getchar();
	}
}
int ptot = 0;
int main()
{
	_read(n), _read(m), _read(p);
	for(int i = 0; i < p; i++) cnt[i] = m / p;
	for(int i = 1; i <= m % p; i++) cnt[i]++;
	P = mat(p, p);
	P.ma[1][1] = cnt[0];
	for(int i = 2; i <= p; i++) P.ma[1][i] = cnt[p - i + 1];
	for(int i = 2; i <= p; i++)
	{
		for(int j = 2; j <= p; j++)
			P.ma[i][j] = P.ma[i - 1][j - 1];
		P.ma[i][1] = P.ma[i - 1][p];
	}
	memset(prime, 1, sizeof prime);
	prime[1] = false;
	for(int i = 2; i <= m; i++)
	{
		if(prime[i] == true)
			pr[++ptot] = i;
		for(int j = 1; i * pr[j] <= m && j <= ptot; ++j)
		{
			prime[i * pr[j]] = false;
			if(i % pr[j] == 0)
				break;
		}
	}
	for(int i = 1; i <= m; i++)
		if(!prime[i]) compo[i % p]++;
	Q = mat(p, p);
	Q.ma[1][1] = compo[0];
	for(int i = 2; i <= p; i++) Q.ma[1][i] = compo[p - i + 1];
	for(int i = 2; i <= p; i++)
	{
		for(int j = 2; j <= p; j++)
			Q.ma[i][j] = Q.ma[i - 1][j - 1];
		Q.ma[i][1] = Q.ma[i - 1][p];
	}
	for(int i = 1; i <= p; i++)
		for(int j = 1; j <= p; j++)
			P.ma[i][j] %= mo, Q.ma[i][j] %= mo;
	V = mat(p, 1);
	W = mat(p, 1);
	for(int i = 1; i <= p; i++) V.ma[i][1] = cnt[i - 1] % mo;
	for(int i = 1; i <= p; i++) W.ma[i][1] = compo[i - 1] % mo;
	P = (P ^ n - 1) * V;
	Q = (Q ^ n - 1) * W;
	printf("%d\n", (P.ma[1][1] - Q.ma[1][1] + mo) % mo);
	return 0;
}