P10548 [THUPC 2024 决赛] 朔望

首先，我们可以直接认为每个行星的轨道是重合的，只是转一圈的时间不同。

我们考虑，两颗行星与中间的恒星共线的条件是什么，显然应该是行星之间隔了半圈或者在同一个位置。那么，如果认为一圈的长度是 $1$，任意两颗行星 $i>j$ 与恒星共线的周期是 $\frac{1}{2(\frac1{t_j}-\frac1{t_i})}=\frac{t_it_j}{2(t_i-t_j)}$（这就与小奥中的时钟问题是相似的）

我们需要钦定一个周期 $T$，$T$ 满足是任意 $t_i$ 的倍数。但是其实 $T$ 的具体值并不重要，因为它在后面的计算中必然会被抵消掉。

由于 $n\le20$，我们考虑去求出对于每个行星集合（至少要 $2$ 个行星）在 $T$ 时间内，集合内所有的行星都共线的次数。

现在我们有了两颗行星与恒星共线的条件。那么，怎么拓展到 $k$ 颗行星呢？这是很容易的，如果设行星们为 $p_1\lt p_2\lt\ldots\lt p_k$，那么我们只需判断对 $2\le i\le k$，是否有 $p_i$ 与 $p_1$ 与恒星共线。

即，对于一个时刻 $t$（不一定是整数），是否对 $\forall 2\le i\le k$，均有 $\frac{t}{\frac{t_{p_1}t_{p_i}}{2(t_{p_i}-t_{p_1})}}$ 是整数。

那么，周期是 $\operatorname{lcm}_{i=2}^k\frac{t_{p_1}t_{p_i}}{2(t_{p_i}-t_{p_i})}$，所以，这些行星在 $T$ 时刻内共线的次数是 $\gcd_{i=2}^k\frac{2(t_{p_i}-t_{p_1})T}{t_{p_i}t_{p_1}}$。

在小奥中，我们了解到两个分数 $\frac{p_1}{q_1}$，$\frac{p_2}{q_2}$ 的 $\gcd$ 是 $\frac{\gcd(p_1,p_2)}{\operatorname{lcm}(q_1,q_2)}$。这可能有点抽象，但是如果我们考虑将两个分数质因数分解，就会比较好理解。我们设 $\frac{p_1}{q_1}=\prod_{i=1}^kx_i^{a_i}$，$\frac{p_2}{q_2}=\prod_{i=1}^kx_i^{b_i}$。其中，$x_i$ 都是质数，并且 $a_i$，$b_i$ 都是整数（不一定为正）。那么

$$\gcd(\frac{p_1}{q_1},\frac{p_2}{q_2})=\prod_{i=1}^kx_i^{\min(a_i,b_i)}$$

这里我们可以带入整数或者上面给出的公式去计算，很容易发现是正确的。

所以，在上面给出的式子中，我们可以直接把 $2T$ 提出来，得到 $2T\gcd_{i=2}^k\frac{t_{p_i}-t_{p_1}}{t_{p_i}t_{p_1}}$。

注意，$\gcd$ 是不能取模的，所以我们只能统计出所有参与 $\gcd$ 的数的质因数分解，然后合并再取模。

于是，考虑预处理 $\frac{t_{p_i}-t_{p_1}}{t_{p_i}t_{p_1}}$ 的质因数分解。这个的时间复杂度大概是 $\mathcal O(n^2\sqrt V)$ 的。

然后对于每个集合，暴力求出在 $T$ 时间内共线的次数除以 $T$ 的值。这个地方可以做到 $\mathcal O(2^nn\omega(V))$ 左右。然后就可以求答案了？

并不是。因为，我们并没有限制集合之外的行星是否与集合中的行星共线。即，如果我们求出的答案是 $a_S$，真正的答案是 $b_S$，那么有

$$a_S=\sum_{S\subseteq T}b_T$$

反过来，

$$b_S=\sum_{T\subseteq S}a_T$$

注意到，枚举子集是 $\mathcal O(3^n)$ 的，不能通过。

那么，我们只能使用 IFWT 变换，具体来讲就是对于每一位做一遍容斥，代码是比较好写的。复杂度为 $\mathcal O(2^nn)$。

至此，思路内容结束。但是上面的内容其它的题解都讲得比较清楚了，接下来我主要讲实现的问题。

首先考虑时间复杂度的问题。这篇题解的思路瓶颈在于枚举子集暴力求出答案，这里需要精细实现，快速幂要预处理，并且注意不要再让复杂度升高。

代码有详细注释。

