高斯消元入门

问题

现在有一个 \(n\) 元一次方程组，类似下面这个样子：

\[\begin{cases} a_{1,1}x_1+a_{1,2}x_2+...+a_{1,n}x_n = b_1\\ a_{2,1}x_1+a_{2,2}x_2+...+a_{2,n}x_n = b_2\\ ...\\ a_{n,1}x_1+a_{n,2}x_2+...+a_{n,n}x_n = b_n \end{cases} \]

我们需要解出 \(x_1, x_2,...x_{n-1},x_n\)。可能有无解或者多组解。

算法流程和代码

首先用一个 \(n \times n\) 的矩阵描述系数，然后用一个 \(1\times n\) 的矩阵描述常数。 也就是下面这个样子：

\[\begin{bmatrix} a_{1,1} & \cdots & a_{1,n}\\ \vdots & \ddots & \vdots\\ a_{n,1} & \cdots & a_{n,n} \end{bmatrix} \begin{bmatrix} b_1\\ \vdots\\ b_n \end{bmatrix} \]

我们的目标是将第一个矩阵化为三角矩阵的形式，也就是类似下面这个样子：

\[\begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n-1} & a_{1,n}\\ 0 & a_{2,2} & \cdots & a_{2,n-1} & a_{2,n}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \dots & a_{n-1,n-1} & a_{n-1,n}\\ 0 & 0 & \dots & 0 & a_{n,n} \end{bmatrix} \begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_{n-1}\\ b_n \end{bmatrix} \]

形式化地，就是需要让第一个矩阵满足对于任何 \(n\ge i > j\ge 1\), \(A_{i,j} = 0\)。

得到这个矩阵之后，我们可以从下往上计算出所有 \(x\)。

举个简单的例子：

\[\begin{bmatrix} 1 & 3 & 2 & 1\\ 3 & -4 & 2 & 2\\ -2 & 3 & 1 & 4\\ 2 & 3 & 1 & 1 \end{bmatrix} \begin{bmatrix} 9\\ 18\\ -10\\ 3 \end{bmatrix} \]

首先使用第一行，也就是第一个方程消去其它行的 \(x_1\) 的系数，接下来方程组就变成了：

\[\begin{bmatrix} 1 & 3 & 2 & 1\\ 0 & -13 & -4 & 1\\ 0 & 9 & 5 & 6\\ 0 & -3 & -3 & -1 \end{bmatrix} \begin{bmatrix} 9\\ -9\\ 8\\ -15 \end{bmatrix} \]

接下来用第二行，消去第三第四行的 \(x_2\) 的系数：

\[\begin{bmatrix} 1 & 3 & 2 & 1\\ 0 & -13 & -4 & 1\\ 0 & 0 & \frac{29}{13} & \frac{69}{13}\\ 0 & 0 & -\frac{27}{13} & -\frac{10}{13} \end{bmatrix} \begin{bmatrix} 9\\ -9\\ \frac{23}{13}\\ -\frac{168}{13} \end{bmatrix} \]

最后用第三行消去第四行的 \(x_3\) 的系数：

\[\begin{bmatrix} 1 & 3 & 2 & 1\\ 0 & -13 & -4 & 1\\ 0 & 0 & \frac{29}{13} & \frac{69}{13}\\ 0 & 0 & 0 & \frac{121}{29} \end{bmatrix} \begin{bmatrix} 9\\ -9\\ \frac{23}{13}\\ -\frac{327}{29} \end{bmatrix} \]

这样就能解出 \(x_1...x_4\) 了。

但是仔细想一下，可以发现其实消元的时候可以不把矩阵消成三角矩阵的形式，而是直接消成对角矩阵的形式，也就是消成下面这个样子：

\[\begin{bmatrix} a_{1,1} & 0 & \cdots & 0 & 0\\ 0 & a_{2,2} & \cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & a_{n-1,n-1} & 0\\ 0 & 0 & \cdots & 0 & a_{n.n} \end{bmatrix} \begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_{n-1}\\ b_n \end{bmatrix} \]

