【洛谷5366】[SNOI2017] 遗失的答案（状压DP）

这个我省三年前的省选真题我肝了将近一天 ……

完了我感觉我省选要凉了怎么办在线等急

虽然我还不知道什么时候省选

这题我写了两种做法。

我这么懒的人怎么会写两种做法呢？

因为这题数据范围太！诡！异！了！

下来再说。

题目

洛谷 5366

分析

显然，如果 (L) 不是 (G) 的倍数、(X>N) 、(X) 不是 (L) 的因数或 (X) 不是 (G) 的倍数则无解，否则可以把题转换成求从 ([1,lfloorfrac{N}{G} floor]) 中选出若干个数，使得最大公因数是 (1) 且最小公倍数是 (frac{L}{G}) 的方案中，有多少种方案选了 (frac{X}{G}) 。以下的 (L) 均为原题中的 (frac{L}{G}) ，(X) 均为原题中的 (frac{X}{G}) 。

令 (x=prod_j p_j^{a_j}) （称为 (x) 的标准分解），其中 (p_j) 是质数，那么有：

[gcd({x_i})=prod_j p_j^{min_i({a_{i,j}})} ]

[mathrm{lcm}({x_i})=prod_j p_j^{max_i({a_{i,j}})} ]

（这个应该非常显然吧）

因此，在最终选出的数的标准分解中，每个质数的最小指数应该是 (0) ，最大指数应该等于 (L) 的标准分解中该质数的指数。

不超过 (10^8) 的数最多只能有 (8) 个不同的质因数，所以考虑设计状态 (S) ，每个质数用两位分别表示是否已经选过标准分解中该质数的指数为 (0) 和指数与 (L) 的指数相等的数。类似地，给每个数确定一个属性表示能满足哪些条件，最终状态就是选出的数的属性或起来。

显然，只有 (L) 的不超过 (2sqrt{L}) 个因数是可以选的，因此可以直接暴力枚举 (L) 的每个质因子来确定这些因数的属性。（我想了一个多小时怎么给 (10^8) 个数确定属性未果后才意识到只有 (2sqrt{L}) 个有用 …… 菜死了。）虽然属性看似有 (2^{8 imes 2}=65536) 种，但事实上由于一个数不可能同时满足一个质数的两个条件（也就是不存在一个质数的两位都是 (1) 的属性），而只有当 (L) 的标准分解中这个质数的指数不小于 (2) 时才可能存在一个数同时不满足两个条件（即一个质数的两位都是 (0) ）。因此事实上不同的属性数是很少的。通过如下 Python 代码发现最多可能只有 (576) 个。

（为什么说可能呢？因为我真的不会写 Python ，怎么循环都是现查的。我知道下面这个代码很丑很菜看起来像刚学两天编程的人写的，甚至怀疑正确性不太对。）

ans = 0
p = [2, 3, 5, 7, 11, 13, 17, 19]
now = 1
pos = 8
N = 1e8

for i in p :
	now = now * i
i = 0
while i < pos :
	ans = max(ans, 3 ** i * 2 ** (pos - i))
	now = now / p[i]
	now = now * p[i] * p[i]
	while now > N :
		pos = pos - 1
		now = now / p[pos]
	i = i + 1
print(ans)

接下来有两条路可以走。

方法一

这个做法思路清晰解法自然，所以我想到就直接写了。

将各种属性按任意顺序排成一个序列，设 (mathrm{pre}_{i,S}) 表示前 (i) 种属性的数组成了集合 (S) 的状态，(mathrm{suf}_{i,S}) 表示后 (i) 种属性的数组成了集合 (S) 的状态。这两个数组可以直接 DP 求出。一种属性贡献的方案数是 (2^k-1) ，其中 (k) 是这种属性的数的个数，减 (1) 是减去这些数都不取的情况。注意特判属性 (0) 即使不取也对 (S) 没有影响，方案数是 (2^k) 。

考虑预处理钦定第 (i) 种属性的某个数必须选的情况的答案。首先把 (mathrm{pre}_{i-1}) 和 (mathrm{suf}_{i+1}) 合并起来。观察发现这一步是一个或卷积，可以用 FWT 优化。然后再考虑第 (i) 种属性的没有被钦定的数，一共有 (2^{k-1}) 种选法。最后，如果所有没有被钦定的数组成的状态或上第 (i) 种属性可成为全 (1) 的状态，则可计入答案。询问时求出 (X) 的属性然后直接输出该属性被钦定一个数时的答案即可。代码见下文。

