好厉害。
首先发现一件事情，就是假如存在一个 \(p | i\)，那么假如某次操作为 \(i | j, i| k\)，那么我们选 \(p|j, p|k\) 是不会更劣的。那么这就意味着，一定存在一些数 \(p\) 满足没有数是它的因数，而这些数很显然也是不可能被删除的。我们找出这样的数的集合 \(T\)，剩下的数的集合 \(S\)，那么容易发现，任何一个操作都可以变成 \(T\) 中的一个数与 \(S\) 中的两个数的操作。我们将倍数关系进行连边，这样我们得到了一张二分图，而题目中的操作相当于每次用 \(T\) 中选一个度数大于等于 \(2\) 的点，在它相邻的所有点中选一个删掉。
正着考虑不好考虑，我们可以反过来考虑，假如最后剩下了某个集合，那么每次操作就相当于将一个与 \(T\) 相连的度数大于等于 \(1\) 的点加入。首先我们把每个弱连通块分开考虑，最后用组合数将每个连通块中的答案合并起来。那么发现，对于一个连通块，只要最后剩下了一个点，就一定能将所有点全部加入，所以最后能剩下的最少的数一定是 \(T\) 的大小加 \(1\)。
那么我们现在就有一个朴素的做法，考虑大力状压，设 \(f_{S}\) 为已经加入 \(S\) 集合中的点的方案数。问题很显然，\(n \le 60\)，且 \(S\) 集合很容易卡满，状态数太多，无法接受。
此时我们考虑的是 \(S\) 集合，假如从 \(T\) 集合的角度来考虑呢？我们记录当前集合中相邻的点的集合为 \(T\)，且当前集合大小为 \(i\)，设 \(f_{T, i}\) 为这种情况时的方案数。转移时考虑是选一个完全包含在 \(T\) 中的点或者选一个能够拓展 \(T\) 集合的点，容易转移。但是这样复杂度看起来更劣了啊？
分析 \(T\) 集合的大小，首先由于 \(T\) 集合一定要有一个倍数，那么首先需要满足这个集合中的数 \(\le 30\)。继续考虑，\(T\) 集合相当于是整个偏序关系的一个最长反链，由 Dilworth 定理，最长反链等于最小链覆盖。我们容易构造 \(\frac{30}{2}\) 条链（考虑对每个奇数 \(k\) 构造一条 \(k \times 2^i\) 的链，类似于 ARC141D），那么我们就能证明最长反链不超过 \(\frac{30}{2} = 15\)，即 \(T\) 集合的大小不超过 \(15\)。那么这个复杂度就能接受了，直接状压即可。
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 66, P = 1000000007;
int n, m, k;
int a[MAXN];
bool in[MAXN];
int id[MAXN];
int to[MAXN];
int f[1 << 15][MAXN];
int cnt[1 << 15];
vector<int> e[MAXN];
bool vis[MAXN];
vector<int> pt;
void dfs(int u) {
vis[u] = 1;
m++, k += in[u];
if (in[u]) id[u] = k - 1;
pt.push_back(u);
for (int v : e[u]) if (!vis[v]) {
dfs(v);
}
}
int C[MAXN][MAXN];
int main() {
scanf("%d", &n);
C[0][0] = 1;
for (int i = 1; i <= n; i++) {
C[i][0] = 1;
for (int j = 1; j <= i; j++) {
C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % P;
}
}
for (int i = 1; i <= n; i++) {
scanf("%d", &a[i]);
}
for (int i = 1; i <= n; i++) {
bool flag = true;
for (int j = 1; j <= n; j++) if (j != i) {
if (a[i] % a[j] == 0) {
flag = false;
break;
}
}
if (flag) {
in[i] = 1;
}
}
for (int i = 1; i <= n; i++) if (!in[i]) {
for (int j = 1; j <= n; j++) if (i != j && in[j]) {
if (a[i] % a[j] == 0) {
e[i].push_back(j), e[j].push_back(i);
}
}
}
int tot = 0, ans = 1;
for (int rt = 1; rt <= n; rt++) if (!vis[rt]) {
m = k = 0;
pt.clear();
dfs(rt);
if (m == k) continue;
for (int s = 0; s < (1 << k); s++) {
cnt[s] = 0;
for (int i = 0; i <= m - k; i++) f[s][i] = 0;
}
for (int i : pt) if (!in[i]) {
for (int j : pt) if (i != j && in[j]) {
if (a[i] % a[j] == 0) {
to[i] |= 1 << id[j];
}
}
cnt[to[i]]++;
}
for (int mid = 1; mid < (1 << k); mid <<= 1) {
for (int l = 0; l < (1 << k); l += (mid << 1)) {
for (int i = 0; i < mid; i++) {
cnt[l + i + mid] += cnt[l + i];
}
}
}
for (int i : pt) if (!in[i]) {
f[to[i]][1]++;
}
for (int s = 1; s < (1 << k); s++) {
for (int i = 1; i <= m - k; i++) if (f[s][i]) {
f[s][i + 1] = (f[s][i + 1] + 1ll * f[s][i] * (cnt[s] - i)) % P;
for (int j : pt) if (!in[j] && (to[j] & s) != 0 && ((to[j] | s) != s)) {
f[s | to[j]][i + 1] = (f[s | to[j]][i + 1] + f[s][i]) % P;
}
}
}
ans = 1ll * ans * f[(1 << k) - 1][m - k] % P * C[tot + m - k - 1][m - k - 1] % P;
tot += m - k - 1;
}
printf("%d\n", ans);
return 0;
}