形式化地，就是这个矩阵满足 \(\forall 1\le i,j\le n, i\ne j, A_{i,j} = 0\)，可以发现，对角矩阵实际上是三角矩阵的子集。有了对角矩阵之后，可以直接分别解出每个\(x\) 的值。

接着考虑无解或者多组解的情况。当我们在消元（注意，这时候构建的是对角矩阵）的时候，如果发现其中一个元 \(x_i\) ，其在所有剩余的方程中的系数都是 \(0\)，那么这个元不可解，直接略过它的消元过程。在消完所有可以解的元之后，在解 \(x\) 的值的时候，对于之前发现不可解的 \(x_i\)，如果它的方程的常数为 \(0\)，那么 \(x_i\) 有多个解（这时候 \(x_i\) 被称为自由元），否则方程组无解。

最后解决编程时候的精度问题，可以发现，当我们消 \(x_i\) 的系数的时候，我们其实可以选择任何没有被选择过的行保留 \(x_i\) 系数。而这一行的 \(x_i\) 系数要作为除数使用，这里可能会产生较大的精度损失。这里有两种方法：

1. 选择没被选择过的行中，\(x_i\) 系数较大的行保留 \(x_i\)系数，将其调换到第 \(i\) 行。这个方法最常用。
2. 选择第 \(i\) 行，第 \(i\) 列开始的右下角的矩阵中所有位置中系数最大的（假设为 \(a_{r,c}\) ） ，将 \(r\) 行和第 \(i\) 行调换，然后将第 \(c\) 列和第 \(i\) 列调换。会使代码长度增大。

现在给出代码，对应下面例题的模板2:

   int gauss(double a[maxn][maxn], int n, double* x) {
       int cnt = 1;
       for (int i = 1; i <= n; i++) {
           int line_id = cnt;
           for (int j = cnt + 1; j <= n; j++)
           if (fabs(a[line_id][i]) < fabs(a[j][i])) line_id = j;
           if (fabs(a[line_id][i]) < EPS) continue;
           for (int j = 1; j <= n + 1; j++) swap(a[cnt][j], a[line_id][j]);
           for (int j = 1; j <= n; j++) if (j != cnt) {
               double t = a[j][i] / a[cnt][i];
               for (int k = 1; k <= n + 1; k++) a[j][k] -= t * a[cnt][k];
           }
           cnt++;
       }
       if (cnt <= n) {
           while (cnt <= n) {
               if (fabs(a[cnt][n + 1]) > EPS) return -1; //no solution
               cnt++;
           }
           return -2; //multiple solutions
       }
       for (int i = 1; i <= n; i++) x[i] = a[i][n + 1] / a[i][i];
       return 0;
   }

模板

模板1(只需要判无解)，模板2(需要判多组解)

应用

计算有后效性的 dp 方程

1. [HNOI2013] 游走

2. [HNOI2011] XOR和路径

3. Broken robot(不是常规的高斯消元)

HNOI2013 游走

给定一个有 \(n\) 个节点， \(m\) 条无向边的无向连通图，保证任何两个能互相到达，没有重边和自环。小 Z 在该图上进行随机游走，初始时小 Z 在 \(1\) 号顶点，每一步小 Z 以相等的概率随机选择当前顶点的某条边，沿着这条边走到下一个顶点，获得等于这条边的编号的分数。当小 Z 到达 \(n\) 号顶点时游走结束，总分为所有获得的分数之和。现在需要给所有边都编一个序号（任意两条边的序号都不能相同，序号范围是 \(1\) 到 \(m\)），使得分数和的期望值最小。

\[1\leq n\leq 500, 1\leq m\leq 125000 \]

可以想到，根据期望的可加性，最优情况一定是按照边的经过次数的期望从小到大扫描边，从大到小分配编号。同时，注意到如果到达了 \(n\) ，就直接停止了，所以要关注特殊的情况。

设 \(f(u)\) 表示从节点 \(1\) 出发，经过节点 \(u\) 的次数的期望值，而 \(g(e)\)表示经过边 \(e\) 的次数的期望值。

