Luogu P6624 「省选联考 2020 A 卷」作业题「矩阵树定理」

Luogu P6624 [省选联考 2020 A 卷] 作业题

题意：

定义一棵生成树的权值：

\[val(T) = \left(\sum_{i = 1}^{n-1} w_{e_i} \right)\times\gcd(w_{e_1},w_{e_2},\ldots,w_{e_{n-1}}) \]

求 \(G\) 所有生成树权值之和。对 \(998244353\) 取模。

题解：反演 + 矩阵树定理

答案就是 \(\sum g[k]\cdot k\)，其中 \(g[k]\) 表示边权 \(\gcd\) 恰好为 \(k\) 的生成树边权和之和。

通过莫比乌斯反演、欧拉反演或者容斥，只需求出 \(f[k]\) 表示边权 \(\gcd\) 是 \(k\) 倍数的生成树边权和之和，这个条件等价于生成树中所有边都是 \(k\) 倍数。

枚举 \(k\)，保留 \(k\mid w\) 的边，删去不连通的情况，只有 \(\mathcal O(n \cdot \max d(w_i))\) 种，其中 \(\max d(w_i) \leq 144\)，表示约数个数最大值。再写个矩阵树定理就是 \(\mathcal O(n^4 \cdot \max d(w_i))\)。

现转换为经典问题，求一张图所有生成树边权和之和。把一条边 \(w\) 变成一次生成函数 \(wx + 1\)（可以理解为两点间有 \(wx + 1\) 条边），对这个图求生成树个数，答案的一次项就是生成树边权和之和。因为对于一棵生成树，其贡献是 \(\prod (w_ix + 1)\)，直接展开一次项就是 \(\sum w_i\) 了。

我们要维护一个一次生成函数类，加减乘都非常基础，对于求逆，直接使用如下方式即可：

\[(a+wx)^{-1} \equiv \dfrac{1}{a^2}(a - wx) \pmod {x^2} \]

上面的算法足以通过 CCF 的数据，但是有点小问题：\(a = 0\) 时无法求逆，但可能有解。

我也不会处理，看 http://uoj.ac/submission/411018。

#include <bits/stdc++.h>
#define rep(i, j, k) for(int i = j; i <= k; ++ i)
#define per(i, j, k) for(int i = j; i >= k; -- i)
using namespace std;
typedef long long ll;
const int D = 152501 + 10, N = 32, M = N * N;
const int mod = 998244353;
inline int add(int x, int y) { return x + y >= mod ? x + y - mod : x + y; }
inline int dec(int x, int y) { return x - y < 0 ? x - y + mod : x - y; }
int n, m, d, f[D], u[M], v[M], w[M];
int qpow(int a, int b) {
   int ans = 1;
   for(; b >= 1; b >>= 1, a = (ll) a * a % mod)
      if(b & 1) ans = (ll) ans * a % mod;
   return ans;
}
struct F {
   int a, b; //a + bx
   void clr() { a = b = 0; }
   F operator + (F rhs) { return (F) {add(a, rhs.a), add(b, rhs.b)}; }
   F operator - (F rhs) { return (F) {dec(a, rhs.a), dec(b, rhs.b)}; }
   F operator * (F rhs) { return (F) {(int)(1ll * a * rhs.a % mod), (int)((1ll * a * rhs.b + 1ll * b * rhs.a) % mod)}; }
   void operator += (F rhs) { *this = *this + rhs; }
   void operator -= (F rhs) { *this = *this - rhs; }
   void operator *= (F rhs) { *this = *this * rhs; }
   F inv() {
      F ans = (F) {a, b ? mod - b : 0};
      if(a != 1) {
         int v = qpow(a, mod - 2);
         v = (ll) v * v % mod;
         ans.a = (ll) ans.a * v % mod;
         ans.b = (ll) ans.b * v % mod;
      }
      return ans;
   }
} mat[N][N];
struct _ufs {
   int f[N];
   void init(int n) { rep(i, 1, n) f[i] = i; }
   int find(int u) { return u == f[u] ? u : f[u] = find(f[u]); }
   void unite(int u, int v) { f[find(u)] = find(v); }
} ufs;
bool check(int d) {
   ufs.init(n);
   rep(i, 1, m) if(w[i] % d == 0) ufs.unite(u[i], v[i]);
   rep(i, 1, n) if(ufs.find(1) != ufs.find(i)) return false;
   return true;
}
int solve(int d) {
   rep(i, 1, n) rep(j, 1, n) mat[i][j].clr();
   rep(i, 1, m) if(w[i] % d == 0) {
      mat[u[i]][u[i]] += (F) {1, w[i]};
      mat[v[i]][v[i]] += (F) {1, w[i]};
      mat[u[i]][v[i]] -= (F) {1, w[i]};
      mat[v[i]][u[i]] -= (F) {1, w[i]};
   }
   F ans = (F) {1, 0};
   rep(i, 2, n) {
      int k = 0;
      rep(j, i, n) if(mat[j][i].a) { k = j; break; }
      if(!k) return 0;
      if(k != i) {
         // puts("swap");
         ans = (F) {0, 0};
         rep(j, i, n) swap(mat[i][j], mat[k][j]);  
      }
      F t = mat[i][i].inv(); ans *= mat[i][i];
      rep(j, i + 1, n) {
         F v = mat[j][i] * t;
         rep(k, i, n) mat[j][k] -= v * mat[i][k];
      }
   }
   return ans.b;
}
int main() {
   scanf("%d%d", &n, &m);
   rep(i, 1, m) { scanf("%d%d%d", u + i, v + i, w + i); d = max(d, w[i]); }
   rep(i, 1, d) f[i] = check(i) ? solve(i) : 0;
   per(i, d, 1) for(int j = i + i; j <= d; j += i) f[i] = dec(f[i], f[j]);
   int ans = 0;
   rep(i, 1, d) ans = (ans + (ll) f[i] * i) % mod;
   printf("%d\n", ans);
   return 0;
}