「SOL」Permanent (Codeforces)
这道题第一个结论都不知道怎么拿部分分啊
题意
一个 \(n\times n\) 的方阵 \(M\)，上面除了 \(k\) 个特殊位置，其他位置都是 \(1\)。第 \(i\) 个特殊位置在 \((x_i,y_i)\)，值为 \(v_i\)。
求 \(M\) 的积和式，即求对所有 \(1\sim n\) 的排列 \(P=(p_1,p_2,\dots,p_n)\) 下式的和：
其实就是行列式去掉 \(\pm 1\) 的系数。
数据规模：\(n\le10^6,k\le50\)。
解析
初步分析
积和式有一个比较易懂的公式，虽然不是很容易想到。记 \(\mathrm{perm}(M)\) 为 \(M\) 的积和式，方阵的加法定义为对应位置相加；对大小为 \(n\) 的方阵 \(M\) 以及全集 \(\{1,2,\dots,n\}\) 的子集 \(S,T\)，简记 \(M_{S,T}\) 表示按相对位置取出 \(M_{i,j}\)（\(i\in S,j\in T\)）的方阵，则
证明则考虑积和式的定义式：
相当于对于每个 \(i\) 要么选 \(A_{i,p_i}\) 要么选 \(B_{i,p_i}\)，枚举集合 \(S\)，对于 \(i\in S\) 选 \(A_{i,p_i}\)，\(i\in \bar S\) 选 \(B_{i,p_i}\)。则式子可以写为：
其中，\(S\) 对应的所有排列 \((p_i\mid i\in S)\) 就是把 \(S\) 中的元素拿来全部排列，\((p_i\mid i\in\bar S)\) 同理。于是可以得到上述结论，在该结论的基础上继续求解。
我们发现一个全 \(1\) 的方阵的积和式是易于计算的，即方阵边长的阶乘。于是考虑把原方阵 \(M\) 拆成全 \(1\) 方阵 \(A\) 和另一个方阵 \(B\) 的和，即
则我们要计算的就是
考虑积和式的实际意义——每行每列恰选一个元素的乘积。这种「每行每列恰选一个」可以想到行列拆点建二分图，则对于所有 \(|S|=i\)，\(\sum_{S}\sum_{T}^{|S|=|T|}\mathrm{perm}(B_{S,T})\) 就是在二分图上匹配 \(i\) 条边。
解法1
由于 \(B\) 上只有 \(k\) 个点非零，可以只考虑它们对应的行和列。
不难设计状压 DP，状压维护哪些行已经匹配过，枚举每一列再枚举一条边进行匹配，复杂度是 \(\mathcal O(2^kk)\)，显然过不了。
注意到 \(k\) 并不大，如果能把 \(2^k\) 变成 \(2^{\frac k2}\) 就能够优化很多。实际上，由于边数只有 \(k\)，所以连通块的大小至多 \(k+1\)；这就意味着每个连通块中，点数较少的一侧的点只有 \(\frac k2\) 个。
于是对每个连通块单独 DP，状压点数较少的一侧。复杂度 \(\mathcal O(2^{\frac k2}k)\)，更具体一些，设二分图点数为 \(n\)，则复杂度为 \(\mathcal O(2^{\frac n2}k)\)。
仍然不能通过所有数据，但在连通块个数较多（意味着点数与边数的差较大）时表现优秀。
解法2
仍然考虑匹配。
如果每个连通块都是一棵树，则可以对每个连通块 \(T\) 做一次 \(\mathcal O(|T|^2)\) 的树形背包。
显然大多数情况下并不是树。我们尝试先生成森林，再决策哪些非树边要选，然后做树形背包。
非树边的决策没什么好方法，只能指数级枚举。但是注意到非树边的数量并不大，只有 \(k-(n-1)\)（\(n\) 是总点数）。复杂度 \(\mathcal O(2^{k-n}n^2)\)。
这个算法在 \(k-n\) 较小时表现更优——也就是连通块个数较少，点数与边数更接近。
正解
两个算法综合，平衡复杂度。主要是平衡指数级枚举的复杂度——
如果 \(n\le \frac 23k\)，则采用算法1；否则采用算法2。最终复杂度 \(\mathcal O(2^{\frac k3}k^2)\)。
源代码
/* Lucky_Glass */
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define CON(typ) const typ &
const int M = 54, MOD = 1e9 + 7, N = 1e6 + 10;
inline int add(int a, CON(int) b) { return (a += b) >= MOD ? a - MOD : a; }
inline int sub(int a, CON(int) b) { return (a -= b) < 0 ? a + MOD : a; }
inline int mul(CON(int) a, CON(int) b) { return int(1ll * a * b % MOD); }
int pPow(int a, int b) {
int r = 1;
while (b) {
if (b & 1) r = mul(r, a);
a = mul(a, a), b >>= 1;
}
return a;
}
#define UPD(a, b, fun) a = fun(a, b)
int rin(int &r) {
int c = getchar(); r = 0;
while (c < '0' || '9' < c) c = getchar();
while ('0' <= c && c <= '9') r = r * 10 + (c ^ '0'), c = getchar();
return r;
}
struct Graph {
int head[M << 1], nxt[M << 1], to[M << 1], len[M << 1];
int cnt_edg;
inline void addEdge(CON(int) u, CON(int) v, CON(int) w) {
int p = ++cnt_edg, q = ++cnt_edg;
to[p] = v, nxt[p] = head[u], head[u] = p, len[p] = w;
to[q] = u, nxt[q] = head[v], head[v] = q, len[q] = w;
}
inline int operator [] (CON(int) u) const { return head[u]; }
};
Graph gr;
int n, m, cnt_x, cnt_y;
int id_x[N], id_y[N], edg[M][3], fac[N];
namespace DP_ON_BLOCK {
int cnt_e, cnt_f;
bool vis[M << 1];
int idx_e[M], idx_f[M], reid[M << 1], f[1 << 25], g[M];
void dfs(CON(int) u) {
vis[u] = true;
if (u <= cnt_x) idx_e[++cnt_e] = u;
else idx_f[++cnt_f] = u;
for (int it = gr[u]; it; it = gr.nxt[it])
if (!vis[gr.to[it]])
dfs(gr.to[it]);
}
/* cnt_a < cnt_b */
void doDP(CON(int) cnt_a, CON(int) cnt_b, int *idx_a, int *idx_b) {
for (int i = 1; i <= cnt_a; ++i) reid[idx_a[i]] = i - 1;
const int LIMITS = 1 << cnt_a;
for (int s = 0; s < LIMITS; ++s) f[s] = 0;
f[0] = 1;
for (int i = 1; i <= cnt_b; ++i) {
for (int s = LIMITS; ~s; --s)
for (int it = gr[idx_b[i]]; it; it = gr.nxt[it]) {
int v = reid[gr.to[it]];
if (!(s >> v & 1))
UPD(f[s | (1 << v)], mul(f[s], gr.len[it]), add);
}
}
int tmp_f[M] = {}, tmp_g[M] = {};
for (int s = 0; s < LIMITS; ++s) {
int pocnt = 0;
for (int i = 0; i < cnt_a; ++i) pocnt += s >> i & 1;
UPD(tmp_f[pocnt], f[s], add);
}
for (int i = 0; i <= cnt_a; ++i)
for (int j = 0; j <= m; ++j)
UPD(tmp_g[i + j], mul(tmp_f[i], g[j]), add);
for (int i = 0; i <= m; ++i) g[i] = tmp_g[i];
}
void main() {
for (int i = 1; i <= m; ++i)
gr.addEdge(edg[i][0], edg[i][1] + cnt_x, edg[i][2]);
g[0] = 1;
for (int i = 1; i <= cnt_x; ++i)
if (!vis[i]) {
cnt_e = cnt_f = 0;
dfs(i);
if (cnt_e < cnt_f) doDP(cnt_e, cnt_f, idx_e, idx_f);
else doDP(cnt_f, cnt_e, idx_f, idx_e);
}
int ans = 0;
for (int i = 0; i <= m; ++i) UPD(ans, mul(fac[n - i], g[i]), add);
printf("%d\n", ans);
}
} /* namespace DP_ON_BLOCK */
namespace DP_ON_TREE {
class Dsu {
private:
int fa[M << 1];
public:
void init(CON(int) siz) {
for (int i = 1; i <= siz; ++i)
fa[i] = i;
}
int find(CON(int) u) { return fa[u] == u ? u : fa[u] = find(fa[u]); }
inline bool merge(int u, int v) {
u = find(u), v = find(v);
if (u == v) return false;
fa[u] = v;
return true;
}
inline bool isRoot(CON(int) u) { return find(u) == u; }
};
int cnt_rem;
int rem_edg[M][3], f[M << 1][M][2], tmpf[M][2], g[M], tmpg[M];
bool isban[M << 1];
Dsu dsu;
int dfs(CON(int) u, CON(int) fa) {
int sizu = 1;
for (int i = 0; i <= m; ++i)
f[u][i][0] = f[u][i][1] = 0;
f[u][0][0] = 1;
for (int it = gr[u]; it; it = gr.nxt[it])
if (gr.to[it] != fa) {
int v = gr.to[it];
int sizv = dfs(v, u);
for (int i = 0, lmt_i = sizu >> 1; i <= lmt_i; ++i)
for (int j = 0, lmt_j = sizv >> 1; j <= lmt_j; ++j) {
UPD(tmpf[i + j][0],
mul(f[u][i][0], add(f[v][j][0], f[v][j][1])), add);
if (!isban[u]) {
UPD(tmpf[i + j][1],
mul(f[u][i][1], add(f[v][j][0], f[v][j][1])), add);
if (!isban[v])
UPD(tmpf[i + j + 1][1],
mul(mul(f[u][i][0], f[v][j][0]), gr.len[it]), add);
}
}
sizu += sizv;
for (int i = 0, lmt_i = sizu >> 1; i <= lmt_i; ++i) {
f[u][i][0] = tmpf[i][0], f[u][i][1] = tmpf[i][1];
tmpf[i][0] = tmpf[i][1] = 0;
}
}
return sizu;
}
void main() {
dsu.init(cnt_x + cnt_y);
for (int i = 1; i <= m; ++i)
if (dsu.merge(edg[i][0], edg[i][1] + cnt_x)) {
gr.addEdge(edg[i][0], edg[i][1] + cnt_x, edg[i][2]);
} else {
rem_edg[cnt_rem][0] = edg[i][0], rem_edg[cnt_rem][1] = edg[i][1] + cnt_x;
rem_edg[cnt_rem][2] = edg[i][2];
++cnt_rem;
}
int ans = 0;
for (int s = 0; s < (1 << cnt_rem); ++s) {
for (int i = 1; i <= cnt_x + cnt_y; ++i) isban[i] = false;
bool iscrash = false;
int pocnt = 0, tot_mul = 1;
for (int i = 0; i < cnt_rem; ++i)
if (s >> i & 1) {
++pocnt, UPD(tot_mul, rem_edg[i][2], mul);
if (isban[rem_edg[i][0]]) { iscrash = true; break; }
if (isban[rem_edg[i][1]]) { iscrash = true; break; }
isban[rem_edg[i][0]] = isban[rem_edg[i][1]] = true;
}
if (iscrash) continue;
for (int i = 0; i <= m; ++i) g[i] = 0;
g[0] = 1;
int cnt_g = 0;
for (int u = 1, lmt_u = cnt_x + cnt_y; u <= lmt_u; ++u)
if (dsu.isRoot(u)) {
int sizu = dfs(u, 0) >> 1;
for (int i = 0; i <= sizu; ++i)
for (int j = 0; j <= cnt_g; ++j)
UPD(tmpg[i + j], mul(add(f[u][i][0], f[u][i][1]), g[j]), add);
cnt_g += sizu;
for (int i = 0; i <= cnt_g; ++i)
g[i] = tmpg[i], tmpg[i] = 0;
}
for (int i = 0; i <= cnt_g; ++i)
UPD(ans, mul(fac[n - pocnt - i], mul(g[i], tot_mul)), add);
}
printf("%d\n", ans);
}
} /* namespace DP_ON_TREE */
void init() {
fac[0] = 1;
for (int i = 1; i < N; ++i) fac[i] = mul(fac[i - 1], i);
}
int main() {
/*
freopen("input.in", "r", stdin);
freopen("debug.out", "w", stdout);
*/
init();
rin(n), rin(m);
for (int i = 1, x, y, w; i <= m; ++i) {
rin(x), rin(y), rin(w);
if (!id_x[x]) id_x[x] = ++cnt_x;
if (!id_y[y]) id_y[y] = ++cnt_y;
edg[i][0] = id_x[x], edg[i][1] = id_y[y], edg[i][2] = sub(w, 1);
}
if (3 * (cnt_x + cnt_y) <= 2 * m) {
/* std::cerr << "block" << std::endl; */
DP_ON_BLOCK::main();
} else {
/* std::cerr << "tree" << std::endl; */
DP_ON_TREE::main();
}
return 0;
}