#include <bits/stdc++.h>
using namespace std;
const int mod = 1e9 + 7;
// 最好加上快速取模
inline void add(int &x, int y) {if ((x += y) >= mod) x -= mod;}
// 求 log(n)
inline int g(int x) {return 31 - __builtin_clz(x);}
// 最低位
inline int lowbit(int x) {return g(x & -x);}
// 这里记录的是对于所有预处理的值，质因数分解的结果
vector <int> v[25][25];
// 跟上面的一样，记录每个质数的计数
int num[25][25][10500];
// 记录的是一开始 ti - tj 每个质数的计数
int tb[25][25][10500];
// 求答案的时候计数每个质数的指数的最小值
int tmp[10500];
int t[25], w[25];
int p[10000], c;
bool ip[40000];
// 记录的是比 sqrt(V) 还大的质数
int nums[500];
int ans[1048576];
// 跟 tmp 一样的辅助数组，表示一个质数是否在所有行星中全部出现，是为了减少复杂度
int co[10500];
// 离散化专用
map <int, int> mp;
void sieve(int N) {
	for (int i = 2; i <= N; i++) {
		if (!ip[i]) p[++c] = i;
		for (int j = 1; j <= c && i * p[j] <= N; j++) {
			ip[i * p[j]] = 1;
			if (i % p[j] == 0) break;
		}
	}
}
// 预处理快速幂，不然过不去
int prepow[10500][2000];
inline int qpow(int x, int y, int res = 1) {
	if (y < 0) return qpow(qpow(x, mod - 2), -y);
	while (y) {
		if (y & 1) res = 1ll * res * x % mod;
		x = 1ll * x * x % mod; y >>= 1;
	}
	return res;
}
signed main() {
	int n, top = 0, cnt = 0, res = 0;
	scanf("%d", &n);
	for (int i = 1; i <= n; i++)
		scanf("%d", &t[i]);
	for (int i = 2; i <= n; i++)
		scanf("%d", &w[i]);
	sieve(sqrt(mod) + 5);
	for (int i = 1; i <= n; i++)
		for (int j = 0; j < i; j++) { // 表示处理 ti - tj 的质因数分解
			int x = t[i] - t[j];
			for (int k = 1; k <= c; k++)
				if (x % p[k] == 0)
					while (x % p[k] == 0)
						x /= p[k], ++num[i][j][k]; // 就是计数
			if (x > 1) nums[++cnt] = num[i][j][c + 1] = x;
		}
	sort(nums + 1, nums + cnt + 1);
	cnt = unique(nums + 1, nums + cnt + 1) - nums - 1;
	for (int i = 1; i <= cnt; i++) mp[nums[i]] = i;
	for (int i = 1; i <= cnt; i++) p[c + i] = nums[i]; // 最方便的写法，直接加到后面
	for (int i = 1; i <= n; i++)
		for (int j = 0; j < i; j++)
			if (num[i][j][c + 1]) { // 有一点点细节
				int t = mp[num[i][j][c + 1]];
				if (t == 1) num[i][j][c + 1] = 1;
				else num[i][j][c + 1] = 0, num[i][j][t + c] = 1;
			}
	for (int i = 2; i <= n; i++)
		for (int j = 1; j < i; j++)
			for (int k = 1; k <= c + cnt; k++)
				if (tb[i][j][k] = num[i][j][k] - num[i][0][k] - num[j][0][k]) // 压行
					v[i][j].push_back(k);
	for (int i = 1; i <= c + cnt; i++) {
		int pw = 1; // 这里我把下标加了 1000（有负数指数）
		prepow[i][1000] = 1;
		for (int j = 1; j <= 999; j++)
			prepow[i][j + 1000] = pw = 1ll * pw * p[i] % mod;
		pw = 1;
		int tmp = qpow(p[i], mod - 2);
		for (int j = -1; j >= -999; j--)
			prepow[i][j + 1000] = pw = 1ll * pw * tmp % mod;
	}
	for (int i = 1; i <= c + cnt; i++) tmp[i] = 1e9;
	for (int i = 3, j; i < 1 << n; i++) {
		if ((1 << (j = lowbit(i))) == i) continue; // 表示只有一颗行星
		int fac = 1, k = lowbit(i - (1 << j++)) + 1;
		for (int x = k; x <= n; x++)
			if (i & (1 << x - 1)) for (auto y : v[x][j])
				tmp[y] = min(tmp[y], tb[x][j][y]), ++co[y]; // 按照正常思路写就行了
		int cnt = __builtin_popcount(i);
		for (int x = k; x <= n; x++)
			if (i & (1 << x - 1)) for (auto y : v[x][j]) if (tmp[y] != 1e9) {
				if (tmp[y] > 0 && co[y] < cnt - 1) tmp[y] = 0; // 要特判，因为如果这个质数没有在所有行星中出现，那么它的指数最多是 0
				fac = 1ll * fac * prepow[y][tmp[y] + 1000] % mod;
				tmp[y] = 1e9, co[y] = 0; // 一定要清零，内存是重复利用的
			}
		ans[i] = fac;
	}
	for (int j = 1; j <= n; j++) // IFWT
		for (int i = (1 << n) - 1; ~i; i--)
			if (i & (1 << j - 1))
				add(ans[i ^ (1 << j - 1)], mod - ans[i]);
	for (int i = 3; i < (1 << n); i++)
		add(res, 1ll * ans[i] * w[__builtin_popcount(i)] % mod);
	return printf("%d", (res <<= 1) % mod), 0; // 我是最后再乘二的
}