这个做法的时间复杂度 (O(Mp2^{2p})) ，其中 (M) 是属性的种类数，(p) 是 (L) 的质因数种类数。最坏情况 (M=576,p=8) ，算下来大概 (6 imes 10^8) 。

虽然网上题解都说常数小可以过，但是我花了一晚上 + 一上午时间使尽浑身解数卡常后洛谷上还！是！过！不！去！

我也不知道某 s 姓选手是怎么带了个火车头就卡进去的。

方法二

这个做法思路更清晰解法更自然，我也不知道我为什么想不到，可能就是太菜了吧。

先介绍两个小技巧

子集反演

若

[f_S=sum_{Tsubset S}g_T ]

则

[g_S=sum_{Tsubset S}(-1)^{(|S|-|T|)}f(T) ]

等等这不就是容斥原理吗 …… 这都有名字系列.jpg 。

子集求和

[f_S=sum_{Tsubset S}g_T ]

在已知 (g) 的情况下求 (f) 。

当然可以直接暴力枚举子集，时间复杂度是 (O(3^n)) 的，其中 (n) 是集合大小。但如下代码可以做到 (O(n2^n)) :

memcpy(f, g, sizeof(int[1 << n]));
for (int i = 0; i < n; i++)
	for (int j = 0; j < (1 << n); j++)
		if (j & (1 << i))
			f[j] += f[j ^ (1 << i)];

我第一次看到这个技巧的时候震惊了，竟然改变两层循环的顺序就能防止重复计算 ……

简单证明一下。设 (f_{i,S}) 表示前 (i) 位是 (S) 的子集，其他位和 (S) 一样的集合的 (g) 之和。那么从 (i-1) 转移到 (i) 时，新增的就是前 (i) 位是 (S) 的子集，第 (i) 位是 (0) （即跟 (S) 的第 (i) 位不一样）且其他位和 (S) 一样的集合之和，也就是 (f_{i,Ssetminus{i}}) 。滚动存储一下就是上面的代码。

回到这道题。仍然是预处理钦定第 (i) 种属性的某个数必须选的情况的答案。用 (f_S) 表示钦定第 (i) 种属性的某个数必须选时选出状态 (S) 的方案数。由于第 (i) 种属性已经钦定选了一个，所以只有当 (S) 是该属性的超集时才有意义。

从所有有意义的 (S) 中删去属性 (i) ，就变成了 (i) 的补集的全体子集，因此可以直接套用子集卷积的结论，即（设属性 (i) 为 (P_i) ）：

[f_{Scup P_i}=sum_{Tsubset S}(-1)^{(|S|-|T|)}g(Tcup P_i) ]

此时 (g_S) 的定义应该是「钦定第 (i) 种属性的某个数必须选时选出状态 (S) 或其子集的方案数」。那么 (g_S=2^{k-1}) ，其中 (k) 是属性为 (S) 的子集的数的个数（显然钦定的那个数也在其中，所以要减掉），可以用上述子集求和的方法预处理出。

处理询问就和方法一一样了。时间复杂度 (O(3^{2p})) ，其中 (p) 是 (L) 的质因数种类数。

代码

方法一

LOJ 上 1372 ms ，洛谷 TLE 一个点。

#include <algorithm>
#include <cstdio>
#include <cctype>
#include <cstring>
using namespace std;

