题解：QOJ9619/洛谷13568 [CCPC 2024 重庆站] 乘积，欧拉函数，求和（数论+状压DP）

首先将 \(\phi(x)\) 拆成 \(\phi(x)= x \prod_{p | x} \frac {p-1}{p}\)，发现我们要求的式子其实可以转化为 \(\sum_{S} (\prod a_i)\prod_{p|\prod a_i} \frac {p-1}{p}\)。

发现其实我们只关心哪些质数 \(p\) 在最终的乘积里出现了，于是考虑记一个集合 \(S\) 表示出现了的质数。令 \(f_{i,S}\) 表示考虑前 \(i\) 个数，最后乘积出现的质数集合为 \(S\) 的 \(\prod a_i\) 之和，答案就是 \(\sum_S f_{n,S} \prod_{p\in S} \frac {p-1}{p}\)。转移的时候考虑 \(a_{i+1}\) 是否选入乘积：

1. 如果选择，则 \(f_{i,S} \times a_{i+1} \to f_{i+1,S\cup T}\)。
2. 如果不选，则 \(f_{i,S}\to f_{i+1,S}\)。

但是 \(3000\) 以内的素数有 \(430\) 个，我们显然不能全都记进 \(S\) 里面。分析一下，发现有很多状态是无效的：一个不超过 \(V\) 的整数，它最多只有一个大于 \(\sqrt V\) 的因子。所以，\(S\) 中大于 \(\sqrt V\) 的数最多只能有一个。

注意到小于 \(\sqrt V\) 的素数只有 \(16\) 个，所以我们不妨将所有数按照他的大于 \(\sqrt V\) 的因子（下方记为 \(F(x)\)）分组，记 \(f_{i,1/0,S}\) 表示考虑前 \(i\) 个数，当前这个数的 \(F(x)\) 有/没有被选进乘积，\(< \sqrt V\) 的素数选了集合 \(S\) 的 \(\prod a_i\) 之和。在转移的过程中，对于 \(F(x)\)，我们其实可以直接把贡献算进 \(f\) 里：

1. 如果 \(F(a_{i+1})=F(a_i)\)，即这两个是一组的，则：
   1. 如果 \(a_{i+1}\) 不选，\(f_{i,0/1,S}\to f_{i+1,0/1,S}\)。
   2. 如果 \(a_{i+1}\) 选，则考虑 \(a_i\) 是否选择：\((f_{i,0,S} \times \frac {F(a_{i+1})-1}{F(a_{i+1})}+f_{i,1,S}) \times a_{i+1} \to f_{i+1,1,S\cup T}\)。
2. 如果 \(F(a_{i+1})\ne F(a_i)\)，则：
   1. 如果 \(a_{i+1}\) 不选，\(f_{i,0,S}+f_{i,1,S}\to f_{i+1,0,S}\)。
   2. 如果 \(a_{i+1}\) 选，则 \((f_{i,0,S}+f_{i,1,S}) \times a_{i+1} \times \frac {F(a_{i+1})-1}{F(a_{i+1})} \to f_{i+1,1,S\cup T}\)。

再判一下如果 \(a_{i+1}\) 没有大于 \(\sqrt V\) 的因子，则按照上面最开始的方式转移即可。

最后时间复杂度 \(O(2^{16}n\times 414)\)。

const int MAXN = 2e3 + 5, MAXV = 3e3 + 5, mod = 998244353;
const int p[16] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53};
int n, f[2][2][1ll<< 17], inv[MAXV];
struct _node {
	int val, S, P;
} a[MAXN];

void add(int &x, int y) {
	x = (x + y) % mod;
}

int quickpow(int x, int y) {
	int ret = 1;
	while (y) {
		if (y & 1) ret = ret * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return ret;
}

void work() {
	cin >> n;
	for (int i = 1; i <= n; i++) {
		cin >> a[i].val;
		a[i].P = a[i].val;
		for (int j = 0; j < 16; j++) {
			if (a[i].P % p[j] == 0) {
				a[i].S |= (1 << j);
				while (a[i].P % p[j] == 0) a[i].P /= p[j];
			}
		}
		if (a[i].P > 1) inv[a[i].P] = quickpow(a[i].P, mod - 2);
		else a[i].P = 0;
	}
	for (int i = 0; i < 16; i++) {
		inv[p[i]] = quickpow(p[i], mod - 2);
	}
	sort(a + 1, a + 1 + n, [](auto x, auto y) {
		return x.P < y.P;
	});
	f[0][0][0] = 1;
	for (int i = 0; i < n; ++i) {
		memset(f[i & 1 ^ 1], 0, sizeof(f[i & 1 ^ 1]));
		for (int S = 0; S < (1ll << 16); ++S) {
			if (a[i + 1].P == 0) {
				add(f[i & 1 ^ 1][0][S], f[i & 1][0][S]);
				add(f[i & 1 ^ 1][0][S | a[i + 1].S], f[i & 1][0][S] * a[i + 1].val % mod);
				continue;
			}
			if (a[i + 1].P == a[i].P) {
				add(f[i & 1 ^ 1][0][S], f[i & 1][0][S]);
				add(f[i & 1 ^ 1][1][S], f[i & 1][1][S]);
				add(f[i & 1 ^ 1][1][S | a[i + 1].S], f[i & 1][0][S] * a[i + 1].val % mod * (a[i + 1].P - 1) % mod * inv[a[i + 1].P] % mod);
				add(f[i & 1 ^ 1][1][S | a[i + 1].S], f[i & 1][1][S] * a[i + 1].val % mod);
			} else {
				add(f[i & 1 ^ 1][0][S], f[i & 1][0][S]);
				add(f[i & 1 ^ 1][0][S], f[i & 1][1][S]);
				add(f[i & 1 ^ 1][1][S | a[i + 1].S], f[i & 1][0][S] * a[i + 1].val % mod * (a[i + 1].P - 1) % mod * inv[a[i + 1].P] % mod);
				add(f[i & 1 ^ 1][1][S | a[i + 1].S], f[i & 1][1][S] * a[i + 1].val % mod * (a[i + 1].P - 1) % mod * inv[a[i + 1].P] % mod);
			}
		}
	}
	int ans = 0;
	for (int S = 0; S < (1ll << 16); ++S) {
		int cof = 1; 
		for (int i = 0; i < 16; ++i) {
			if (S >> i & 1) cof = cof * (p[i] - 1) % mod * inv[p[i]] % mod;
		}
		add(ans, f[n & 1][0][S] * cof);
		add(ans, f[n & 1][1][S] * cof);
	}
	cout << ans << endl;
}