这道题第一个结论都不知道怎么拿部分分啊
题意
一个 (n imes n) 的方阵 (M)，上面除了 (k) 个特殊位置，其他位置都是 (1)。第 (i) 个特殊位置在 ((x_i,y_i))，值为 (v_i)。
求 (M) 的积和式，即求对所有 (1sim n) 的排列 (P=(p_1,p_2,dots,p_n)) 下式的和：
[prod_{i=1}^nM_{i,p_i}
]
其实就是行列式去掉 (pm 1) 的系数。
数据规模：(nle10^6,kle50)。
解析
初步分析
积和式有一个比较易懂的公式，虽然不是很容易想到。记 (mathrm{perm}(M)) 为 (M) 的积和式，方阵的加法定义为对应位置相加；对大小为 (n) 的方阵 (M) 以及全集 ({1,2,dots,n}) 的子集 (S,T)，简记 (M_{S,T}) 表示按相对位置取出 (M_{i,j})（(iin S,jin T)）的方阵，则
[mathrm{perm}(A+B)=sum_{Ssubseteq{1,2,dots,n}}sum_{Tsubseteq{1,2,dots,n}}^{|S|=|T|}mathrm{perm}(A_{S,T})cdotmathrm{perm}(B_{ar S,ar T})
]
证明则考虑积和式的定义式：
[mathrm{perm}(A+B)=sum_Pprod_{i=1}^n(A_{i,p_i}+B_{i,p_i})
]
相当于对于每个 (i) 要么选 (A_{i,p_i}) 要么选 (B_{i,p_i})，枚举集合 (S)，对于 (iin S) 选 (A_{i,p_i})，(iin ar S) 选 (B_{i,p_i})。则式子可以写为：
[mathrm{perm}(A+B)=sum_{P}sum_{S}prod_{iin S}A_{i,p_i}prod_{jin ar S}B_{i,p_i}
]
其中，(S) 对应的所有排列 ((p_imid iin S)) 就是把 (S) 中的元素拿来全部排列，((p_imid iinar S)) 同理。于是可以得到上述结论，在该结论的基础上继续求解。
我们发现一个全 (1) 的方阵的积和式是易于计算的，即方阵边长的阶乘。于是考虑把原方阵 (M) 拆成全 (1) 方阵 (A) 和另一个方阵 (B) 的和，即
[B_{i,j}=egin{cases}
v_k-1&i=x_k,j=y_k\
0& ext{otherwise}
end{cases}
]
则我们要计算的就是
[sum_{S}sum_{T}^{|S|=|T|}mathrm{perm}(B_{S,T})cdot(|ar S|!)
]
考虑积和式的实际意义——每行每列恰选一个元素的乘积。这种「每行每列恰选一个」可以想到行列拆点建二分图，则对于所有 (|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) 较小时表现更优——也就是连通块个数较少，点数与边数更接近。
正解
两个算法综合，平衡复杂度。主要是平衡指数级枚举的复杂度——
如果 (nle 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
", 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
", 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;
}