首先将 \(\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}\) 是否选入乘积：
- 如果选择，则 \(f_{i,S} \times a_{i+1} \to f_{i+1,S\cup T}\)。
- 如果不选，则 \(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\) 里：
- 如果 \(F(a_{i+1})=F(a_i)\)，即这两个是一组的，则：
- 如果 \(a_{i+1}\) 不选，\(f_{i,0/1,S}\to f_{i+1,0/1,S}\)。
- 如果 \(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}\)。
- 如果 \(F(a_{i+1})\ne F(a_i)\)，则：
- 如果 \(a_{i+1}\) 不选，\(f_{i,0,S}+f_{i,1,S}\to f_{i+1,0,S}\)。
- 如果 \(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;
}