首先想办法算出 \(f(u)\)：

\[f(u)=[u = 1] + \sum\limits_{(u,v)\in E} \frac{f(v)}{\mathrm{deg}(v)}, f(n) = 0 \]

稍微解释一下，\(\Sigma\) 中的分数表示从 \(v\) 到 \(u\) 的贡献，如果 \(u\) 就是起始点 \(1\)，那么根据我们的思路，需要额外加上 \(1\) 。这个式子显然是有后效性的，但因为 \(n\) 比较小，可以使用高斯消元求出 \(f\)。

接着解决 \(g(e)\) 的问题：

\[g[(u,v)] = \frac{f(u)}{\mathrm{deg}(u)}\times [u\ne n] + \frac{f(v)}{\mathrm{deg}(v)}\times [v\ne n] \]

经过边 \((u,v)\) 有两种方式，一种是经过 \(u\), 然后到 \(v\)，这要求 \(u\) 不能是 \(n\)，另一种是从 \(v\) 到 \(u\) ，同理。

完成上述计算之后，我们可以比较轻松地计算出最终的答案了。设 \(a\) 表示将 \(g\) 从小到大排序得到的长度为 \(m\) 的序列，那么最终答案为：

\[\mathrm{ANS} = \sum\limits_{i=1}^m (m-i+1)\times a_i \]

CF24D Broken robot

给定一个 \(n\times m\) 的网格图，有一个机器人在 \((x,y)\) 。现在可以等概率地选择向下，向左，向右，不动四个方向中可行的（就是不走出网格图的）方向构成的集合中选取一个方向，移动 \(1\) 个单位。问走到最后一行，也就是第 \(n\) 行的期望步数。

\[1\leq n, m\leq 1\times 10^3 \]

设定起始行为实际的第一行， 可以看出是一个期望 \(dp\)，依照反着来的套路，设 \(f(i,j)\) 为到达 \((i,j)\) 之后，走到最后一行所需要的期望步数，那么可以得到下面的转移方程组：

\[\begin{cases} f(n, i) &= 0\\ f(i,1) &= \frac{1}{3}\left[ f(i,1) + f(i,2) + f(i+1,1) \right] + 1\\ f(i,m) &= \frac{1}{3}\left[ f(i,m) + f(i,m-1) + f(i+1,m) \right] + 1\\ f(i,j) &= \frac{1}{4}\left[ f(i,j) + f(i,j+1) + f(i,j-1) + f(i+1,j) \right] + 1 \end{cases} \]

最终答案是 \(f(x,y)\)。

这个方程组在第一维没有后效性，但是在第二维却有后效性，所以不能单单用一个 for 循环解决。我们可以先将第 \(i+1\) 行的所有 \(dp\) 值求出来，然后想办法求第 \(i\) 行的 \(dp\) 值。

对方程进行变形，可以得到下面的三个方程：

\[\begin{cases} 2f(i,1) - f(i,2) &= f(i+1,1) + 3\\ -f(i,j - 1) + 3f(i,j) - f(i,j+1) &= f(i+1,j) + 4\\ -f(i,m-1) + 2f(i,m) &= f(i+1,m) + 3 \end{cases} \]

这样实际上就得到了 \(m\) 个方程。但是高斯消元是 \(O(n^3)\) 的，会 TLE。需要再次进行变形，把前 \(m - 1\) 个方程变成只有两个未知数，最有一个方程只有一个未知数的方程组（实际上就是手动构建三角矩阵）。

假设第 \(j-1(1 < j < n)\) 行变化之后的形式是 \(af(i,j-1) + bf(i,j) = c\)，可以得到：

\[f(i,j-1) = \frac{c - b \cdot f(i,j)}{a} \]

将其带第 \(j\) 行，可以得到：

\[\begin{aligned} \frac{b \cdot f(i,j) - c}{a} + 3f(i,j) - f(i, j+1) &= f(i+1,j) + 4\\ \left(3 + \frac{b}{a}\right) f(i,j) - f(i,j + 1) &= f(i+1,j) + 4 + \frac{c}{a} \end{aligned} \]