namespace zyt
{
	template<typename T>
	inline bool read(T &x)
	{
		char c;
		bool f = false;
		x = 0;
		do
			c = getchar();
		while (c != EOF && c != '-' && !isdigit(c));
		if (c == EOF)
			return false;
		if (c == '-')
			f = true, c = getchar();
		do
			x = x * 10 + c - '0', c = getchar();
		while (isdigit(c));
		if (f)
			x = -x;
		return true;
	}
	template<typename T>
	inline void write(T x)
	{
		static char buf[20];
		char *pos = buf;
		if (x < 0)
			putchar('-'), x = -x;
		do
			*pos++ = x % 10 + '0';
		while (x /= 10);
		while (pos > buf)
			putchar(*--pos);
	}
	typedef long long ll;
	const int N = 600, S = 1 << 16, PR = 10, P = 1e9 + 7, INV2 = (P + 1) / 2;
	int n, g, l, q, pow[N], p[PR], pk[PR], pcnt, val[N], id[S], num[N], cnt;
	int pre[N][S], suf[N][S], tot[N][S], ans[N];
	int power(int a, int b)
	{
		int ans = 1;
		while (b)
		{
			if (b & 1)
				ans = (ll)ans * a % P;
			a = (ll)a * a % P;
			b >>= 1;
		}
		return ans;
	}
	void get_prime(int n)
	{
		for (int i = 2; i * i <= n; i++)
			if (n % i == 0)
			{
				p[pcnt] = i, pk[pcnt] = 1;
				while (n % i == 0)
					n /= i, pk[pcnt] *= i;
				++pcnt;
			}
		if (n > 1)
			p[pcnt] = pk[pcnt] = n, pcnt++;
	}
	int get_type(const int n)
	{
		int ans = 0;
		for (int i = 0; i < pcnt; i++)
		{
			if (n % pk[i] == 0)
				ans |= (1 << i);
			else if (n % p[i])
				ans |= (1 << (i + pcnt));
		}
		return ans;
	}
	namespace FWT
	{
		void fwt(int *const a, const int n)
		{
			for (int l = 1; l < n; l <<= 1)
				for (int i = 0; i < n; i += (l << 1))
					for (int k = i; k < i + l; k++)
						a[l + k] = (a[l + k] + a[k]) % P;
		}
		void ifwt(int *const a, const int n)
		{
			for (int l = 1; l < n; l <<= 1)
				for (int i = 0; i < n; i += (l << 1))
					for (int k = i; k < i + l; k++)
						a[l + k] = (a[l + k] - a[k] + P) % P;
		}
	}
	int work()
	{
		using namespace FWT;
		read(n), read(g), read(l), read(q);
		if (l % g)
		{
			while (q--)
				puts("0");
			return 0;
		}
		l /= g, n /= g;
		get_prime(l);
		for (int i = 1; i * i <= l && i <= n; i++)
			if (l % i == 0)
			{
				int t = get_type(i);
				if (!id[t])
					id[t] = cnt, num[cnt] = t, val[cnt] = 1, cnt++;
				val[id[t]] = val[id[t]] * 2 % P;
				if (i * i < l && l / i <= n)
				{
					t = get_type(l / i);
					if (!id[t])
						id[t] = cnt, num[cnt] = t, val[cnt] = 1, cnt++;
					val[id[t]] = val[id[t]] * 2 % P;
				}
			}
		pre[0][0] = 1, pre[0][num[0]] += val[0] - 1;
		fwt(pre[0], 1 << (pcnt * 2));
		for (int i = 1; i < cnt; i++)
		{
			if (!num[i])
				pre[i][0] = val[i];
			else
				pre[i][0] = 1, pre[i][num[i]] = val[i] - 1;
			fwt(pre[i], 1 << (pcnt * 2));
			for (int j = 0; j < (1 << (pcnt * 2)); j++)
				pre[i][j] = (ll)pre[i][j] * pre[i - 1][j] % P;
		}
		suf[cnt - 1][0] = 1, suf[cnt - 1][num[cnt - 1]] += val[cnt - 1] - 1;
		fwt(suf[cnt - 1], 1 << (pcnt * 2));
		for (int i = cnt - 2; i >= 0; i--)
		{
			if (!num[i])
				suf[i][0] = val[i];
			else
				suf[i][0] = 1, suf[i][num[i]] = val[i] - 1;
			fwt(suf[i], 1 << (pcnt * 2));
			for (int j = 0; j < (1 << (pcnt * 2)); j++)
				suf[i][j] = (ll)suf[i][j] * suf[i + 1][j] % P;
		}
		for (int i = 0; i < cnt; i++)
		{
			static int tmp[S];
			if (!i)
				memcpy(tmp, suf[i + 1], sizeof(int[1 << (pcnt * 2)]));
			else if (i == cnt - 1)
				memcpy(tmp, pre[i - 1], sizeof(int[1 << (pcnt * 2)]));
			else
				for (int j = 0; j < (1 << (pcnt * 2)); j++)
					tmp[j] = (ll)pre[i - 1][j] * suf[i + 1][j] % P; 
			ifwt(tmp, 1 << (pcnt * 2));
			for (int j = 0; j < (1 << (pcnt * 2)); j++)
			{
				tot[i][j | num[i]] = (tot[i][j | num[i]] + tmp[j]) % P;
				if ((j | num[i]) == (1 << (pcnt * 2)) - 1)
					ans[i] = (ans[i] + tot[i][j]) % P;
			}
			ans[i] = (ll)ans[i] * val[i] % P * INV2 % P;
		}
		while (q--)
		{
			int x;
			read(x);
			if (x % g || (x / g) > n || l % (x / g))
				write(0);
			else
			{
				int t = get_type(x / g);
				write(ans[id[t]]);
			}
			putchar('
');
		}
		return 0;
	}
}
int main()
{
	return zyt::work();
}

方法二

LOJ 上 58 ms ，洛谷上 201 ms 。

#include <algorithm>
#include <cstdio>
#include <cctype>
#include <cstring>
using namespace std;

namespace zyt
{
	template<typename T>
	inline bool read(T &x)
	{
		char c;
		bool f = false;
		x = 0;
		do
			c = getchar();
		while (c != EOF && c != '-' && !isdigit(c));
		if (c == EOF)
			return false;
		if (c == '-')
			f = true, c = getchar();
		do
			x = x * 10 + c - '0', c = getchar();
		while (isdigit(c));
		if (f)
			x = -x;
		return true;
	}
	template<typename T>
	inline void write(T x)
	{
		static char buf[20];
		char *pos = buf;
		if (x < 0)
			putchar('-'), x = -x;
		do
			*pos++ = x % 10 + '0';
		while (x /= 10);
		while (pos > buf)
			putchar(*--pos);
	}
	typedef long long ll;
	const int N = 600, S = 1 << 16, PR = 10, P = 1e9 + 7, SQRT= 1e4 + 10;
	int n, g, l, q, p[PR], pk[PR], pcnt, val[S], id[S], num[N], cnt, pow[SQRT], count[S], ans[N];
	void get_prime(int n)
	{
		for (int i = 2; i * i <= n; i++)
			if (n % i == 0)
			{
				p[pcnt] = i, pk[pcnt] = 1;
				while (n % i == 0)
					n /= i, pk[pcnt] *= i;
				++pcnt;
			}
		if (n > 1)
			p[pcnt] = pk[pcnt] = n, pcnt++;
	}
	int get_type(const int n)
	{
		int ans = 0;
		for (int i = 0; i < pcnt; i++)
		{
			if (n % pk[i] == 0)
				ans |= (1 << i);
			else if (n % p[i])
				ans |= (1 << (i + pcnt));
		}
		return ans;
	}
	void init()
	{
		pow[0] = 1;
		for (int i = 1; i < SQRT; i++)
			pow[i] = (ll)pow[i - 1] * 2 % P;
		for (int i = 0; i < S; i++)
			count[i] = count[i >> 1] + (i & 1);
	}
	int work()
	{
		init();
		read(n), read(g), read(l), read(q);
		if (l % g)
		{
			while (q--)
				puts("0");
			return 0;
		}
		l /= g, n /= g;
		get_prime(l);
		for (int i = 1; i * i <= l && i <= n; i++)
			if (l % i == 0)
			{
				int t = get_type(i);
				if (!id[t])
					id[t] = cnt, num[cnt] = t, cnt++;
				++val[t];
				if (i * i < l && l / i <= n)
				{
					t = get_type(l / i);
					if (!id[t])
						id[t] = cnt, num[cnt] = t, cnt++;
					++val[t];
				}
			}
		for (int i = 0; i < pcnt * 2; i++)
			for (int j = 0; j < (1 << (pcnt * 2)); j++)
				if (j & (1 << i))
					val[j] += val[j ^ (1 << i)];
		for (int i = 0; i < cnt; i++)
		{
			int t = ((1 << (pcnt * 2)) - 1) ^ num[i];
			for (int j = t; j; j = ((j - 1) & t))
				ans[i] = ((ll)ans[i] + (count[j | num[i]] & 1 ? -1 : 1) * pow[val[j | num[i]] - 1] + P) % P;
			ans[i] = ((ll)ans[i] + (count[num[i]] & 1 ? -1 : 1) * pow[val[num[i]] - 1] + P) % P;
		}
		while (q--)
		{
			int x;
			read(x);
			if (x % g || (x / g) > n || l % (x / g))
				write(0);
			else
			{
				int t = get_type(x / g);
				write(ans[id[t]]);
			}
			putchar('
');
		}
		return 0;
	}
}
int main()
{
	return zyt::work();
}