【洛谷5366】[SNOI2017] 遗失的答案（状压DP）
这个我省三年前的省选真题我肝了将近一天 ……
完了我感觉我省选要凉了怎么办在线等急
虽然我还不知道什么时候省选
这题我写了两种做法。
我这么懒的人怎么会写两种做法呢？
因为这题数据范围太！诡！异！了！
下来再说。
题目
分析
显然，如果 \(L\) 不是 \(G\) 的倍数、\(X>N\) 、\(X\) 不是 \(L\) 的因数或 \(X\) 不是 \(G\) 的倍数则无解，否则可以把题转换成求从 \([1,\lfloor\frac{N}{G}\rfloor]\) 中选出若干个数，使得最大公因数是 \(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\) 是质数，那么有：
（这个应该非常显然吧）
因此，在最终选出的数的标准分解中，每个质数的最小指数应该是 \(0\) ，最大指数应该等于 \(L\) 的标准分解中该质数的指数。
不超过 \(10^8\) 的数最多只能有 \(8\) 个不同的质因数，所以考虑设计状态 \(S\) ，每个质数用两位分别表示是否已经选过标准分解中该质数的指数为 \(0\) 和指数与 \(L\) 的指数相等的数。类似地，给每个数确定一个属性表示能满足哪些条件，最终状态就是选出的数的属性或起来。
显然，只有 \(L\) 的不超过 \(2\sqrt{L}\) 个因数是可以选的，因此可以直接暴力枚举 \(L\) 的每个质因子来确定这些因数的属性。（我想了一个多小时怎么给 \(10^8\) 个数确定属性未果后才意识到只有 \(2\sqrt{L}\) 个有用 …… 菜死了。）虽然属性看似有 \(2^{8\times 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\times 10^8\) 。
虽然网上题解都说常数小可以过，但是我花了一晚上 + 一上午时间使尽浑身解数卡常后洛谷上还！是！过！不！去！
我也不知道某 s 姓选手是怎么带了个火车头就卡进去的。
方法二
这个做法思路更清晰解法更自然，我也不知道我为什么想不到，可能就是太菜了吧。
先介绍两个小技巧
子集反演
若
则
等等这不就是容斥原理吗 …… 这都有名字系列.jpg 。
子集求和
在已知 \(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,S\setminus\{i\}}\) 。滚动存储一下就是上面的代码。
回到这道题。仍然是预处理钦定第 \(i\) 种属性的某个数必须选的情况的答案。用 \(f_S\) 表示钦定第 \(i\) 种属性的某个数必须选时选出状态 \(S\) 的方案数。由于第 \(i\) 种属性已经钦定选了一个，所以只有当 \(S\) 是该属性的超集时才有意义。
从所有有意义的 \(S\) 中删去属性 \(i\) ，就变成了 \(i\) 的补集的全体子集，因此可以直接套用子集卷积的结论，即（设属性 \(i\) 为 \(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('\n');
}
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('\n');
}
return 0;
}
}
int main()
{
return zyt::work();
}