可以发现，无论是第几个方程，\(b = -1\)，那么再次化简，可以得到：

\[\left(3-\frac{1}{a}\right) f(i,j) - f(i,j+1) = f(i+1,j) + 4 + \frac{c}{a} \]

也就是说:

\[\begin{aligned} a' &= 3 - \frac{1}{a}\\ c' &= f(i+1,j) + 4 + \frac{c}{a} \end{aligned} \]

对于第 \(m\) 行，情况稍微特殊一点:

\[\begin{aligned} \frac{b \cdot f(i,m) - c}{a} + 2f(i,m) &= f(i+1,m) + 3\\ \left(2 - \frac{1}{a}\right) f(i,m) &= f(i+1, m) + 3 + \frac{c}{a} \end{aligned} \]

这样就能写出代码了，但是还有一种特殊情况没有讨论，就是当 \(m=1\) 的时候，转移方程是完全不一样的，这个时候：

\[f(i) = \frac{1}{2}[f(i) + f(i + 1)] + 1 \]

可以推出：

\[f(i) = f(i+1) + 2 \]

线性代数相关

线性相关的判定和线性基的计算

线性组合：一堆向量 \(B = (\vec v_1, \vec v_2,...\vec v_n)\) 再配合一堆系数 \((A_1, A_2,...A_n)\) 可以得到一个新的向量：

\[\vec v=\sum\limits_{i=1}^n A_i \vec v_i \]

我们称 \(\vec v\) 是 \((\vec v_1, \vec v_2,...\vec v_n)\) 的一个线性组合。

线性空间：\(S = (\vec v_1, \vec v_2,...\vec v_n)\) 组成的 \(\vec v\) 的集合 \(\mathrm{span}(S)\) 是 \(S\) 的线性空间。可以理解为我们常说的平面，直线，空间...

线性相关和线性无关： \(S\) 中的向量是线性相关的，当且仅当 \(S\) 中存在一个向量是其余向量的线性组合。如果不存在这个这种向量，就称 \(S\) 中的向量无关。可以理解一下，如果 \(S\) 中的向量是线性相关的，那么去掉其中一个向量，线性空间不会变。

线性基：\(S\) 的最大的线性无关的子集，这里称为 \(B\)，如果向量的维度数为 \(m\)，那么线性基的大小至多为 \(m\)。可以想到，\(B' \subset B, B' \neq B\) ，\(\mathrm{span}(B') \ne \mathrm{span}(B)\)。

当然，我们不可能枚举其中一个子集，然后判断是否线性无关。

考虑使用增量法，逐个加入向量，作为矩阵中新的一行，试构建三角矩阵（和之前一样是倒三角）。考虑一下这里的消元是在干什么。对于已经在线性基中的向量 \(\vec v_n\) ，它实际上是：

\[\vec{v_n} = \vec{v'_n} - \sum\limits_{i=1}^{n-1} a_i \vec{v_i} \]

其中 \(\vec{v'}_i\) 表示 \(v_i\) 还没有被消元时候的样子。

当我们用 \(\vec v_n\) 消去 \(\vec {nv}\) 的时候（这个时候已经用 \(\vec v_{1...n-1}\) 消了 \(\vec {nv}\)），实际上就是在调整 \(\vec {v'}_{1...n}\) 的参数 \(a\)。整个过程就是在尝试找一组解 \(a_{1...n}\) 使得:

\[\sum\limits_{i=1}^{n} a_i\vec{v_i} = \vec{nv} \]

如果 \(\vec{nv}\) 的被消成 \(\vec 0\) ，那么 \(\vec{nv}\) 就和之前的向量线性相关了，不需要加入（当然也可以将 \(nv\)
加入，然后取出之前向量中的其中一个），否则就将 \(\vec{nv}\) 加入。

看一道模板题：[JLOI2015]装备购买

给定 \(n\) 个向量，向量的维度数为 \(m\)，每个向量有一个选择的代价，要选择一个代价和最小的线性基，输出线性基大小和最小代价。

可以想到，像 MST 一样，将向量按从小到大尝试加入，如果当前向量和之前的向量线性相关，那么就不能加入，否则加入。

  int n, m, ans1, ans2, ls[maxn];
  struct Note {
      D v[maxn]; int c;
      friend bool operator < (Note a, Note b) {
          return a.c < b.c;
      }
  } a[maxn];
  sort(a + 1, a + 1 + n);
  for (int i = 1; i <= n; i++) {
      for (int j = 1; j <= m; j++) {
          if (fabs(a[i].v[j]) < EPS) continue;
          if (!ls[j]) {
              ls[j] = i, ans1++, ans2 += a[i].c;
              break;
          }
          D t = a[i].v[j] / a[ls[j]].v[j];
          for (int k = 1; k <= m; k++) a[i].v[k] -= t * a[ls[j]].v[k];
      }
  }

行列式计算

学过线性代数的 dalao 应该知道，对于一个 \(n\times n\) 的矩阵，其行列式可以表示为：

\[\det (A) = |A| = \sum\limits_{P} (-1) ^ {\tau(P)} \prod\limits_{i=1}^n a_{i,p_i} \]

其中 \(P\) 是 \(1\) 到 \(n\) 的排列， \(\tau(A)\) 表示序列 \(A\) 的逆序对个数。

这个东西有一些不错的性质：

1. 交换矩阵的两行，其行列式值取反。这个比较显然。

2. 当矩阵中存在一行完全为 \(0\) ，那么矩阵的行列式为 \(0\)。（废话）

3. 某一行的值都乘上常数 \(a\)，矩阵的行列式值乘上 \(a\)。

4. 用公式表示如下：

   \[\begin{vmatrix}a+a' & b+b'\\c & d\end{vmatrix} = \begin{vmatrix}a & b \\ c & d\end{vmatrix} + \begin{vmatrix} a' & b'\\ c & d \end{vmatrix} \]

   这个用分配律可以证明，就不说了。

5. 若矩阵中的每一行都相同，那么矩阵的行列式为 \(0\) 。

6. 矩阵中的一行减去另一行的常数倍，矩阵的行列式不变。形式化表示如下：

   \[\begin{vmatrix} a & b \\ c - ka & d - kb \end{vmatrix} = \begin{vmatrix} a & b \\ c & d \end{vmatrix} \]

   证明：首先可以得到下面的推论：

   \[\begin{vmatrix} a & b\\ka & kb \end{vmatrix} = \frac{1}{k}\cdot \begin{vmatrix} ka & ka \\ ka & kb \end{vmatrix} = 0 \]

   接着证明原式：

   \[\mathrm{Left}= \begin{vmatrix} a & b\\c & d \end{vmatrix} - \begin{vmatrix} a & b\\ka & kb \end{vmatrix} = \mathrm{Right} \]

7. 三角矩阵的行列式为对角线元素之和。证明也比较简单。

联合上面的所有性质，可以发现，我们可以用高斯消元的方法，将整个矩阵化成三角矩阵的形式，将最终得到的矩阵的行列式乘上消元过程产生的系数得到原矩阵的行列式。

这样就可以写出常规的求行列式代码了。

 double det(double a[maxn][maxn], int n) {
     int fl = 1; double dt = 1;
     for (int i = 1; i <= n; i++) {
         int line_id = i;
         for (int j = i + 1; j <= n; j++)
         if (fabs(a[line_id][i]) < fabs(a[j][i])) line_id = j;
         if (line_id != i) {
             for (int j = 1; j <= n; j++) swap(a[line_id][j], a[i][j]);
             fl = -fl;
         }
         if (fabs(a[i][i]) < EPS) return 0;
         for (int j = i + 1; j <= n; j++) {
             double t = a[j][i] / a[i][i];
             for (int k = 1; k <= n; k++) a[j][k] -= t * a[i][k];
         }
     }
     for (int i = 1; i <= n; i++) dt *= a[i][i];
     return dt * fl;
 }

但是这个代码的精度并不高，而且有的时候我们还需要对计算结果进行取模，那么上面的代码就不能用了，但是我们并不喜欢使用逆元，因此我们需要不使用实数除法的消元方法：辗转消元法。该方法和求最大公约数用的辗转相除法过程有共通之处，在此直接给出代码。

 #define LL long long
  LL det(LL a[maxn][maxn], int n, LL MOD) {
      int fl = 0; LL dt = 1;
      for (int i = 1; i <= n; i++) {
          int line_id = -1;
          for (int j = i; j <= n; j++) if (a[j][i] && (line_id == -1 || a[j][i] > a[line_id][i])) { 
                  line_id = j; break; 
              }
          if (line_id == -1) return 0;
          if (i != line_id) swap(a[line_id], a[i]), fl ^= 1;
          for (int j = i + 1; j <= n; j++) {
              if (a[j][i] > a[i][i]) swap(a[i], a[j]), fl ^= 1;
              while (a[j][i]) {
                  int t = a[i][i] / a[j][i];
                  for (int k = i; k <= n; k++) a[i][k] = (a[i][k] - t * a[j][k] % MOD + MOD) % MOD;
                  swap(a[i], a[j]), fl ^= 1;
              }
          }
      }
      for (int i = 1; i <= n; i++) dt = dt * a[i][i] % MOD;
      if (fl) return (MOD - dt) % MOD;
      return dt;
  }

矩阵树定理

矩阵树定理可以用于计算生成树个数，权值积。在这里只介绍无向图生成树的版本。

生成树个数

定义一个度数矩阵 \(D\)，其构造方式为 \(D_{i,i} = \mathrm{deg}(i), D_{i,j}=0 (i \ne j)\)。定义一个邻接矩阵\(A\)，满足 \(\forall (u,v)\in E, A_{u,v} = A_{v,u} = c(u,v)\), \(c(u,v)\) 表示连接 \(u,v\) 的边数。最后定义 Kirchhoff（基尔霍夫）矩阵 \(K = D - A\)。

对于任何 \(i \in [1,n]\)，都有$ \mathrm{ANS}=|C_{i}[K]|$, 其中 \(C_i[K]\) 表示矩阵 \(K\) 去除第 \(i\) 行，第 \(i\) 列之后的矩阵，一般来说都会去除最后一行，最后一列。

生成树边权积

这个才是矩阵数的本质在计算的东西：

\[\mathrm{ANS} = \sum\limits_{T\in TS} \prod_{e\in T} w_e \]

其中 \(TS\) 是原图的生成树的集合。其中 \(\prod_{e\in T} w_e\) 可以认为是 \(w_e\) 的卷积，而实数的卷积就是乘法。

更改一下矩阵的定义就行了。

\[\begin{aligned} D_{u,u} &= \sum\limits_{(u,v) \in E} w(u,v)\\ A_{u,v} &= w(u,v)\\ K &= D - A\\ \mathrm{ANS} &= |C_i[K]| (1\leq i\leq n) \end{aligned} \]

其中 \(w(u,v)\) 表示所有 \(u,v\) 之间的边的权值和（重边的权值和）。

生成树边权和

对于原图中边权为 \(w\) 的边，将其边权更改为多项式 \(1+wx\)。那么
\(\mathrm{ANS}\) 的一次项系数就是边权和。

可以发现，在计算的时候，结构体只需要记录下一次项系数和常数就好了。多项式除法可以使用下面的式子：

\[\frac{ax+b}{cx+d}=\frac{ad-bc}{d^2} x+\frac{b}{d} \]

例题列表

1. [JSOI2008]最小生成树计数

2. [省选联考 2020 A 卷]作业题

3. [北京省选集训2019]生成树计数

[JSOI2008]最小生成树计数 (PLUS)

给定一个有 \(n\) 个点，\(m\) 条带权边的无向图，相同边权的边的数量不超过 \(50\) 条（原题是 \(10\) 条），求出有多少颗不同的最小生成树，两个生成树 \(T_1, T_2\) 不相同，当且仅当 \(\exists e\in T_1, e\notin T_2\)。

\[1\leq n\leq 100, 1\leq m\leq 1000 \]

MST 有两个比较重要的结论：

1. 对于任意边权, 在两颗不同的最小生成树中的数量相同。

2. 不同的生成树中，某一种权值的边任意加入需要的数量后，形成的联通块状态是一样的。

可以发现，每一种边权对于答案的贡献可以分开计算，然后乘起来。

首先构建一颗最小生成树，记录一下使用的边。 对于某一边权 \(w\) 的答案统计，可以首先将先前得到的最小生成树中除边权为 \(w\) 以外的所有边加入，将连通块缩点（并查集）。在缩点之后的图上，用原图中所有的边权为 \(w\) 的边跑矩阵树，得到的答案就是边权 \(w\) 的贡献。

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

给定一颗 \(n\) 个点，\(m\) 条带权边的无向图，计算下面的式子，对 \(998244353\) 取模：

\[\sum\limits_{T\in TS} \left(\sum\limits_{e\in T} w_e\right) \gcd\left[w_e(e\in T)\right] \]

\[1\leq n\leq 30 \]

化式子：

\[=\sum\limits_{d=1}^n \varphi(d) \sum\limits_{\forall e\in T, d | w_e} \sum\limits_{e\in T} w_e \]

可以发现，可以使用线性筛求出 \(\varphi\)，然后对于枚举的 \(d\)，可以使用矩阵数求一下边权和。

[北京省选集训2019]生成树计数

给定一个 \(n\) 个点的带权无向完全图，求其所有生成树权值的 \(k\)
次方之和，答案对 \(998244353\) 取模。形式化地，就是求出下面这个式子：

\[\sum\limits_{T\in TS}\left(\sum\limits_{e \in T} w_e\right) ^ k \]

\[1\leq n, k\leq 30 \]

做这道题之前，我们必须知道 \(e^x\) 的泰勒展开，也就是它的生成函数：

\[e^x = \sum\limits_{i=0}^{\infty} \frac{x^i}{i!} \]

接着要知道，指数型生成函数的积就是：

\[\begin{aligned} F(x)G(x) & = \sum\limits_{i=0}^{\infty}\sum\limits_{j=0}^{\infty} \frac{a_i x^i}{i!}\frac{b_j x_j}{j!}\\ & = \sum\limits_{i=0}^{\infty}\sum\limits_{j=0}^i \binom{i}{j}a_jb_{i-j} \frac{x_i}{i!} \end{aligned} \]

和就是说，新的生成函数的第 \(i\) 项（从零开始标号），就是 \(\sum_{j=0}^i \binom{i}{j}a_jb_{i-j}\)，长得和二项式定理一样。考虑可合并的方式计算和的 \(k\) 次方，也可以想到使用二项式定理，也就是这个样子: (\(W\)
为先前的边权和， \(w\) 为新加入的边权)

\[(W+w)^k = \sum\limits_{i=0}^k \binom{k}{i} W^i w^{k-i} \]

一次边权的合并，可以对应一次生成函数的积。

对于一条原图中边权为 \(w\) 的边，将其边权改为 \(e^{wx}\)，它的生成函数是这样的：

\[e^{wx} = \sum\limits_{i=0}^{\infty} w^i\frac{x^i}{i!} \]

另外，我们需要求出一个生成函数的逆元的前 \(k\) 项，也就是求出\(F^{-1}(x)\):

\[\begin{aligned} F^{-1}(x)F(x)&=1\\ \sum\limits_{i=0}^{\infty}\sum\limits_{j=0}^i \binom{i}{j}a'_ja_{i-j} \frac{x_i}{i!} &= 1\\ \end{aligned} \]

可以令 \(a'_0a_0=1\)，剩余的(\(i>0\)) 满足 \(\sum_{j=0}^i \binom{i}{j}a'_ja_{i-j} = 0\)。可以推出式子：

\[a'_i=-\frac{1}{a_0}\sum\limits_{j=0}^{i-1} \binom{i}{j} a'_j a_{i-j} \]

跑一遍矩阵树，答案就是生成函数的第 \(k\) 项。