「算法笔记」线性代数基础

逐渐修改于后续不知道哪年不知道哪月。

2020 年写的（已折叠）

注：这是我什么都不会的时候写的（东抄西抄拼起来），有很多锅，建议不要看了 QAQ。

一、行列式

1. 定义

二阶行列式：\( \begin{vmatrix} a_{11} & a_{12}\\ a_{21} & a_{22} \end{vmatrix} =a_{11}a_{22}-a_{21}a_{12} \)。即主对角线的乘积减去副对角线的乘积。

三阶行列式：三阶行列式的计算方法与二阶行列式类似。

\( \begin{vmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{vmatrix}=a_{11}a_{22}a_{33}+a_{12}a_{23}a_{31}+a_{13}a_{21}a_{32}-a_{13}a_{22}a_{31}-a_{12}a_{21}a_{33}-a_{11}a_{23}a_{32} \)

\(n\) 阶行列式：

由 \(n^2\) 个元素构成的 \(n\) 阶行列式为：

\( \left| a_{ij}\right|_n= \begin{vmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{vmatrix} =\sum\limits_{j_1 j_2 \cdots j_n}(-1)^{r(j_1 j_2 \cdots j_n)}a_{1j_1}a_{2j_2}\cdots a_{nj_n} \)

其中 \(j_1 j_2 \cdots j_n\) 是 \(1,2,\cdots,n\) 的一个排列。

\(r(j_1 j_2 \cdots j_n)\) 是排列的逆序对数量。

2. 性质

设 \(n\) 阶行列式为 \(D=\left| a_{ij}\right| _n\)，称行列式 \(\left| a_{ji}\right| _n\)（行与列交换一下，即 \(a_{ij}\) 变为 \(a_{ji}\)）为 \(D\) 的转置行列式，记作 \(D^T\)。

一些性质：

• \(D^T=D\)。

• 交换行列式 \(D\) 中任意两行得到 \(D_1\)，则 \(D_1=-D\)。

• 行列式 \(D\) 中某一行都乘上 \(k\) 得到 \(D_1\)，则 \(D_1=kD\)。

• 行列式 \(D\) 中的某一行加上另一行的 \(c\) 倍（对应相加，其中 \(c\) 为一个常数）得到 \(D_1\)，则 \(D_1=D\)。

3. 代数余子式

余子式与代数余子式：

在 \(n\)（\(n>1\)） 阶行列式 \(D=|a_{ij}|\) 中，划去元素 \(a_{ij}\) 所在的第 \(i\) 行和第 \(j\) 列，余下的元素按原来顺序组成的 \(n-1\) 阶行列式，称为 \(D\) 中元素 \(a_{ij}\) 的余子式，即为 \(M_{ij}\)。

\(a_{ij}\) 的余子式 \(M_{ij}\)，在它前面添加符号 \((-1)^{i+j}\) 后，称为 \(a_{ij}\) 在 \(D\) 中的代数余子式，记为 \(A_{ij}\)，即 \(A_{ij}=(-1)^{i+j} M_{ij}\)。

\(n\) 阶行列式 \(D=|a_{ij}|\) 等于它的任意一行（列）中各元素与其对应的代数余子式乘积的和，即 \(D=a_{i1}A_{i1}+a_{i2}A_{i2}+\cdots+a_{in}A_{in}\) 或 \(D=a_{1j}A_{1j}+a_{2j}A_{2j}+\cdots+a_{nj}A_{nj}\)（\(i,j=1,2,\cdots,n\)）。

举个栗子：

\(D=\begin{vmatrix}1 & 2 & 3\\0 & 2 & 0\\2 & 3 & 1\end{vmatrix}=0\cdot A_{21}+2\cdot A_{22}+0\cdot A_{23}\)
\(=2\cdot A_{22}=2\cdot (-1)^{2+2}\begin{vmatrix}1 & 3\\2 & 1\end{vmatrix}=2\times (-5)=-10\)

 \(k\) 阶子式的余子式与代数余子式：

在 \(n\) 阶行列式 \(D\) 中划去任意选定的 \(k\) 行、\(k\) 列后（\(0<k<n\)），位于这些行和列交叉处的 \(k^2\) 个元素按原来顺序组成的 \(k\) 阶行列式 \(M\)，称为 \(M\) 是行列式 \(D\) 的一个 \(k\) 阶子式。

在 \(n\) 阶行列式 \(D\) 中划去任意选定的 \(k\) 行、\(k\) 列后（\(0<k<n\)），余下的元素按原来顺序组成的 \(n-k\) 阶行列式 \(N\)，称为 \(N\) 是 \(k\) 阶子式 \(M\) 的余子式。

如果 \(k\) 阶子式 \(M\) 在行列式 \(D\) 中的行和列的标号分别为 \(i_1,i_2,\cdots,i_k\) 和 \(j_1,j_2,\cdots,j_k\)，则在 \(M\) 的余子式 \(N\) 前面添加符号 \((-1)^{(i_1+i_2+\cdots+i_k)+(j_1+j_2+\cdots+j_k)}\) 后，所得到的 \(n-k\) 阶行列式，称为行列式 \(D\) 的 \(k\) 阶子式 \(M\) 的代数余子式，记作 \(A\)。即 \(A=(-1)^{(i_1+i_2+\cdots+i_k)+(j_1+j_2+\cdots+j_k)}N\)。

（主子式：去掉的行和列编号相同。）

4. 范德蒙德行列式

\(D_n=\begin{vmatrix} 1 & 1 & \cdots & 1 \\ x_1 & x_2 & \cdots & x_n \\ x_1^2 & x_2^2 & \cdots & x_n^2\\ \vdots & \vdots & \ddots & \vdots \\ x_1^{n-1} & x_2^{n-1} & \cdots & x_n^{n-1} \end{vmatrix} =\prod\limits_{1\leq i<j\leq n}(x_j-x_i)\)

用数学归纳法证明。

当 \(n=2\) 时，\(D_2=\begin{vmatrix} 1 & 1 \\ x_1 & x_2\end{vmatrix}=x_2-x_1\)，显然结论成立。

假设该结论对 \(n-1\) 阶范德蒙行列式成立，即：

\(D_{n-1}=\begin{vmatrix} 1 & 1 & \cdots & 1 \\ x_1 & x_2 & \cdots & x_{n-1} \\ x_1^2 & x_2^2 & \cdots & x_{n-1}^2\\ \vdots & \vdots & \ddots & \vdots \\ x_1^{n-2} & x_2^{n-2} & \cdots & x_{n-1}^{n-2} \end{vmatrix}=\prod\limits_{1\leq i<j\leq {n-1}}(x_j-x_i)\)

考虑 \(n\) 阶范德蒙行列式的情形。

从第 \(n\) 行开始，自下而上地依次让每一行减去它上一行的 \(x_n\) 倍。

\(D_n=\begin{vmatrix} 1 & 1 & \cdots & 1 \\ x_1-x_n & x_2-x_n & \cdots & 0 \\ x_1(x_1-x_n) & x_2(x_2-x_n) & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots \\ x_1^{n-2}(x_1-x_n) & x_2^{n-2}(x_2-x_n) & \cdots & 0 \end{vmatrix}\)
\(=(-1)^{n+1}\begin{vmatrix}x_1-x_n & x_2-x_n & \cdots & x_{n-1}-x_n \\ x_1(x_1-x_n) & x_2(x_2-x_n) & \cdots & x_{n-1}(x_{n-1}-x_n)\\ \vdots & \vdots & \ddots & \vdots \\ x_1^{n-2}(x_1-x_n) & x_2^{n-2}(x_2-x_n) & \cdots & x_{n-1}^{n-2}(x_{n-1}-x_n) \end{vmatrix}\)

注意到行列式的一个性质：行列式 \(D\) 中某一行都乘上 \(k\) 得到 \(D_1\)，则 \(D_1=kD\)。列也是如此。那么就可以提取每一列的公因式。

\(D_n=(-1)^{n+1}(x_1-x_n)(x_2-x_n)\cdots (x_{n-1}-x_n)\begin{vmatrix}1 & 1 & \cdots & 1 \\ x_1 & x_2 & \cdots & x_{n-1}\\ \vdots & \vdots & \ddots & \vdots \\ x_1^{n-2} & x_2^{n-2} & \cdots & x_{n-1}^{n-2} \end{vmatrix}\)
\(=(-1)^{n+1}(x_1-x_n)(x_2-x_n)\cdots (x_{n-1}-x_n)D_{n-1}\)
\(=(-1)^{n+1}(x_1-x_n)(x_2-x_n)\cdots (x_{n-1}-x_n)\prod\limits_{1\leq i<j\leq {n-1}}(x_j-x_i)\)
\(=\prod\limits_{1\leq i<j\leq n}(x_j-x_i)\)

二. 矩阵

1. 定义

由 \(m\times n\) 个数排成的 \(m\) 行 \(n\) 列的矩阵数表

\(\begin{bmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{m1}&a_{m2}&\cdots&a_{mn}\end{bmatrix}\)

称为一个 \(m\times n\) 的矩阵。

矩阵的运算：加法、数乘、乘法、转置。此处略。

2. 伴随矩阵

设 \(A=(a_{ij})_n\) 为 \(n\) 阶矩阵，\(A_{ij}\) 为 \(A\) 中元素 \(a_{ij}\) 的代数余子式，\(i,j=1,2,\cdots,n\)，则称矩阵

\(\begin{bmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{bmatrix}\)

为 \(A\) 的伴随矩阵，记作 \(A^*\)。可以证明 \(AA^*=A^*A=\left|A\right|E\)。

3. 求逆矩阵

设 \(A\)&是&<span class="math inline">\(n\)&阶矩阵，若存在&<span class="math inline">\(n\)&阶矩阵&<span class="math inline">\(B\)，使得&<span class="math inline">\(AB=BA=E\)，则称矩阵&<span class="math inline">\(A\)&可逆，<span class="math inline">\(B\)&是&<span class="math inline">\(A\)&的逆矩阵，记作&<span class="math inline">\(B=A^{-1}\)。

\(A^{-1}=\frac{1}{|A|}A^*\)

构造一个 增广矩阵 \(W=[A\mid E]\)，对 \(W\) 初等变换，变成 \(W=[E\mid B]\)，可以确定 \(A^{-1}=B\)。

初等变换包括：

• 交换的矩阵的两行（列）。

• 矩阵某一行（列）的元素乘同一个不等于 \(0\) 的数。

• 矩阵某一行（列）的元素加上另一行（列）对应元素的相同倍数。

三、矩阵加速递推

HDU 4549 M斐波那契数列

题目大意：\(M\) 斐波那契数列 \(F_n\) 是一种整数数列。\(F_0=a,F_1=b,F_n=F_{n-1}\times F_{n-2}\)（\(n>1\)）。给定 \(a,b,n\)，求 \(F_n\) 的值。对 \(10^9+7\) 取模。

Solution：

\(F_0=a,F_1=b,F_n=F_{n-1}\times F_{n-2}\)，那么 \(F_2=ab,F_3=ab^2,F_4=a^2b^3,F_5=a^3b^5,F_6=a^5b^8,\cdots,F_n=a^{f_{n-1}\ \ }b^{f_n}\)。其中 \(f_i\) 为斐波那契数列的第 \(i\) 项。问题转化为如何快速求斐波那契数列。

考虑构造出转移矩阵 \(\begin{bmatrix}a_1 & a_2\\a_3 & a_4\end{bmatrix}\)，使得：

\(\begin{bmatrix}f_{n-1} & f_n\end{bmatrix}\times \begin{bmatrix}a_1 & a_2\\a_3 & a_4\end{bmatrix}=\begin{bmatrix}f_n & f_{n+1}\end{bmatrix}\)

\(f_{n-1}\times a_1+f_n\times a_3=f_n\)，\(f_{n-1}\times a_2+f_n\times a_4=f_{n+1}\)，显然 \(a_1=0,a_2=1,a_3=1,a_4=1\)。

即转移矩阵为 \(\begin{bmatrix}0 & 1\\1 & 1\end{bmatrix}\)。

\(\begin{bmatrix}f_n & f_{n+1}\end{bmatrix}=\begin{bmatrix}f_{n-1} & f_n\end{bmatrix}\times \begin{bmatrix}0 & 1\\1 & 1\end{bmatrix}=\begin{bmatrix}f_{n-2} & f_{n-1}\end{bmatrix}\times \begin{bmatrix}0 & 1\\1 & 1\end{bmatrix}^2=\begin{bmatrix}f_0 & f_1\end{bmatrix}\times \begin{bmatrix}0 & 1\\1 & 1\end{bmatrix}^n\)

\(F_n=(a^{f_{n-1}\ \ }\bmod m)\times(b^{f_n}\bmod m)\)，其中 \(m=10^9+7\)。则 \(f_n\) 和 \(f_{n-1}\) 会非常大，需要用费马小定理降幂。

\(a^p\bmod m=a^{p\bmod (m-1)}\)。所以 \(F_n=a^{(f_{n-1}\ \bmod (m-1))}b^{(f_n\bmod (m-1))}\)。

#include<bits/stdc++.h>
#define int long long 
using namespace std;
const int N=2,mod=1e9+7;
int x,y,n,f[N][N],k1,k2;
void mul(int x[N][N],int y[N][N]){    //矩阵乘法 
    int c[N][N];
    memset(c,0,sizeof(c));
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            for(int k=0;k<N;k++)
                c[i][j]=(c[i][j]+x[i][k]*y[k][j])%(mod-1);
    memcpy(x,c,sizeof(c));
}
int ksm(int x,int n,int mod){    //快速幂 
    int ans=mod!=1;
    for(x%=mod;n;n>>=1,x=x*x%mod)
        if(n&1) ans=ans*x%mod;
    return ans;
}
signed main(){
    while(~scanf("%lld%lld%lld",&x,&y,&n)){
        int a[N][N]={{0,1},{1,1}};    //转移矩阵 
        memset(f,0,sizeof(f));
        for(int i=0;i<N;i++) f[i][i]=1;    //构造单位矩阵。单位矩阵起着特殊的作用，如同数的乘法中的 1。
        for(;n;n>>=1,mul(a,a))    //矩阵快速幂 
            if(n&1) mul(f,a);
        k1=f[0][0],k2=f[0][1];    //k1=f[n],k2=f[n+1]。其中 f[i] 为斐波那契数列第 i 项。 
        printf("%lld\n",ksm(x,k1,mod)*ksm(y,k2,mod)%mod);
    }
    return 0;
}

多阶线性递推

例如 \(f_1=f_2=0\)，\(f_n=7f_{n-1}+6f_{n-2}+5n+4\times 3^n\)。这里直接给出转移：

\(\begin{bmatrix}f_n & f_{n-1} & n & 3^n &1\end{bmatrix}\times \begin{bmatrix}7 & 1 & 0 & 0 & 0\\6 & 0 & 0 & 0 & 0\\5 & 0 & 1 & 0 & 0\\12 & 0 & 0 & 3 & 0\\5 & 0 & 1 & 0 & 1\end{bmatrix}=\begin{bmatrix}f_{n+1} & f_n & n+1 & 3^{n+1} &1\end{bmatrix}\)

四、高斯消元

1. 基本思想

高斯消元是一种求解线性方程组的方法。求解这种方程组的步骤可概括成对 增广矩阵 的三类操作：

• 用一个非零的数乘某一行。

• 把其中一行的若干倍加到另一行上。

• 交换两行的位置。

我们把这三类操作成为矩阵的“初等行变换”。同理，我们也可以定义矩阵的“初等列变换”。

通过初等行变换把增广矩阵变为简化阶梯形矩阵的线性方程组求解算法就是高斯消元算法。高斯消元的算法思想就是，对于每个未知量 \(x_i\)，找到一个 \(x_i\) 的系数非零，但 \(x_1\sim x_{i-1}\) 的系数都是零的方程（当然也有可能找不到），然后用初等行变换把其他方程的 \(x_i\) 的系数全部消成零。

另外，在高斯消元的过程中，可能会遇到各种各样的特殊情形。

当消元后出现了 \(0=b\) 这样的方程，其中 \(b\) 是一个非零常数，则说明原方程组无解。

当消元后某一方程有不只一个系数非零，那么原方程组有无数解。

2. 举个栗子

\(\left\{\begin{aligned} x_1+2x_2-x_3&=-6 \notag \\ 2x_1+x_2-3x_3&=-9 \notag \\ -x_1-x_2+2x_3&=7 \notag \end{aligned}\right.\)

先把它写成增广矩阵。

\(\left\{\begin{aligned} x_1+2x_2-x_3&=-6 \notag \\ 2x_1+x_2-3x_3&=-9 \notag \\ -x_1-x_2+2x_3&=7 \notag \end{aligned}\right. \Rightarrow \begin{bmatrix}\begin{array}{ccc|c}1 & 2 & -1 & 3\\0 & -3 & -1 & 3\\-1 & -1 & 2 & 7\end{array}\end{bmatrix}\)

然后用若干次初等行变换求解上面的方程组，过程如下：

\(\begin{bmatrix}1 & 2 & -1 & -6\\2 & 1 & -3 & -9\\-1 & -1 & 2 & 7\end{bmatrix}\stackrel{r_2-2r_1}{\Longrightarrow}\begin{bmatrix}1 & 2 & -1 & -6\\0 & -3 & -1 & 3\\-1 & -1 & 2 & 7\end{bmatrix}\stackrel{r_3+r_1}{\Longrightarrow}\begin{bmatrix}1 & 2 & -1 & -6\\0 & -3 & -1 & 3\\0 & 1 & 1 & 1\end{bmatrix}\)

\(\stackrel{swap(r_2,r_3)}{\Longrightarrow}\begin{bmatrix}1 & 2 & -1 & -6\\0 & 1 & 1 & 1\\0 & -3 & -1 & 3\end{bmatrix} \stackrel{r_3+3r_2}{\Longrightarrow}\begin{bmatrix}1 & 2 & -1 & -6\\0 & 1 & 1 & 1\\0 & 0 & 2 & 6\end{bmatrix} \stackrel{r_3\times 0.5}{\Longrightarrow}\begin{bmatrix}1 & 2 & -1 & -6\\0 & 1 & 1 & 1\\0 & 0 & 1 & 3\end{bmatrix}\)

最后得到的矩阵被称为“阶梯形矩阵”，它的系数矩阵部分被称为“上三角矩阵”。这个矩阵表达的信息是：

\(\begin{bmatrix}\begin{array}{ccc|c}1 & 2 & -1 & -6\\0 & 1 & 1 & 1\\0 & 0 & 1 & 3\end{array}\end{bmatrix} \Rightarrow \left\{\begin{aligned} x_1+2x_2-x_3&=-6 \notag \\ x_2+x_3&=1 \notag \\ x_3&=3 \notag \end{aligned}\right.\)

因此，我们已经知道了最后一个未知量的值，从下往上依次代回方程组，即可得到每个未知量的解。事实上，该矩阵也可以进一步化简：

\(\begin{bmatrix}1 & 2 & -1 & -6\\0 & 1 & 1 & 1\\0 & 0 & 1 & 3\end{bmatrix}\stackrel{r_1+r_3,r_2-r_3}{\Longrightarrow} \begin{bmatrix}1 & 2 & 0 & -3\\0 & 1 & 0 & -2\\0 & 0 & 1 & 3\end{bmatrix}\stackrel{r_1-2r_2}{\Longrightarrow} \begin{bmatrix}\begin{array}{ccc|c}1 & 0 & 0 & 1\\0 & 1 & 0 & -2\\0 & 0 & 1 & 3\end{array}\end{bmatrix}\)

最后得到的矩阵被称为“简化阶梯形矩阵”，它的系数矩阵部分是一个“对角矩阵”。该矩阵已经直接给出了方程的解。

//Luogu P3389
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=110;
int n;
double a[N][N];
signed main(){
    scanf("%lld",&n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            scanf("%lf",&a[i][j]);    //增广矩阵 
    for(int i=1;i<=n;i++){    //消第 i 个元
        int p=i;
        for(int j=i;j<=n;j++)    //找到 x[i] 的系数不为 0 的一个方程
            if(fabs(a[j][i])>1e-8){p=j;break;}
        for(int j=1;j<=n+1;j++)
            swap(a[i][j],a[p][j]);
        if(fabs(a[i][i])<1e-8) puts("No Solution"),exit(0);    //不存在唯一解 
        for(int j=1;j<=n;j++){    //消去其他方程的 x[i] 的系数
            if(i==j) continue;
            double x=a[j][i]/a[i][i];
            for(int k=i;k<=n+1;k++)
                a[j][k]-=a[i][k]*x;
        }
    }
    for(int i=1;i<=n;i++)
        printf("%.2lf\n",a[i][n+1]/a[i][i]);
    return 0;
} 

五、生成树计数

1. 前置概念

\(G\) 的度数矩阵 \(D[G]\) 是一个 \(n\times n\) 的矩阵，并且满足：当 \(i\neq j\) 时，\(d_{ij}=0\)；当 \(i=j\) 时，\(d_{ij}\) 等于 \(v_i\) 的度数。

\(G\) 的邻接矩阵 \(A[G]\) 也是一个 \(n\times n\) 的矩阵，并且满足：如果 \(v_i\)、\(v_j\) 之间有边直接相连，则 \(a_{ij}=1\)，否则为 \(0\)。

定义 \(G\) 的 Kirchhoff 矩阵（也称为拉普拉斯算子）\(C[G]\) 为 \(C[G]=D[G]-A[G]\)。

2. Matrix-Tree 定理

\(G\) 的所有不同的生成树的个数等于其 Kirchhoff 矩阵 \(C[G]\) 任何一个 \(n-1\) 阶主子式的行列式的绝对值。

证明省略。但在证明中我们用到了一个矩阵，称为关联矩阵。

关联矩阵是一个 \(n\) 行 \(m\) 列的矩阵，行对应点而列对应边。

关联矩阵满足，如果存在一条边 \(e=\{v_i,v_j\}\)，那在 \(e\) 所对应的列中，\(v_i\) 和 \(v_j\) 所对应的那两行，一个为 \(1\)、另一个为 \(-1\)，其他的均为 \(0\)。至于哪个是 \(1\) 哪个是 \(-1\) 并不重要。

我们令关联矩阵为 \(B\)，考虑 \(BB^T\)，可以发现 \(C=BB^T\)。

//SP104 HIGH - Highways
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=20;
int t,n,m,x,y;
double a[N][N],ans; 
void Gauss(int n){    //高斯消元 
    for(int i=1;i<=n;i++){
        int p=i;
        for(int j=i;j<=n;j++)
            if(fabs(a[j][i])>1e-8){p=j;break;}
        for(int j=1;j<=n+1;j++)
            swap(a[i][j],a[p][j]);
        if(fabs(a[i][i])<1e-8) return ;
        for(int j=1;j<=n;j++){
            if(i==j) continue;
            double x=a[j][i]/a[i][i];
            for(int k=i;k<=n+1;k++)
                a[j][k]-=a[i][k]*x;
        }
    }
}
signed main(){
    scanf("%lld",&t);
    while(t--){
        memset(a,0,sizeof(a)),ans=1;
        scanf("%lld%lld",&n,&m);
        for(int i=1;i<=m;i++){
            scanf("%lld%lld",&x,&y);
            a[x][x]++,a[y][y]++;
            a[x][y]--,a[y][x]--;    //Kirchhoff 矩阵 
        }
        n--,Gauss(n);    //利用高斯消元将 Kirchhoff 矩阵消为对角矩阵
        for(int i=1;i<=n;i++)
            ans=ans*a[i][i]; 
        printf("%.0lf\n",fabs(ans));    //高斯消元后，由于已经消成了对角矩阵，所以对角线乘积的绝对值就是答案
    } 
    return 0;  
}  

六、习题

矩阵加速递推：

• BZOJ 4887「TJOI 2007」可乐
• LOJ 6208 树上询问（这个代码 过了是什么情况 ¿）

高斯消元：

• HDU 2408 String Equations
• BZOJ 3143「HNOI 2013」游走

生成树计数：

• HDU 4408 Minimum Spanning Tree
• HDU 6836 Expectation

一、高斯消元

for(int i=1;i<=n;i++){
	int p=i;
	for(int j=i;j<=n;j++) if(fabs(a[j][i])>fabs(a[p][i])) p=j;	//避免精度误差
	if(p!=i) swap(a[i],a[p]),swap(b[i],b[p]);
	for(int j=i+1;j<=n;j++){
		double x=-a[j][i]/a[i][i];
		for(int k=i;k<=n;k++) a[j][k]+=a[i][k]*x; b[j]+=b[i]*x;	//注意 b 也要消
	}
}
for(int i=n;i>=1;i--){
	b[i]/=a[i][i];
	for(int j=i-1;j>=1;j--) b[j]-=a[j][i]*b[i];
}

辗转相除：

避免 a[i][i] 没有逆元的情况，也不会被 double 卡没。对于两行，不断将主元较大的那行减去主元较小的，最终一定有一行主元为 \(0\)。

int Gauss(int n){
	int ans=1;
	for(int i=1;i<=n;i++){
		for(int j=i+1;j<=n;j++)
			while(a[j][i]){
				int x=a[i][i]/a[j][i];
				for(int k=1;k<=n;k++) a[i][k]=(a[i][k]-a[j][k]*x%mod+mod)%mod;
				swap(a[i],a[j]),ans=mod-ans;
			}
		ans=ans*a[i][i]%mod;
	}
	return (ans+mod)%mod;
}

行列式：消成上三角矩阵，对角线乘起来。交换两行，行列式取反。

解决有后效性 DP：

• 设计完状态、写出转移方程后，发现状态间不满足 DAG 的性质，无法找到一个顺序来依次计算 DP 值。
• 转移只涉及到状态的一次项。

1. P3232 [HNOI2013]游走

2022.7.21

给出一张 \(n\) 个点 \(m\) 条边的无向连通图，小 Z 从 \(1\) 出发，每次随机选择某条边走到下一顶点，并将 \(ans\) 加上这条边的编号，走到 \(n\) 时结束。

要求对 \(n\) 进行编号，使得总分的期望值最小。求最小值。

\(2\leq n\leq 500\)，\(1\leq m\leq 125000\)，无重边无自环。

肯定是给期望经过次数最多的边编号为 \(1\)，第二多编号为 \(2\)，以此类推。

边的期望可以由点的期望转化而来。设 \(f_x\) 表示点 \(x\) 的期望经过次数，\(f_x=[x=1]+\sum_{(x,y)\in E,y\neq n}\frac{f_y}{deg_y}\)。可以将所有 \(f\) 移到左边，常数项移到右边，高斯消元。

怎么将点的期望转化为边的期望？设 \(g_i\) 表示第 \(i\) 条边的期望经过次数，\(g_i=\large\frac{f_{x_i}}{deg_{x_i}}\normalsize+\large\frac{f_{y_i}}{deg_{y_i}}\normalsize\)。

#include<bits/stdc++.h>
using namespace std;
const int N=510,M=2e5+5;
int n,m,x,y,u[M],v[M],d[N];	//某 sb 开了 u[N],v[N] wa 了调不出来
double a[N][N],b[N],f[M],ans;
vector<int>g[N];
void Gauss(int n){
	for(int i=1;i<=n;i++){
		int p=i;
		for(int j=i;j<=n;j++) if(fabs(a[j][i])>fabs(a[p][i])) p=j;
		if(p!=i) swap(a[i],a[p]),swap(b[i],b[p]);
		for(int j=i+1;j<=n;j++){
			double x=-a[j][i]/a[i][i];
			for(int k=i;k<=n;k++) a[j][k]+=a[i][k]*x; b[j]+=b[i]*x;
		}
	}
	for(int i=n;i>=1;i--){
		b[i]/=a[i][i];
		for(int j=i-1;j>=1;j--) b[j]-=a[j][i]*b[i];
	}
}
signed main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=m;i++){
		scanf("%d%d",&x,&y),d[x]++,d[y]++,u[i]=x,v[i]=y;
		g[x].push_back(y),g[y].push_back(x);
	}
	for(int x=1;x<n;x++){
		a[x][x]=1;
		for(int y:g[x]) if(y!=n) a[x][y]-=1.0/d[y];
	}
	b[1]=1,Gauss(n-1);
	for(int i=1;i<=m;i++) f[i]=b[u[i]]/d[u[i]]+b[v[i]]/d[v[i]];
	sort(f+1,f+1+m,greater<double>());
	for(int i=1;i<=m;i++) ans+=f[i]*i;
	printf("%.3lf\n",ans); 
	return 0;
}

2. P4035 [JSOI2008]球形空间产生器

2022.7.21

给出一个 \(n\) 维球体球面上的 \(n+1\) 个点，求这个球体球心的位置。

\(1\leq n\leq 10\)，\(|a_{i,j}|\leq 20000\)。

设球心为 \((x_1,x_2,\cdots,x_n)\)，到球面上点的距离为 \(d\)。

\[\begin{cases} (x_1-a_{1,1})^2+(x_2-a_{1,2})^2+\cdots+(x_n-a_{1,n}^2)=d^2\\ (x_1-a_{2,1})^2+(x_2-a_{2,2})^2+\cdots+(x_n-a_{2,n}^2)=d^2\\ \cdots\\ (x_1-a_{n+1,1})^2+(x_2-a_{n+1,2})^2+\cdots+(x_n-a_{n+1,n}^2)=d^2 \end{cases} \]

二次方程不好做。拆开：

\[\begin{cases} x_1-2x_1a_{1,1}+a_{1,1}^2+x_2^2-2x_2a_{1,2}+a_{1,2}^2+\cdots+x_n^2-2x_na_{1,n}+a_{1,n}^2=d^2\\ x_1-2x_1a_{2,1}+a_{2,1}^2+x_2^2-2x_2a_{2,2}+a_{2,2}^2+\cdots+x_n^2-2x_na_{2,n}+a_{2,n}^2=d^2\\ \cdots\\ x_1-2x_1a_{n+1,1}+a_{n+1,1}^2+x_2^2-2x_2a_{n+1,2}+a_{n+1,2}^2+\cdots+x_n^2-2x_na_{n+1,n}+a_{n+1,n}^2=d^2 \end{cases} \]

显然每个方程左边都有 \(x_1^2+\cdots+x_n^2\)，右边都有 \(d^2\)。考虑将第 \(1\sim n\) 个方程分别减去第 \(n+1\) 个方程，就得到了一次的方程组。

\[\begin{cases} 2(a_{n+1,1}-a_{1,1})x_1+\cdots+2(a_{n+1,n}-a_{1,n})x_n=a_{n+1,1}^2-a_{1,1}^2+\cdots+a_{n+1,n}^2-a_{1,n}^2\\ \cdots\\ 2(a_{n+1,1}-a_{n,1})x_1+\cdots+2(a_{n+1,n}-a_{n,n})x_n=a_{n+1,1}^2-a_{n,1}^2+\cdots+a_{n+1,n}^2-a_{n,n}^2 \end{cases} \]

由于 \(a\) 是给出的，就可以直接高斯消元了。

#include<bits/stdc++.h>
using namespace std;
const int N=15;
int n;
double a[N][N],b[N];
signed main(){
	scanf("%d",&n);
	for(int i=1;i<=n+1;i++)
		for(int j=1;j<=n;j++) scanf("%lf",&a[i][j]);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			b[i]+=a[n+1][j]*a[n+1][j]-a[i][j]*a[i][j],a[i][j]=2*(a[n+1][j]-a[i][j]);
	for(int i=1;i<=n;i++){
		int p=i;
		for(int j=i;j<=n;j++) if(fabs(a[j][i])>fabs(a[p][i])) p=j;
		if(p!=i) swap(a[i],a[p]),swap(b[i],b[p]);
		for(int j=i+1;j<=n;j++){
			double x=-a[j][i]/a[i][i];
			for(int k=i;k<=n;k++) a[j][k]+=a[i][k]*x; b[j]+=b[i]*x;	//注意 b 也要消
		}
	}
	for(int i=n;i>=1;i--){
		b[i]/=a[i][i];
		for(int j=i-1;j>=1;j--) b[j]-=a[j][i]*b[i];
	}
	for(int i=1;i<=n;i++) printf("%.3lf ",b[i]);
	return 0;
}

3. P6125 [JSOI2009] 有趣的游戏（结合 AC 自动机）

2022.7.21 概率转化为期望

给出 \(n\) 个长度为 \(l\) 的字符串 \(a_{1\sim n}\)，由前 \(m\) 个大写字母组成。现在随机生成一个字母序列，每次第 \(i\) 个字母有 \(\frac{p_i}{q_i}\) 的概率生成。求每个 \(a_i\) 求它作为首个出现在字母序列中的字符串的概率。

\(1\leq n,l,m\leq 10\)，\(0\leq p_i\leq q_i\leq 10\) 且 \(\gcd(p,q)=1\)。

对 \(a_{1\sim n}\) 建立 AC 自动机，考虑在 AC 自动机上 DP，变为求 \(a_i\) 对应节点作为它们中首个出现的节点的概率。

由于一个点可能被经过多次，每个节点被经过的概率不好做。所以考虑转化为求每个点经过次数的期望，且到达任意一个 \(a_i\) 对应节点就不能再走。由于到达 \(a_i\) 对应节点就不能走了，即这个节点只会经过一次，因此在 \(a_i\) 对应节点停下的期望就是在这里停下的概率。

类似 P3232 [HNOI2013]游走，列出方程然后高斯消元。若 \(x\) 不是 \(a_i\) 对应节点，从节点 \(x\) 走到它第 \(i\) 个儿子的概率是 \(\frac{p_i}{q_i}\)。且只有根节点等式右边为 \(1\)，其余都为 \(0\)。

#include<bits/stdc++.h>
using namespace std;
const int N=110;
int n,len,m,x,y,tot=1,tg[N],id[N],ch[N][11],fail[N];
double p[N],a[N][N],b[N];
char s[N];
queue<int>q; 
void getfail(){
	for(int i=0;i<m;i++) ch[0][i]=1;
	q.push(1),fail[1]=0;
	while(q.size()){
		int x=q.front(),y;q.pop();
		for(int i=0;i<m;i++){
			if((y=ch[x][i])) fail[y]=ch[fail[x]][i],q.push(y);
			else ch[x][i]=ch[fail[x]][i];
		}
	}
}
void Gauss(int n){
	for(int i=1;i<=n;i++){
		int p=i;
		for(int j=i;j<=n;j++) if(fabs(a[j][i])>fabs(a[p][i])) p=j;
		if(p!=i) swap(a[i],a[p]),swap(b[i],b[p]);
		for(int j=i+1;j<=n;j++){
			double x=-a[j][i]/a[i][i];
			for(int k=i;k<=n;k++) a[j][k]+=a[i][k]*x; b[j]+=b[i]*x;
		}
	}
	for(int i=n;i>=1;i--){
		b[i]/=a[i][i];
		for(int j=i-1;j>=1;j--) b[j]-=a[j][i]*b[i];
	}
}
signed main(){
	scanf("%d%d%d",&n,&len,&m);
	for(int i=0;i<m;i++)
		scanf("%d%d",&x,&y),p[i]=1.0*x/y;
	for(int i=1;i<=n;i++){
		scanf("%s",s+1);
		int p=1;
		for(int j=1,k;j<=len;j++){
			if(!ch[p][k=s[j]-'A']) ch[p][k]=++tot;
			p=ch[p][k];
		}
		tg[id[i]=p]=1;
	}
	getfail();
	for(int i=1;i<=tot;i++){
		a[i][i]+=1;
		if(!tg[i]) for(int j=0;j<m;j++) a[ch[i][j]][i]-=p[j];
	}
	b[1]=1,Gauss(tot);
	for(int i=1;i<=n;i++) printf("%.2lf\n",b[id[i]]);
	return 0;
}

二、行列式

• 定义：

  对于一个 \(n\times n\) 的矩阵 \(A\)，定义其行列式为：

  \[\det A= \begin{vmatrix} A_{1,1} & A_{1,2} & \cdots & A_{1,n} \\ A_{2,1} & A_{2,2} & \cdots & A_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ A_{n,1} & A_{n,2} & \cdots & A_{n,n} \end{vmatrix} =\sum_{p_1,p_2,\cdots,p_n}(-1)^{r(p)}\prod_{i=1}^n A_{i,p_i} \]

  其中 \(p\) 为一个 \(1\sim n\) 的排列，\(r(p)\) 表示 \(p\) 的逆序对数。

• 性质：略。

• 范德蒙德行列式

  \[D_n= \begin{vmatrix} 1&1&\cdots&1\\ x_1&x_2&\cdots&x_n\\ x_1^2&x_2^2&\cdots&x_n^2\\ \vdots&\vdots&\ddots&\vdots\\ x_1^{n-1}&x_2^{n-1}&\cdots&x_n^{n-1}\\ \end{vmatrix} =\prod_{1\leq i<j\leq n}(x_j-x_i) \]

  证明：数学归纳法。

• 求法：

  用高斯消元将矩阵消成上三角。由于交换两行，行列式只有符号会变化，记录一下交换了几次即可。

  消元时如果开 double 可能会爆精度，如果模数不是质数，有一个辗转消元法：假如要消掉 \(A_{j,i}\)，令 \(d=\lfloor\large\frac{A_{i,i}}{A_{j,i}}\normalsize\rfloor\)，不断将第 \(j\) 行乘以 \(-d\) 倍后加到第 \(i\) 行，并交换 \(i,j\) 两行，直到 \(A_{j,i}=0\) 为止。每个元素最多被操作 \(\log p\) 次，复杂度 \(n^2\log p+n^3\)。

  //Luogu P7112
  #include<bits/stdc++.h>
  using namespace std;
  const int N=610;
  int n,mod,a[N][N],ans=1;
  signed main(){
  	scanf("%d%d",&n,&mod);
  	for(int i=1;i<=n;i++)
  		for(int j=1;j<=n;j++) scanf("%d",&a[i][j]); 
  	for(int i=1;i<=n;i++)
  		for(int j=i+1;j<=n;j++)
  			while(a[j][i]){
  				int x=a[i][i]/a[j][i];
  				for(int k=i;k<=n;k++) a[i][k]=(a[i][k]+1ll*(mod-a[j][k])*x)%mod;
  				//注意减少取模避免 TLE 
  				swap(a[i],a[j]),ans*=-1;
  			}
  	for(int i=1;i<=n;i++) ans=1ll*ans*a[i][i]%mod;
  	printf("%d\n",(ans+mod)%mod);
  	return 0;
  } 
  

  模数是质数：

  int ans=1;
  for(int i=1;i<=n;i++){
  	int p=i;
  	for(int j=i;j<=n;j++) if(a[j][i]){p=j;break;}
  	if(p!=i) swap(a[i],a[p]),ans=mod-ans;
  	for(int j=i+1,iv=qpow(a[i][i],mod-2);j<=n;j++){
  		int x=1ll*a[j][i]*iv%mod;
  		for(int k=i;k<=n;k++) a[j][k]=(a[j][k]+1ll*(mod-a[i][k])*x%mod)%mod;
  	}
  	ans=1ll*ans*a[i][i]%mod;
  }
  

三、矩阵树定理

• 给出一个无权无向图 \(G\)，设 \(D\) 为度数矩阵（\(d_{i,i}\) 等于节点 \(i\) 的度数），\(A\) 为邻接矩阵。

  定义 \(G\) 的基尔霍夫矩阵 \(K=D-A\)，然后令 \(K'\) 表示 \(K\) 去掉第 \(k\) 行与第 \(k\) 列（\(k\) 任意）的结果（\(n-1\) 阶主子式），则 \(\text{det}(K')\) 记为 \(G\) 的生成树个数。

• 边带权：

  度数矩阵变成了相邻边的权值和。加入一条无向边 \((x,y,z)\) 时，a[x][x]+=z,a[y][y]+=z,a[x][y]-=z,a[y][x]-=z。

• 边有向：

  构造基尔霍夫矩阵时，如果外向树（从根向外），那么 \(D\) 就改成度入矩阵，如果是内向树（从外向根），则 \(D\) 改成出度矩阵。 （同样可以加权！）

  如果要以 \(k\) 为根，那么要把第 \(k\) 行和第 \(k\) 列删掉。

  加入一条有向边 \((x,y,z)\) 时（以外向树为例），a[y][y]+=z,a[x][y]-=z。

• 强制某条边不能出现在树形图中，求方案数：算树形图总数 - 去掉这条边后的树形图个数。

• 算一条边出现在多少树形图中：将这条边缩起来，或容斥转成不在。

1. LOJ#6271. 「长乐集训 2017 Day10」生成树求和 加强版

2024.4.23 矩阵树定理

给出一张 \(n\) 个点 \(m\) 条边的无向图 \(G\)，边有边权。

定义一棵生成树的权值为：它所包含的所有边的边权按三进制不进位加法相加得到的数。

求 \(G\) 的所有生成树的权值在十进制下的和 \(\bmod 10^9+7\)。

\(n\leq 100\)，\(m\leq \frac{n(n-1)}{2}\)，\(0\leq w_i\leq 10^4\)，保证无重边无自环。

令 \(1\) 的边权值为 \(1\)，\(2\) 的边权值为 \(x\)，\(3\) 的边权值为 \(x^2\)，把多项式运算放在 \(\bmod x^3-1\) 的意义下进行。

逆元：手搓。三个未知数（\(b_0,b_1,b_2\)）三个方程。

\[(a_0x^0+a_1x^1+a_2x^2)(b_0x^0+b_1x^1+b_2x^2)\equiv 1\pmod{x^3-1}\\ \Rightarrow \begin{cases} a_0b_0+a_1b_2+a_2b_1=1\\ a_0b_1+a_1b_0+a_2b_2=0\\ a_0b_2+a_1b_1+a_2b_0=0 \end{cases} \\ \Rightarrow \begin{cases} b_0=\dfrac{a_0^2-a_1a_2}{a_0^3+a_1^3+a_2^3-3a_0a_1a_2}\\ b_1=\dfrac{a_2^2-a_0a_1}{a_0^3+a_1^3+a_2^3-3a_0a_1a_2}\\ b_2=\dfrac{a_1^2-a_0a_2}{a_0^3+a_1^3+a_2^3-3a_0a_1a_2} \end{cases} \]

特判：\(a_0=a_1=a_2\) 时无逆元。

2. 求和

2024.4.23

给出一张 \(n\) 个点 \(m\) 条边的无向图 \(G\)，边有边权。

定义一棵生成树的权值为：所有边权在三进制下，每一位做不进位相乘后的三进制数的值。

求 \(G\) 的所有生成树的权值在十进制下的和 \(\bmod 10^9+7\)。

\(n\leq 100\)，\(1\leq w_i\leq 10^4\)，保证无重边无自环。

拆位，考虑三进制下每一位的贡献。

• 权值非 \(0\) 的生成树一定不包含边权为 \(0\) 的边。把边权为 \(0\) 的边扔掉，对剩下的图求生成树权值和。

• 边权为 \(1\)：等于没有。

  边权为 \(2\)：\(2^1\equiv 2\pmod 3\)，\(2^2\equiv 1\pmod 3\)。

  权值只取决于边权为 \(2\) 的边数量的奇偶性（若为奇，则权值为 \(2\)；否则权值为 \(1\)）。

做法 1：

• 令 \(1\) 的边权值为 \(1\)，\(2\) 的边权值为 \(x\)，跑矩阵树定理，把得到的多项式按奇偶次项统计一下。时间复杂度 \(\mathcal O(n^5\log_3 V)\) 或 \(\mathcal O(n^4\log n\log_3 V)\)，可以获得 \(50\) 分。

  取 \(n\) 个点值后拉格朗日插值。时间复杂度 \(\mathcal O(n^4\log V)\)，可以获得 \(70\) 分。

  不插值，把多项式运算放在 \(\bmod x^2-1\) 的意义下进行，时间复杂度 \(\mathcal O(n^3\log_3 V)\)。可以获得 \(100\) 分。

  求逆元：\((a+bx)(c+dx)\equiv 1\pmod{x^2-1}\)，即 \((ac+bd)+(ad+bc)x=1\)，解 \(\begin{cases}ac+bd=1\\ad+bc=0\end{cases}\) 得 \(d=\frac{b}{b^2-a^2},c=-\frac{ad}{b}\)。特判 \(b=0\)（\(c=\frac 1 a,d=0\)）和 \(a=b\)（\(a=b\) 时方程组无解，故无逆元）即可。

  如果高消时一整列都没有逆元怎么办？说明一整列都是 \(a=b\) 的形式，直接当成一个数消就好了。特判掉。

做法 2：

• 设偶数条边权为 \(2\) 的边的方案数为 \(a\)，奇数条边权为 \(2\) 的边的方案数为 \(b\)。

  令 \(1\) 的边权值为 \(1\)，\(2\) 的边权值为 \(1\)，跑矩阵树定理，可以得到 \(a+b\)。

  令 \(1\) 的边权值为 \(1\)，\(2\) 的边权值为 \(-1\)，跑矩阵树定理，可以得到 \(a-b\)。

  解方程就能求出 \(a,b\) 了，答案为 \(a+2b\)。

  时间复杂度 \(\mathcal O(n^3\log_3 V)\)。可以获得 \(100\) 分。

做法 1：

#include<bits/stdc++.h>
using namespace std;
const int N=110,mod=1e9+7;
int n,m,x,y,z,w[N][N],ans;
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
struct num{
	int x,y;
	num operator+(num a){return {(x+a.x)%mod,(y+a.y)%mod};}
	num operator-(num a){return {(x-a.x)%mod,(y-a.y)%mod};}
	num operator*(num a){return {(1ll*x*a.x+1ll*y*a.y)%mod,(1ll*x*a.y+1ll*y*a.x)%mod};}
	bool ok(){return (x+mod)%mod!=(y+mod)%mod;}
	num inv(){
		if(!y) return {qpow(x,mod-2),0};
		int d=1ll*y*qpow((1ll*y*y-1ll*x*x)%mod,mod-2)%mod;
		return {-1ll*x*d%mod*qpow(y,mod-2)%mod,d};
	}
}a[N][N];
void add(int x,int y,num z){
	a[x][x]=a[x][x]+z,a[x][y]=a[x][y]-z;
}
num det(int n){
	num ans={1,0};
	for(int i=1;i<=n;i++){
		int p=0;
		for(int j=i;j<=n;j++)
			if((a[j][i].x+mod)%mod!=(a[j][i].y+mod)%mod){p=j;break;}
		if(p){
			if(p!=i) swap(a[i],a[p]),ans={-ans.x,-ans.y};
			num iv=a[i][i].inv();
			for(int j=i+1;j<=n;j++){ 
				num x=a[j][i]*iv;
				for(int k=i;k<=n;k++) a[j][k]=a[j][k]-a[i][k]*x;
			} 
		}
		else{	//一整列都没有逆元
			for(int j=i;j<=n;j++) if(a[j][i].x){p=j;break;}
			if(!p) return {0,0};
			if(p!=i) swap(a[i],a[p]),ans={-ans.x,-ans.y};
			for(int j=i+1,iv=qpow(a[i][i].x,mod-2);j<=n;j++){
				num x={1ll*a[j][i].x*iv%mod,0};
				for(int k=i;k<=n;k++) a[j][k]=a[j][k]-a[i][k]*x;
			}
		}
		ans=ans*a[i][i];
	}
	return ans;
}
signed main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=m;i++)
		scanf("%d%d%d",&x,&y,&z),
		w[x][y]=w[y][x]=z;
	for(int k=0,pw=1;k<9;k++,pw=3ll*pw%mod){
		for(int i=1;i<=n;i++)
			for(int j=1;j<=n;j++) a[i][j]={0,0};
		for(int i=1;i<=n;i++)
			for(int j=1;j<=n;j++){
				int x=w[i][j]%3; w[i][j]/=3;
				if(x==1) add(i,j,{1,0});
				if(x==2) add(i,j,{0,1});
			}
		num x=det(n-1);
		ans=(ans+(x.x+2ll*x.y)*pw)%mod;
	}
	printf("%d\n",(ans+mod)%mod);
	return 0;
}

做法 2：

#include<bits/stdc++.h>
using namespace std;
const int N=110,mod=1e9+7;
int n,m,x,y,z,w[N][N],iv2=(mod+1)/2,ans;
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
struct MT{
	int a[N][N];
	void clear(){memset(a,0,sizeof(a));}
	void add(int x,int y,int z){a[x][x]+=z,a[x][y]-=z;}
	int det(int n){
		int ans=1;
		for(int i=1;i<=n;i++){
			int p=i;
			for(int j=i;j<=n;j++) if(a[j][i]){p=j;break;}
			if(i!=p) swap(a[i],a[p]),ans=mod-ans;
			int iv=qpow(a[i][i],mod-2);
			for(int j=i+1;j<=n;j++)
				for(int k=i,x=1ll*a[j][i]*iv%mod;k<=n;k++)
					a[j][k]=(a[j][k]-1ll*x*a[i][k])%mod;
			ans=1ll*ans*a[i][i]%mod;
		}
		return ans;
	}
}A,B;
signed main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=m;i++)
		scanf("%d%d%d",&x,&y,&z),
		w[x][y]=w[y][x]=z;
	for(int k=0,pw=1;k<9;k++,pw=3ll*pw%mod){
		A.clear(),B.clear();
		for(int i=1;i<=n;i++)
			for(int j=1;j<=n;j++){
				int x=w[i][j]%3; w[i][j]/=3;
				if(x==1) A.add(i,j,1),B.add(i,j,1);
				if(x==2) A.add(i,j,1),B.add(i,j,-1);
			}
		int x=A.det(n-1),y=B.det(n-1);
		tie(x,y)=make_pair(1ll*(x+y)*iv2%mod,1ll*(x-y)*iv2%mod);
		ans=(ans+(x+2ll*y)*pw)%mod;	//!!! 2ll*
	}
	printf("%d\n",(ans+mod)%mod);
	return 0;
}

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

给出一张无向图，求它的所有生成树 \(T\) 的价值之和对 \(998244353\) 取模后的结果。

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

其中 \(e_1,e_2,\cdots,e_{n-1}\) 为 \(T\) 包含的边的编号。

\(1\leq n\leq 30\)，\(1\leq m\leq \frac{n(n-1)}{2}\)，\(1\leq w_i\leq 152501\)。

根据 \(\sum_{d\mid n}\varphi(d)=n\)，改写为

\[\sum_{T}\left(\sum_{i=1}^{n-1}w_{e_i}\right)\times\left(\sum_{d|\gcd(w_{e_1},w_{e_2},\dots,w_{e_{n-1}})}\varphi(d)\right)=\sum_{d=1}^{\max w}\varphi(d)\times\left(\sum_{T,d\mid\gcd(w_{e_1},w_{e_2},\dots,w_{e_{n-1}})}\sum_{i=1}^{n-1}w_{e_i}\right) \]

后面括号里的部分，相当于原图边权为 \(d\) 的倍数的边组成的子图中，所有生成树的边权之和。考虑 Matrix-Tree 定理。一般来说它解决的是求所有生成树边权 乘积 之和，而这里要求生成树边权 和 之和。

一个朴素做法是考虑每条边的贡献，也就是包含这条边的生成树数量，可以用原图生成树数量减去去掉这条边后数量求得，复杂度 \(\mathcal O(w\cdot mn^3)\)，无法通过本题。

有个挺妙的套路，我们将所有边的边权看作一个一次函数 \(y=1+wx\)，那么边权和就是生成树上所有一次函数乘积的一次项。考虑怎么定义矩阵行列式中涉及加减乘除运算。首先由于我们最终只关系函数的一次项，可以把所有运算放在 \(\bmod x^2\) 的意义下进行。加减直接加，乘法 \((a+bx)(c+dx)=ac+(ad+bc)x\)，除法稍微有点麻烦，首先 \(\large\frac{a+bx}{c+dx}\) 的常数项必须是 \(\large\frac a c\)，因为只有 \(\large\frac a c\) 乘上 \(c+dx\) 后常数项才能得到 \(a\)，待定系数法算一下也即可得到一次项是 \(\large\frac{bc-ad}{c^2}\)，即 \(\large\frac{a+bx}{c+dx}=\frac{a}{c}+\frac{bc-ad}{c^2}\normalsize x\)，写个结构体维护一下即可。

时间复杂度 \(\mathcal O(wn^3)\)。为避免卡常，可以加一个小优化：如果边权为 \(i\) 的倍数的边数 \(<n-1\) 就直接令边权和为 \(0\)，优化后复杂度变为 \(\mathcal O(\large\frac{md(w_i)}{n-1}\normalsize n^3)\)。

#include<bits/stdc++.h>
using namespace std;
const int N=460,W=152505,mod=998244353;
int n,m,w,x[N],y[N],z[N],cnt,p[W],phi[W],ans;
bool vis[W];
vector<int>e[W];
int mul(int x,int n,int mod){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod)
		if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
struct node{
	int x,y;
	friend node operator+(node a,node b){return {(a.x+b.x)%mod,(a.y+b.y)%mod};}
	friend node operator-(node a,node b){return {(a.x-b.x+mod)%mod,(a.y-b.y+mod)%mod};}
	friend node operator*(node a,node b){return {1ll*a.x*b.x%mod,(1ll*a.x*b.y%mod+1ll*a.y*b.x%mod)%mod};}
	friend node operator/(node a,node b){int inv=mul(b.x,mod-2,mod); return {1ll*a.x*inv%mod,1ll*(1ll*a.y*b.x%mod-1ll*a.x*b.y%mod+mod)%mod*inv%mod*inv%mod};} 
}a[N][N];
node det(){
	node ans={1,0};
	for(int i=2;i<=n;i++){
		int p=i;
		for(int j=i;j<=n;j++) if(a[j][i].x){p=j;break;}
		if(p!=i) swap(a[i],a[p]),ans={mod-ans.x,mod-ans.y};
		node inv=(node){1,0}/a[i][i];
		for(int j=i+1;j<=n;j++){ 
			node x=a[j][i]*inv;
			for(int k=i;k<=n;k++) a[j][k]=a[j][k]-a[i][k]*x;
		} 
		ans=ans*a[i][i];
	}
	return ans;
}
signed main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=m;i++){
		scanf("%d%d%d",&x[i],&y[i],&z[i]),w=max(w,z[i]);
		for(int j=1;j*j<=z[i];j++) if(z[i]%j==0){
			e[j].push_back(i);
			if(j*j!=z[i]) e[z[i]/j].push_back(i);
		}
	}
	vis[0]=vis[1]=1,phi[1]=1;
	for(int i=2;i<=w;i++){
		if(!vis[i]) p[++cnt]=i,phi[i]=i-1;
		for(int j=1;j<=cnt&&i*p[j]<=w;j++){
			vis[i*p[j]]=1;
			if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;} 
			phi[i*p[j]]=phi[i]*phi[p[j]];
		}
	}
	for(int i=1;i<=w;i++) if((int)e[i].size()>=n-1){
		for(int j=1;j<=n;j++)
			for(int k=1;k<=n;k++) a[j][k]={0,0};
		for(int j:e[i]){ 
			a[x[j]][x[j]]=a[x[j]][x[j]]+(node){1,z[j]},a[y[j]][y[j]]=a[y[j]][y[j]]+(node){1,z[j]};
			a[x[j]][y[j]]=a[x[j]][y[j]]-(node){1,z[j]},a[y[j]][x[j]]=a[y[j]][x[j]]-(node){1,z[j]};
		} 
		ans=(ans+1ll*phi[i]*det().y%mod)%mod;
	}
	printf("%d\n",ans);
	return 0;
} 

4. P5406 [THUPC2019]找树

FWT + Matrix-Tree 定理

给出一张 \(n\) 个点 \(m\) 条边的无向图。定义一棵生成树的权值为，边权每一位 and/or/xor 的结果，即对于每个二进制位，给出这一位进行哪种运算。

求所有生成树的最大权值。如果图不连通输出 -1。

\(1\leq n\leq 70\)，\(1\leq m\leq 5000\)，二进制位的个数 \(1\leq w\leq 12\)。

首先，这不太好通过最优化问题的思路解决，只好按照计数题的思路，对于每个 \(i\)，求出权值为 \(i\) 的生成树的数量 \(cnt_i\)，答案就是最大的 \(cnt_i\neq 0\) 的 \(i\)。现在问题就转化为求 \(cnt_i\)。

我们知道，一般来说矩阵树定理解决的是求所有生成树边权乘积之和。假如是权值是边权和：把一条边权为 \(w\) 的边看作一个幂级数 \(x^w\)，矩阵树定理求出一个幂级数 \(F(x)\)，则 \(cnt_i=[x^i]F(x)\)。

将加法换成二进制运算后，原来加法卷积变成了位运算卷积。如果运算全部都是 and/or/xor，就可以 FWT 了。不过虽然每一位的运算规律都不同，但由于它们彼此之间互相独立且都是二进制运算，FWT 对该运算依然使用。

时间复杂度 \(\mathcal O(2^2n^3)\)，记得模一个大质数。

#include<bits/stdc++.h>
using namespace std;
const int N=80,M=(1<<12)+5,mod=998244353;
int n,m,w,x,y,v,a[N][N][M],tmp[N][N],b[M],inv=(mod+1)/2;
char s[N]; 
void FWT(int *f,int n,int op){
	for(int k=2,lg=0;k<=n;k<<=1,lg++)
		for(int i=0,m=k>>1;i<n;i+=k)
			for(int j=i;j<i+m;j++){
				if(s[lg]=='|') f[j+m]=(f[j+m]+f[j]*op)%mod;
				else if(s[lg]=='&') f[j]=(f[j]+f[j+m]*op)%mod;
				else{
					x=f[j],y=f[j+m],f[j]=(x+y)%mod,f[j+m]=(x-y+mod)%mod;
					if(op==-1) f[j]=1ll*f[j]*inv%mod,f[j+m]=1ll*f[j+m]*inv%mod;
				} 
			}
	for(int i=0;i<n;i++) f[i]=(f[i]+mod)%mod;
}
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
int det(int a[N][N]){
	int ans=1;
	for(int i=1;i<n;i++){
		int p=i;
		for(int j=i;j<n;j++) if(a[j][i]){p=j;break;}
		if(p!=i) swap(a[i],a[p]),ans=mod-ans;
		for(int j=i+1,iv=qpow(a[i][i],mod-2);j<n;j++){
			int x=1ll*a[j][i]*iv%mod;
			for(int k=i;k<n;k++) a[j][k]=(a[j][k]+1ll*(mod-a[i][k])*x%mod)%mod;
		}
		ans=1ll*ans*a[i][i]%mod;
	}
	return ans;
}
signed main(){
	scanf("%d%d%s",&n,&m,s),w=strlen(s);
	for(int i=1;i<=m;i++){
		scanf("%d%d%d",&x,&y,&v);
		a[x][x][v]++,a[y][y][v]++,a[x][y][v]--,a[y][x][v]--;
		//a[x][y][v] 实际上是矩阵 (x,y) 对应幂级数 x^v 前的系数
	} 
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++) FWT(a[i][j],1<<w,1); 
	for(int x=0;x<(1<<w);x++){
		for(int i=1;i<=n;i++)
			for(int j=1;j<=n;j++) tmp[i][j]=a[i][j][x];
		b[x]=det(tmp);
	}
	FWT(b,1<<w,-1);
	for(int i=(1<<w)-1;i>=0;i--) if(b[i]) printf("%d\n",i),exit(0);
	puts("-1");
	return 0;
}

5. 图

2022.7.30 组合意义 + 矩阵树定理

对于一张无向图 \(G\)，定义它的权值 \(f(G)\)：

1. 若 \(G\) 的边数 \(>\) 点数，则 \(f(G)=0\)。
2. 设 \(G\) 有 \(n\) 个点 \(k\) 条边，若不存在一种加入 \(n-k\) 条边的方案使图连通，则同样 \(f(G)=0\)。
3. 否则，\(f(G)\) 为这张图所有连通块大小的乘积。

现给出一张 \(n\) 个点 \(m\) 条边的无重边无自环的无向图，求共 \(2^m\) 个边集的子集的导出子图的权值之和 \(\bmod 998244353\)。

\(1\leq n\leq 16\)。

\(f(G)\neq 0\)，当且仅当：至多一个连通块是基环树，其余连通块是树。

考虑组合意义：连通块大小的乘积，相当于从每个连通块里选一个点的方案数。

第一种情况：所有连通块是树。

• 答案为 \((T,S)\) 的数量，其中 \(T\) 是原图的生成森林，\(S\) 表示每个连通块选出的点的集合。

• 构造辅助图 \(G'\)：建一个虚点，与所有点之间连一条边。

  发现 \((T,S)\) 与 \(G'\) 的生成树之间存在一一映射（一个点选入 \(S\)，相当于把它与虚点之间的边选入 \(G'\) 的生成树），因此答案等于 \(G'\) 的生成树数量。矩阵树定理即可。

  时间复杂度 \(\mathcal O(n^3)\)。

第二种情况：一个连通块是基环树，其余连通块是树。

• 枚举环的点集 \(S\)，将环缩成一个点后，转化为第一种情况。

  注意，虚点要向缩的点连 \(|S|\) 条边。贡献为缩完后的答案 \(\times\) 点集 \(S\) 连成环的方案数。

  时间复杂度 \(\mathcal O(2^nn^3)\)。

#include<bits/stdc++.h>
using namespace std;
const int N=20,M=310,S=1<<16,mod=998244353;
int n,m,u[M],v[M],e[N][N],a[N][N],cir[S],f[N][S],tot,id[N],ans;
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
int det(int n){
	int ans=1;
	for(int i=1;i<=n;i++){
		int p=i;
		for(int j=i;j<=n;j++) if(a[j][i]) p=j;
		if(p!=i) swap(a[i],a[p]),ans=mod-ans;
		for(int j=i+1,iv=qpow(a[i][i],mod-2);j<=n;j++){
			int x=mod-1ll*a[j][i]*iv%mod;
			for(int k=1;k<=n;k++) a[j][k]=(a[j][k]+1ll*a[i][k]*x%mod)%mod;
		}
		ans=1ll*ans*a[i][i]%mod;
	}
	return ans;
}
int calc(int s){	//环为 s
	tot=s!=0,memset(a,0,sizeof(a));
	for(int i=1;i<=n;i++) id[i]=s>>(i-1)&1?1:++tot;
	auto add=[&](int x,int y){
		if(x!=y) a[x][x]++,a[y][y]++,(a[x][y]+=mod-1)%=mod,(a[y][x]+=mod-1)%=mod;	//注意 x!=y
	}; 
	for(int i=1;i<=m;i++) add(id[u[i]],id[v[i]]);
	for(int i=1;i<=n;i++) add(id[i],tot+1);	//连向虚点
	return det(tot);
}
signed main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=m;i++)
		scanf("%d%d",&u[i],&v[i]),e[u[i]][v[i]]=e[v[i]][u[i]]=1;
	for(int p=1;p<=n;p++){	//枚举环的起点
		memset(f,0,sizeof(f)),f[p][1<<(p-1)]=1;
		for(int s=0;s<(1<<n);s++)
			for(int i=1;i<=n;i++) if((s>>(i-1)&1)&&f[i][s]){
				if(e[i][p]) (cir[s]+=f[i][s])%=mod;
				for(int j=1;j<=n;j++)
					if(!(s>>(j-1)&1)&&e[i][j]) (f[j][s|(1<<(j-1))]+=f[i][s])%=mod;
			}
	}
	for(int s=1;s<(1<<n);s++){	//枚举环
		int sz=__builtin_popcount(s);
		if(sz>2) ans=(ans+1ll*calc(s)*cir[s]%mod*qpow(sz*2,mod-2)%mod)%mod;	//一个大小为 sz 的环会被算 sz*2 次（环上每个点都能作为起点，有两种方向）
	}
	printf("%d\n",(ans+calc(0))%mod);	//加上没有环的情况
	return 0;
}

也有一个 \(\mathcal O(3^n)\) 做法。

求出 \(f_{S}\) 表示将 \(S\) 里的点连成一棵树的方案数。

怎么避免算重？记 \(x=lowbit(S)\)，一般来说都是 \(f_{S}\gets f_{S-\{x\}}\) 来避免算重，但问题是 \(S-\{x\}\) 不一定连通。

考虑 \(S\) 形成的树把 \(x\) 删掉会裂成若干子树，记 \(y=lowbit(S-\{x\})\)，枚举 \(y\) 所在的子树的点集 \(T\)，\(f_S\gets f_T\times f_{S-T}\times |N(x)\cap T|\)。

类似地，设 \(g_S\) 表示将 \(S\) 里的点连成一棵基环树的方案数。记 \(x=lowbit(S)\)，\(y=lowbit(S-\{x\})\)，\(g_S\gets (f_T\times g_{S-T}+g_T\times f_{S-T})\times |N(x)\cap T|+f_T\times f_{S-T}\times \binom{|N(x)\cap T|}{2}\)。

6. LOJ#6044. 「雅礼集训 2017 Day8」共

矩阵树定理 + 手推行列式

求有多少棵 \(n\) 个节点的树（以 \(1\) 为根），恰好有 \(k\) 个深度为奇数的点，答案对 \(p\) 取模。

\(1<k<n\leq 5\times 10^5\)，\(p\) 为质数。

首先有个小套路，将整棵树按深度奇偶性转化为一张二分图，显然我们只能在左右部点之间连边。

\(1\) 只能在左部，组成左部的方案数为 \(\large\binom{n-1}{k-1}\)。不妨设左部的另外 \(k-1\) 个点为 \([2,k]\)，问题转化为，有一张二分图，左部有 \(k\) 个点，右部有 \(n-k\) 个点，要在它们之间连 \(n-1\) 条边使其构成一棵生成树。

考虑矩阵树定理（也可以直接 Prufer 序列），由于 \(n\) 较大，不能暴力求。注意到这张图很特殊，考虑手推行列式。基尔霍夫矩阵（度数矩阵 - 邻接矩阵）：

\[\begin{bmatrix} n-k&0&\cdots&0&-1&-1&\cdots&-1\\ 0&n-k&\cdots&0&-1&-1&\cdots&-1\\ 0&0&\cdots&0&-1&-1&\cdots&-1\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&n-k&-1&-1&\cdots&-1\\ -1&-1&\cdots&-1&k&0&\cdots&0\\ -1&-1&\cdots&-1&0&k&\cdots&0\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\ -1&-1&\cdots&-1&0&0&\cdots&k\\ \end{bmatrix} \]

其中左边有 \(k\) 列，右边有 \(n-k\) 列。去掉第一行第一列后，左边有 \(k-1\) 列，右边有 \(n-k\) 列（上面有 \(k-1\) 行，下面有 \(n-k\) 行）。考虑倍加相消，从第 \(k\) 行开始，自上而下地依次让每一行加上下一行的 \(-1\) 倍。

\[\begin{bmatrix} n-k&0&\cdots&0&-1&-1&-1&\cdots&-1&-1\\ 0&n-k&\cdots&0&-1&-1&-1&\cdots&-1&-1\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&n-k&-1&-1&-1&\cdots&-1&-1\\ 0&0&\cdots&0&k&-k&0&\cdots&0&0\\ 0&0&\cdots&0&0&k&-k&\cdots&0&0\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&0&0&0&0&\cdots&k&-k\\ -1&-1&\cdots&-1&0&0&0&\cdots&0&k\\ \end{bmatrix} \]

然后依次将前 \(k-1\) 行的 \(\large\frac{1}{n-k}\) 倍加到最后一行：

\[\begin{bmatrix} n-k&0&\cdots&0&-1&-1&-1&\cdots&-1&-1\\ 0&n-k&\cdots&0&-1&-1&-1&\cdots&-1&-1\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&n-k&-1&-1&-1&\cdots&-1&-1\\ 0&0&\cdots&0&k&-k&0&\cdots&0&0\\ 0&0&\cdots&0&0&k&-k&\cdots&0&0\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&0&0&0&0&\cdots&k&-k\\ 0&0&\cdots&0&-\frac{k-1}{n-k}&-\frac{k-1}{n-k}&-\frac{k-1}{n-k}&\cdots&-\frac{k-1}{n-k}&k-\frac{k-1}{n-k}\\ \end{bmatrix} \]

记 \(A=k,B=\large\frac{k-1}{n-k}\)：

\[\begin{bmatrix} n-k&0&\cdots&0&-1&-1&-1&\cdots&-1&-1\\ 0&n-k&\cdots&0&-1&-1&-1&\cdots&-1&-1\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&n-k&-1&-1&-1&\cdots&-1&-1\\ 0&0&\cdots&0&A&-A&0&\cdots&0&0\\ 0&0&\cdots&0&0&A&-A&\cdots&0&0\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&0&0&0&0&\cdots&A&-A\\ 0&0&\cdots&0&-B&-B&-B&\cdots&-B&A-B\\ \end{bmatrix} \]

然后从第 \(k\) 行开始到倒数第二行，自上而下地依次将第 \(i\) 行的 \((i-k+1)\times \large\frac{B}{A}\) 倍加到最后一行：

\[\begin{bmatrix} n-k&0&\cdots&0&-1&-1&-1&\cdots&-1&-1\\ 0&n-k&\cdots&0&-1&-1&-1&\cdots&-1&-1\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&n-k&-1&-1&-1&\cdots&-1&-1\\ 0&0&\cdots&0&A&-A&0&\cdots&0&0\\ 0&0&\cdots&0&0&A&-A&\cdots&0&0\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&0&0&0&0&\cdots&A&-A\\ 0&0&\cdots&0&0&0&0&\cdots&0&A-(n-k)B\\ \end{bmatrix} \]

我们发现，\(A-(n-k)B=1\)，所以该矩阵的行列式就是 \((n-k)^{k-1}\times A^{n-k-1}=(n-k)^{k-1}\times k^{n-k-1}\)。\(ans=(n-k)^{k-1}\times k^{n-k-1}\times \large\binom{n-1}{k-1}\)。

#include<bits/stdc++.h>
using namespace std;
const int N=5e5+5;
int n,k,mod,fac[N],inv[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
int C(int n,int m){return n<m?0:1ll*fac[n]*inv[m]%mod*inv[n-m]%mod;}
signed main(){
	scanf("%d%d%d",&n,&k,&mod),fac[0]=inv[0]=inv[1]=1;
	for(int i=2;i<=n;i++) inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
	for(int i=1;i<=n;i++)
		fac[i]=1ll*fac[i-1]*i%mod,inv[i]=1ll*inv[i-1]*inv[i]%mod;
	printf("%lld\n",1ll*qpow(n-k,k-1)*qpow(k,n-k-1)%mod*C(n-1,k-1)%mod);
	return 0; 
} 

四、代数余子式、伴随矩阵

余子式：对于一个矩阵 \(A=(a_{i,j})_{n\times n}\)，\(a_{i,j}\) 的余子式 \(M_{i,j}\) 为矩阵 \(a\) 删去第 \(i\) 行和第 \(j\) 列后的行列式。

代数余子式：\(a_{i,j}\) 的代数余子式为 \(A_{i,j}=(-1)^{i+j}M_{i,j}\)。

一个矩阵的行列式是它任意一行或任意一列的 \(\sum a_{i,j}A_{i,j}\)。也就是所谓的对某一行/列展开。

伴随矩阵：将 \(a_{i,j}\) 换成 \(A_{i,j}\) 再转置，记作 \(A^{*}\)。

\[\begin{bmatrix} A_{1,1}&A_{2,1}&\cdots&A_{n,1}\\ A_{1,2}&A_{2,2}&\cdots&A_{n,2}\\ \vdots&\vdots&\ddots&\vdots\\ A_{1,n}&A_{2,n}&\cdots&A_{n,n} \end{bmatrix} \]

一个公式：\(AA^*=A^*A=\text{det}(A)I\)。那么可以得到 \(A^*=\frac{\text{det}(A)I}{A}\)（逆矩阵每个位置的值 \(\times \text{det}(A)\)），可以矩阵求逆算。而算出伴随矩阵就能知道代数余子式（记得转置！）。

矩阵求逆：\(AA^{-1}=I\)。初始时左边放 \(A\) 右边放 \(I\)，放在一起消。目标是将左边变成 \(I\)，消的时候右边和左边做一样的操作。当左边变成 \(I\) 时右边就是 \(A^{-1}\)。若左边无法变成 \(I\) 则无解。

1. CF736D Permutations（*2800）

2022.8.18

给出一张两侧各 \(n\) 个点，共 \(m\) 条边的二分图，保证其完美匹配个数为奇数。对于二分图的每一条边，询问将这条边删去后，剩下的二分图的完美匹配数是否仍为奇数。

\(1\leq n\leq 2000\)，\(n\leq m\leq \min(n^2,5\times 10^5)\)。

令 \(a_{i,j}=[(i,j')\in E]\)，那么二分图完美匹配数为奇数等价于 \(\sum_ p\prod_{i=1}^n a_{i,p_i}\equiv 1\pmod 2\)，由于 \(-1\equiv 1\pmod 2\)，所以也等价于 \(\text{det}(A)\equiv 1\pmod 2\)。

删除一条边后的完美匹配数，就是强制这条边不在完美匹配中。容斥转化为强制这条边在完美匹配中，这恰好是 \(a_{i,j}\) 的余子式 \(M_{i,j}\)。我们只关心 \(M_{i,j}\) 的奇偶性，同样由于 \(-1\equiv 1\pmod 2\)，所以只需计算 \(a_{i,j}\) 的代数余子式 \(A_{i,j}=(-1)^{i+j}M_{i,j}\)。

求出 \(A\) 的伴随矩阵 \(A^*\) 就能求出代数余子式，而 \(AA^*=A^*A=\text{det}(A)I\)，因而 \(A^*=\frac{\text{det}(A)I}{A}\)，矩阵求逆即可。由于运算在 \(\mathbb F_2\) 下进行，可以 bitset 优化。

时间复杂度 \(\mathcal O(\frac{n^3}{w})\)。

#include<bits/stdc++.h>
using namespace std;
const int N=2e3+5,M=5e5+5;
int n,m,x[M],y[M];
bitset<N<<1>a[N];
signed main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=m;i++)
		scanf("%d%d",&x[i],&y[i]),a[x[i]][y[i]]=1;
	for(int i=1;i<=n;i++) a[i][i+n]=1;	//矩阵求逆
	for(int i=1;i<=n;i++){
		if(!a[i][i])
			for(int j=i+1;j<=n;j++)
				if(a[j][i]){swap(a[i],a[j]);break;}
		for(int j=1;j<=n;j++)
			if(j!=i&&a[j][i]) a[j]^=a[i];
	}
	for(int i=1;i<=m;i++) puts(a[y[i]][x[i]+n]?"NO":"YES");	//记得转置。强制不在 = 总 - 强制在，若 a[y[i]][x[i]+n]!=0 则强制在的方案数为奇数，所以强制不在的方案数为偶数
	return 0;
} 

2. 生成树计数

2022.8.18 生日攻击 + 线性代数

给出 \(k\)，要求构造一个点数 \(\leq 100\) 的有向图和无向图，不能存在重边 和自环（有向图中 \((u,v)\) 和 \((v,u)\) 不算重边），使其以 \(n\) 为根的生成树（有向图中则是内向生成树）个数是 \(k\) 的倍数（不能为 \(0\)）。无解输出 \(-1\)。

\(1\leq k<2^{31}\)。

考虑一张无向连通图 \(G\)，其中 \((a,b)\in E\)，记 \(f(G)\) 为 \(G\) 中强制选择 \((a,b)\) 的生成树个数，\(g(G)\) 为 \(G\) 中强制不选 \((a,b)\) 的生成树个数（假定 \(f(G),g(G)\neq 0\)）。

对于两张图 \(G_1,G_2\)，删掉它们中的边 \((a,b)\)，并令 \(a,b\) 分别重合，则合并后新图的生成树个数为 \(f(G_1)g(G_2)+f(G_2)g(G_1)\)，我们希望它是 \(k\) 的倍数。

设 \(k\) 是质数，且 \(f(G_1),f(G_2),g(G_1),g(G_2)\) 都不是 \(k\) 的倍数，则 \(f(G_1)g(G_2)\equiv -f(G_2)g(G_1)\pmod k\)，\(\frac{f(G_1)}{g(G_1)}\equiv -\frac{f(G_2)}{g(G_2)}\pmod k\)，因此使用生日悖论冲撞 \(\frac{f(G)}{g(G)}\) 即可。图的大小只需要让 \(n^{n-2}\) 一定程度上超过 \(k\)。

注意到一次高斯消元才能求出一个生成树个数很浪费。考虑每次取基尔霍夫矩阵的前 \(n-1\) 行 \(n-1\) 列做行列式，删掉一条 \((x,n)\) 的边只会把 \(a_{x,x}\) 减 \(1\)。对第 \(x\) 行展开，由于 \(\text{det}(A)=\sum_{i=1}^{n-1}a_{x,i}A_{x,i}\)，而只有 \(a_{x,x}\) 减了 \(1\)，所以 \(\text{det}(A)\) 减了 \(A_{x,x}\)。要算出所有 \(A_{x,x}\)，所以只要根据 \(A^*=\text{det}(A)I\cdot A^{-1}\) 算出伴随矩阵即可。

#include<bits/stdc++.h>
using namespace std;
const int N=110;
int k,n,m,a[N][N],b[N][N<<1],mod,d,to[N];
struct G{
	int n;
	vector<int>v[N]; 
	G friend operator+(G a,G b){	//直接合并
		G c; c=a,c.n+=b.n-1;
		for(int i=1;i<=b.n;i++)
			for(int j:b.v[i]) c.v[i+a.n-1].push_back(j+a.n-1);
		return c;
	}
}tmp,ans;
struct E{
	int ban,v[N];
	G friend operator+(E a,E b){	//把边 a.ban 和 b.ban 都删掉，并令边 a.ban 和 b.ban 的端点分别重合，合并成一张新的图
		G ans; int tot=n;
		for(int i=1;i<=n;i++){
			to[i]=i==b.ban?a.ban:(i==n?n:++tot);
			for(int j=1;j<i;j++)
				if((a.v[i]>>(j-1)&1)&&!(i==n&&j==a.ban)) ans.v[i].push_back(j),ans.v[j].push_back(i);
		}
		ans.n=tot;
		for(int i=1;i<=n;i++)
			for(int j=1;j<i;j++)
				if((b.v[i]>>(j-1)&1)&&!(i==n&&j==b.ban)) ans.v[to[i]].push_back(to[j]),ans.v[to[j]].push_back(to[i]);
		return ans;
	}
}e;
map<int,E>mp;
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
int det(int n){	//计算行列式的同时计算逆矩阵
	int ans=1;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			b[i][j]=(0ll+a[i][j]%mod+mod)%mod,b[i][j+n]=i==j;	//注意 0ll+，因为 mod 可能很大，直接加会爆 int（调了半天/kk）。右半边拿来算逆矩阵
	for(int i=1;i<=n;i++){
		int p=i;
		for(int j=i;j<=n;j++) if(b[j][i]){p=j;break;}
		if(p^i) swap(b[i],b[p]),ans=mod-ans;
		ans=1ll*ans*b[i][i]%mod;	//要在 b[i][i] 消成 1 前算行列式
		for(int j=i,iv=qpow(b[i][i],mod-2);j<=2*n;j++) b[i][j]=1ll*b[i][j]*iv%mod;	//注意要消成 I！
		for(int j=1;j<=n;j++) if(i^j){
			int x=mod-b[j][i];
			for(int k=i;k<=2*n;k++) b[j][k]=(b[j][k]+1ll*b[i][k]*x%mod)%mod;
		}
	}
	return ans;
}
void add(int x,int y){a[x][x]++,a[y][y]++,a[x][y]--,a[y][x]--;}
void rnd(int n){	//随机一张图
	for(int i=1;i<=n;i++) fill(a[i]+1,a[i]+1+n,0);
	for(int i=2;i<n;i++) add(rand()%(i-1)+1,i);
	for(int i=1;i<n;i++)
		for(int j=1;j<i;j++)
			if(!a[i][j]&&rand()%3==0) add(i,j);
	for(int i=1;i<n;i++) add(i,n);	//保证生成树个数不为 0
}
G calc(int p){
	mp.clear(),mod=p,n=2;
	if(mod<5) n=4;
	else{while(pow(n,n-2)<mod) n++; n+=2;}
	G ans; ans.n=n;
	while(1){
		rnd(n);
		if(!(d=det(n-1))){
			for(int i=1;i<=n;i++)
				for(int j=1;j<i;j++)
					if(a[i][j]) ans.v[i].push_back(j),ans.v[j].push_back(i);
			return ans;
		}
		for(int i=1;i<=n;i++){
			e.v[i]=0;
			for(int j=1;j<=n;j++)
				if(i^j&&a[i][j]) e.v[i]|=1<<(j-1);
		}
		for(int x=1;x<n;x++){
			int g=(d-1ll*d*b[x][x+(n-1)]%mod+mod)%mod,f=(0ll+d-g+mod)%mod;	//注意 0ll+。用 A*=det(A)*A^{-1} 计算 a_{x,x} 的代数余子式，即 d*b[x][x+(n-1)]。g 表示强制不选 (x,n)，f 表示强制选 (x,n)
			if(!g){
				for(int i=1;i<=n;i++)
					for(int j=1;j<i;j++)
						if(a[i][j]&&!(i==n&&j==x)) ans.v[i].push_back(j),ans.v[j].push_back(i);
				return ans;
			}
			int u=1ll*f*qpow(g,mod-2)%mod,v=(mod-u)%mod;	//生日悖论冲撞 f/g
			e.ban=x,mp[u]=e;
			if(mp.find(v)!=mp.end()) return e+mp[v];
		}
	}
}
signed main(){
	scanf("%d",&k),ans.n=1;
	for(int i=2;1ll*i*i<=k;i++) if(k%i==0){	//质因数分解转化成 k 为质数的情况
		tmp=calc(i);
		while(k%i==0) ans=ans+tmp,k/=i;
	}
	if(k>1) ans=ans+calc(k);
	for(int i=1;i<=ans.n;i++) m+=ans.v[i].size();
	printf("%d %d\n",ans.n,m);
	for(int i=1;i<=ans.n;i++)
		for(int j:ans.v[i]) printf("%d %d\n",i,j);
	printf("%d %d\n",ans.n,m/2);
	for(int i=1;i<=ans.n;i++)
		for(int j:ans.v[i])
			if(i<j) printf("%d %d\n",i,j);
	return 0;
}

3. 树形图求和

2022.9.8

给出一张 \(n\) 个点 \(m\) 条边的有向图，求其所有以 \(n\) 为根的内向图的边权和 \(\bmod 10^9+7\)。

\(2\leq n\leq 300\)，\(0\leq m\leq 10^5\)，\(1\leq w_i\leq 10^9\)，2s。

考虑每条边的贡献，要算它出现在多少树形图中。容斥转成不在，算树形图总数 - 去掉这条边后的树形图个数。

按上一题的方法，每次取基尔霍夫矩阵的前 \(n-1\) 行 \(n-1\) 列做行列式，删掉一条边 \(x\to y\) 会把 \(a_{x,x}\) 减 \(1\)，\(a_{x,y}\) 加 \(1\)。对第 \(x\) 行展开，由于 \(\text{det}(A)=\sum_{i=1}^{n-1}a_{x,i}A_{x,i}\)，所以 \(\text{det}(A)\) 减了 \(A_{x,x}\)，加了 \(A_{x,y}\)。

要算出所有 \(A_{x,x},A_{x,y}\)，所以只要根据 \(A^*=\text{det}(A)I\cdot A^{-1}\) 算出伴随矩阵即可。

#include<bits/stdc++.h>
using namespace std;
const int N=310,M=1e5+5,mod=1e9+7;
int n,m,u[M],v[M],w[M],a[N][N],b[N][N<<1],f[N][N],d,ans;
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
int det(int n){
	int ans=1;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			b[i][j]=(a[i][j]+mod)%mod,b[i][j+n]=i==j;
	for(int i=1;i<=n;i++){
		int p=i;
		for(int j=i;j<=n;j++) if(b[j][i]){p=j;break;}
		if(p^i) swap(b[i],b[p]),ans=mod-ans;
		ans=1ll*ans*b[i][i]%mod;
		for(int j=i,iv=qpow(b[i][i],mod-2);j<=2*n;j++) b[i][j]=1ll*b[i][j]*iv%mod;
		for(int j=1;j<=n;j++) if(i^j){
			int x=mod-b[j][i];
			for(int k=i;k<=2*n;k++) b[j][k]=(b[j][k]+1ll*b[i][k]*x%mod)%mod;
		}
	}
	return ans;
}
signed main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=m;i++)
		scanf("%d%d%d",&u[i],&v[i],&w[i]),a[u[i]][u[i]]++,a[u[i]][v[i]]--;
	d=det(n-1);
	for(int i=1;i<n;i++)
		for(int j=1;j<n;j++) f[i][j]=1ll*d*b[i][j+(n-1)]%mod;
	for(int i=1;i<=m;i++)
		ans=(ans+1ll*w[i]*(d-(d-f[u[i]][u[i]]+f[v[i]][u[i]])%mod)%mod)%mod;
	printf("%d\n",(ans+mod)%mod);
	return 0;
}

五、BEST 定理

目的：

• 求有向图欧拉回路。注意不是无向图。

使用前提：

• 是欧拉图，即存在欧拉回路。check 每个点“入度 = 出度”且所有边弱连通。

  需要提前特判掉不是欧拉图的情况。

定理内容：

• 如果需要去除边序列循环同构的欧拉回路，\(ans=Tree(S)\times \prod_{i=1}^n(out_i-1)!\)。

  如果钦定 \(S\) 为起点（钦定边序列以 \(S\) 的出边为开头），前一种情况的欧拉回路的边序列，把每个 \(S\) 的出边循环位移到第一个位置，就对应了 \(out_S\) 种这种情况的边序列。\(ans=Tree(S)\times \prod_{i=1}^n(out_i-1)!\times out_S\)。

  其中，\(Tree(S)\) 表示以 \(S\) 为根的 内向树/外向树 个数，\(out_i\) 为 \(i\) 的 出度/入度（因为出度 = 入度）。

细节：

• 注意计算行列式时忽略孤立点（不然把 \(0\) 乘进去了）。

感性理解：

• 假如钦定了 \(S\) 为起点。考虑对于每个点 \(x\neq S\)，提出它出边中走的最后一条，这些边构成一棵以 \(S\) 为根的内向树，其他非树边走的时间顺序任意排列。

  注意由于是有向图，出边和入边是不一样的，\(S\) 的所有出边都是非树边。

注意计算行列式时忽略孤立点（不然把 \(0\) 乘进去了）。

1. P5807 Which Dreamed It /【模板】BEST 定理

求从 \(1\) 出发的有向图欧拉回路个数 \(\bmod 10^6+3\)。

\(1\leq T\leq 15\)，\(1\leq n\leq 100\)，\(0\leq \sum m\leq 2\times 10^5\)。

#include<bits/stdc++.h>
using namespace std;
const int N=110,M=2e5+5,mod=1e6+3;
int t,n,x,a[N][N],f[N],in[N],out[N],fac[M],flg,ans;
int find(int x){return x==f[x]?x:f[x]=find(f[x]);}
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
} 
int det(){
	int ans=1;
	for(int i=2;i<=n;i++){
		int p=i;
		for(int j=i;j<=n;j++) if(a[j][i]){p=j;break;}
		if(i^p) ans=mod-ans,swap(a[i],a[p]);
		for(int j=i+1,iv=qpow(a[i][i],mod-2);j<=n;j++) if(a[j][i]){
			int x=mod-1ll*a[j][i]*iv%mod;
			for(int k=i;k<=n;k++)
				a[j][k]=(a[j][k]+1ll*a[i][k]*x%mod)%mod;
		}
		if(out[i]) ans=1ll*ans*a[i][i]%mod;
	}
	return ans;
}
signed main(){
	scanf("%d",&t),fac[0]=1;
	for(int i=1;i<M;i++) fac[i]=1ll*fac[i-1]*i%mod;
	while(t--){
		scanf("%d",&n),flg=1;
		for(int i=1;i<=n;i++) f[i]=i,fill(a[i]+1,a[i]+1+n,0),in[i]=0;
		for(int i=1;i<=n;i++){
			scanf("%d",&out[i]);
			for(int j=1;j<=out[i];j++)
				scanf("%d",&x),f[find(i)]=find(x),a[x][x]++,(a[i][x]+=mod-1)%=mod,in[x]++;
		}
		for(int i=1;i<=n;i++)
			if(in[i]^out[i]||(out[i]&&find(i)^find(1))){puts("0");goto End;}
		for(int i=1;i<=n;i++) flg&=!out[i];
		if(flg){puts("1");goto End;}
		ans=det();
		for(int i=1;i<=n;i++)
			if(out[i]) ans=1ll*ans*fac[out[i]-1]%mod;
		printf("%lld\n",1ll*ans*out[1]%mod); 
		End: ;
	}
	return 0;
}

六、LGV 引理

• 给出一个 DAG，设 \(w(P)\) 为有向路径 \(P\) 上所有边权的乘积，\(f(a,b)\) 为 \(a\to b\) 的所有有向路径边权乘积之和，即 \(f(a,b)=\sum_{P:a\to b}w(P)\)。

  给出一个起点集合 \(A\) 和终点集合 \(B\)。设 \(p\) 为某个 \(1\sim n\) 的排列，\(P_i\) 为 \(a_i\to b_{p_i}\) 的一条路径，且 \(P_{1\sim n}\) 两两互不相交，记 \(P_{1\sim n}\) 为 \(A\to B\) 的一个不交路径集合。则满足：

  \[M=\begin{bmatrix} f(a_1,b_1)&f(a_1,b_2)&\cdots&f(a_1,b_n)\\ f(a_2,b_1)&f(a_2,b_2)&\cdots&f(a_2,b_n)\\ \vdots&\vdots&\ddots&\vdots\\ f(a_n,b_1)&f(a_n,b_2)&\cdots&f(a_n,b_n) \end{bmatrix} \]
  \[\text{det}(M)=\sum_{P:A\to B}(-1)^{\tau(p)}\prod_{i=1}^n w(P_i) \]

  即 \(\text{det}(M)\) 为所有 \(A\to B\) 不相交路径 的 带符号和。

• 路径不交问题：

  求从每个 \(a_i\to b_i\)，所有路径不交方案 \(P_{1\sim n}\) 的 \(\prod_{i=1}^n w(P_i)\) 之和。

  根据行列式的计算式：

  \[\text{det}(M)=\sum_p(-1)^{\tau(p)}\prod_{i=1}^n f(a_i,b_{p_i}) \]

  若两个人的路径相交，可以将相交后的部分交换，相当于两人分别走到对方的终点，逆序对数 \(-1\)。本质是容斥，\(p\) 表示每个人走到了谁的终点去，容斥系数是 \((-1)^{\tau(p)}\)，最后只会把没有逆序对的方案算进答案，答案就是 \(\text{det}(M)\)。

1. P6657 【模板】LGV 引理

给出一张 \(n\times n\) 的棋盘，每次只能向右或向下走一步。

有 \(m\) 个棋子，分别要从 \((a_i,1)\) 走到 \((b_i,n)\)，求路径不交的方案数 \(\bmod 998244353\)。

\(T\leq 5\)，\(2\leq n\leq 10^6\)，\(1\leq m\leq 100\)，\(1\leq a_1\leq a_2\leq \cdots\leq a_m\leq n\)，\(1\leq b_1\leq b_2\cdots\leq b_m\leq n\)。

每条路径边权都是 \(1\)，\(f(a_i,b_j)\) 实际上就是 \((a_i,1)\to (b_j,n)\) 的方案数 \(\binom{b_j-a_i+n-1}{n-1}\)。

#include<bits/stdc++.h>
using namespace std;
const int N=2e6+5,M=110,mod=998244353;
int t,n=2e6,m,a[N],b[N],fac[N],inv[N],e[M][M];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
int C(int n,int m){return n<m?0:1ll*fac[n]*inv[m]%mod*inv[n-m]%mod;}
int det(){
	int ans=1;
	for(int i=1;i<=m;i++){
		int p=i;
		for(int j=i;j<=m;j++) if(e[j][i]){p=j;break;}
		if(p^i) swap(e[i],e[p]),ans=mod-ans;
		for(int j=i+1,iv=qpow(e[i][i],mod-2);j<=m;j++){
			int x=mod-1ll*e[j][i]*iv%mod;
			for(int k=i;k<=m;k++) e[j][k]=(e[j][k]+1ll*e[i][k]*x%mod)%mod;
		}
		ans=1ll*ans*e[i][i]%mod;
	}
	return ans;
}
signed main(){
	scanf("%d",&t),fac[0]=inv[0]=inv[1]=1;
	for(int i=2;i<=n;i++) inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
	for(int i=1;i<=n;i++)
		fac[i]=1ll*fac[i-1]*i%mod,inv[i]=1ll*inv[i-1]*inv[i]%mod;
	while(t--){
		scanf("%d%d",&n,&m);
		for(int i=1;i<=m;i++) scanf("%d%d",&a[i],&b[i]);
		for(int i=1;i<=m;i++)
			for(int j=1;j<=m;j++) e[i][j]=C(b[j]-a[i]+n-1,n-1);
		printf("%d\n",det());
	}
	return 0;
}

2. CF348D Turtles（*2500）

给出一张 \(n\times m\) 的网格图，有一些障碍格子。求有多少对 \((1,1)\to (n,m)\) 的路径，满足除起点和终点外不相交。对 \(10^9+7\) 取模。

\(2\leq n,m\leq 3000\)。

LGV 引理要求的不相交是包括起点和终点的，所以 \(a_1=a_2=(1,1)\)，\(b_1=b_2=(n,m)\) 算出来的答案是 \(0\)。

可以 \(a_1=(1,2)\)，\(a_2=(2,1)\)，\(b_1=(n-1,m)\)，\(b_2=(n,m-1)\)。

设 \(f1_{i,j},f2_{i,j}\) 分别表示从 \((1,2),(2,1)\) 走到 \((i,j)\) 的方案数，\(ans=\text{det}\begin{bmatrix}f1_{n-1,m}&f1_{n,m-1}\\f2_{n-1,m}&f2_{n,m-1}\end{bmatrix}\)，直接算。

3. P7736 [NOI2021] 路径交点

有一张 \(k\) 层的图，第 \(i\) 层有 \(n_i\) 个点，相邻层之间有一些有向边（从第 \(i\) 层连向 \(i+1\) 层）。每层的点从上到下排列，层从左到右排列。

以第 \(1\) 层的 \(n_1\) 个点作为起点，第 \(k\) 层的 \(n_k\) 个点作为终点（保证 \(n_1=n_k\)），将起点和终点一一匹配，求所有路径不交方案中，交点数量为偶数的减去为奇数的方案有多少个 。对 \(98244353\) 取模。

\(t\leq 5\)，\(2\leq k\leq 100\)，\(2\leq n_1\leq 100\)，\(n_1\leq n_i\leq 2\times n_1\)。

发现对于一组路径方案 \(P\)，交点个数的奇偶性 = \(\tau(p)\) 的奇偶性。

所以答案就是 \(\text{det}(M)\)。\(f(a_i,b_j)\) 可以矩阵快速幂或者直接 bfs。

#include<bits/stdc++.h>
using namespace std;
const int N=210,mod=998244353;
int t,k,n[N],m[N],x,y;
struct mat{
	int n,m,x[N][N];
	friend mat operator*(mat a,mat b){
		mat c; c.n=a.n,c.m=b.m;
		for(int i=1;i<=c.n;i++)
			for(int j=1;j<=c.m;j++){
				c.x[i][j]=0;
				for(int k=1;k<=a.m;k++) c.x[i][j]=(c.x[i][j]+1ll*a.x[i][k]*b.x[k][j]%mod)%mod;
			}
		return c;
	} 
}tmp,ans;
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
int det(int e[N][N],int n){
	int ans=1;
	for(int i=1;i<=n;i++){
		int p=i;
		for(int j=i;j<=n;j++) if(e[j][i]){p=j;break;}
		if(p^i) swap(e[i],e[p]),ans=mod-ans;
		for(int j=i+1,iv=qpow(e[i][i],mod-2);j<=n;j++){
			int x=mod-1ll*e[j][i]*iv%mod;
			for(int k=i;k<=n;k++) e[j][k]=(e[j][k]+1ll*e[i][k]*x%mod)%mod;
		}
		ans=1ll*ans*e[i][i]%mod;
	}
	return ans;
}
signed main(){
	scanf("%d",&t);
	while(t--){
		scanf("%d",&k);
		for(int i=1;i<=k;i++) scanf("%d",&n[i]);
		for(int i=1;i<k;i++) scanf("%d",&m[i]);
		for(int i=1;i<k;i++){
			tmp.n=n[i],tmp.m=n[i+1];
			for(int x=1;x<=n[i];x++)
				for(int y=1;y<=n[i+1];y++) tmp.x[x][y]=0;
			for(int j=1;j<=m[i];j++) scanf("%d%d",&x,&y),tmp.x[x][y]=1;
			ans=i==1?tmp:ans*tmp;
		}
		printf("%d\n",det(ans.x,n[1]));
	}
	return 0;
}

七、特征多项式

注：还不太懂的时候写的，目的仅为学习下面那道例题。

1. 定义

对于一个 \(n\times n\) 的矩阵 \(A\)，若存在一个数 \(\lambda\)、一个 \(n\times 1\) 的非零向量 \(\vec v\)，满足 \(\lambda\vec v=A\vec v\)，则称 \(\lambda\) 与 \(\vec v\) 为一组对应的 特征值 与 特征向量。

\[\begin{bmatrix} \lambda v_1\\ \lambda v_2\\ \vdots\\ \lambda v_n \end{bmatrix}=\begin{bmatrix} A_{1,1}&A_{1,2}&\cdots&A_{1,n}\\ A_{2,1}&A_{2,2}&\cdots&A_{2,n}\\ \vdots&\vdots&\ddots&\vdots\\ A_{n,1}&A_{n,2}&\cdots&A_{n,n} \end{bmatrix} \times \begin{bmatrix} v_1\\ v_2\\ \vdots\\ v_n \end{bmatrix} \]

\(\lambda I\times \vec v=A\times \vec v\Leftrightarrow(\lambda I-A)\vec v=0\)，其中 \(I\) 表示单位矩阵。

固定 \(\lambda\)，存在对应的非零向量 \(\vec v\) \(\Leftrightarrow\) \(\text{det}(\lambda I-A)=0\)。也就是说，\(\lambda\) 是一个特征值 \(\Leftrightarrow\) \(\text{det}(\lambda I-A)=0\)，且一个特征值 \(\lambda\) 能对应无穷多个 \(\vec v\)。

我们将 \(\text{det}(\lambda I-A)=0\) 称为矩阵 \(A\) 的 特征方程，特征值 \(\lambda\) 就相当于是特征方程的解。

根据行列式的计算公式，\(\text{det}(\lambda I-A)\) 是一个关于 \(\lambda\) 的 \(n\) 次多项式 \(p_A(\lambda)\)，称 \(p_A(\lambda)\) 为矩阵 \(A\) 的 特征多项式。\(p_A(\lambda)=0\) 的解就是特征值 \(\lambda\)。

\[p_A(\lambda)=\text{det}(\lambda I-A)=\det( \begin{bmatrix} \lambda -A_{1,1}&-A_{1,2}&\cdots&-A_{1,n}\\ -A_{2,1}&\lambda-A_{2,2}&\cdots&-A_{2,n}\\ \vdots&\vdots&\ddots&\vdots\\ -A_{n,1}&-A_{n,2}&\cdots&\lambda-A_{n,n} \end{bmatrix} ) \]

设 \(p_A(\lambda)\) 的 \(n\) 个零点为 \(\lambda_1,\lambda_2,\cdots,\lambda_n\)（可能重复。这恰好是 \(A\) 的 \(n\) 个特征值），则 \(p_A(\lambda)=(\lambda-\lambda_1)(\lambda-\lambda_2)\cdots(\lambda-\lambda_n)\)。那么：

• \(p_A(0)=(-1)^n\prod_{i=1}^n\lambda_i\)

  而 \(p_A(0)=\text{det}(0I-A)=(-1)^n\text{det}(A)\)，

  所以 \(\text{det}(A)=\prod_{i=1}^n\lambda_i\)。

• \([\lambda^1]p_A(\lambda)=(-1)^{n-1}\sum_{i=1}^n(\prod_{j\neq i}\lambda_j)\)

  如果 \(p_A(0)=0\)（意味着 \(\text{det}(A)=0\)），则存在某个 \(\lambda_j=0\)，这里的 \(\sum\) 就可以去掉变成只有一项。

• \([\lambda^n]p_A(\lambda)=1\)

• \([\lambda^{n-1}]p_A(\lambda)=\sum_{i=1}^n-\lambda _i\)

  并且 \([\lambda^{n-1}]\text{det}(\lambda I-A)=[\lambda^{n-1}](\lambda-A_{1,1})(\lambda-A_{2,2})\cdots(\lambda-A_{n,n})\)\(=\sum_{i=1}^n-A_{i,i}\)（因为 \(p\) 不选主对角线至多只能 \(\lambda^{n-2}\)），

  所以 \(tr(A)=\sum_{i=1}^n\lambda_i\)。其中 \(tr(A)=\sum_{i=1}^n A_{i,i}\) 被称为矩阵的迹。

2. 性质

\[\begin{aligned} \text{det}(\lambda I-A)&=\sum_p(-1)^{\tau(p)}\prod_{i=1}^n(\lambda I-A)_{i,p_i}\\ &=\sum_p(-1)^{\tau(p)}\prod_{i=1}^n \begin{cases} -A_{i,p_i}&(i\neq p_i)\\ \lambda-A_{i,p_i}&(i=p_i) \end{cases}\\ &=\sum_p(-1)^{\tau(p)}\sum_{S\subseteq\{i\mid i=p_i\}}\prod_{i=1}^n \begin{cases} \lambda&(i\in S)\\ -A_{i,p_i}&(i\notin S) \end{cases}\\ &=\sum_{S\subseteq\{1,2,\cdots,n\}}\lambda^{|S|}(-1)^{n-|S|}\text{det}(A\ 删去\ S\ 中的行列) \end{aligned} \]

\(S\) 表示，\(i=p_i\) 的 \(i\) 里，乘起来时选了 \(\lambda\) 的 \(i\) 的集合。

则：

• \([\lambda^0]p_A(\lambda)=(-1)^n\text{det(A)}\)
• \([\lambda^1]p_A(\lambda)=(-1)^{n-1}\sum_{i=1}^n\text{det}(A\ 删去第\ i\ 行第\ i\ 列)\)
• \([\lambda^k]p_A(\lambda)=(-1)^{n-k}\times\) \(A\) 的所有 \(n-k\) 阶主子式的行列式之和

3. 循环矩阵的特征值

对于循环矩阵 \(A\)：

\[A= \begin{bmatrix} a_0&a_1&\cdots&a_{n-1}\\ a_{n-1}&a_0&\cdots&a_{n-2}\\ \vdots&\vdots&\ddots&\vdots\\ a_1&a_2&\cdots&a_0 \end{bmatrix} \]

取 \(\omega\) 为某个 \(n\) 次单位根，利用 \(\omega\) 的幂次循环移位的性质，令：

\[\vec v= \begin{bmatrix} \omega^0\\ \omega^1\\ \omega^2\\ \cdots\\ \omega^{n-1} \end{bmatrix}\\ \lambda=a_0\omega^0+a_1\omega^1+a_2\omega^2+\cdots+a_{n-1}\omega^{n-1} \]

则 \(A\vec v=\lambda\vec v\)。由于 \(\omega\) 有 \(n\) 种取值，故 \(A\) 的 \(n\) 个特征值都找到了。

或者说，设第一行的生成函数为 \(f(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1}\)，则 \(\lambda=f(\omega)\)。

又因为 \(\text{det}(A)=\prod_{i=1}^n\lambda_i\)，可以 FFT 算循环矩阵的行列式。

4. 例题

2023.8.13 特征多项式 + 循环矩阵 + 单位根反演

有 \(A\times B\) 个节点，编号为 \(0\sim AB-1\)，\((x,y)\in E\) 当且仅当 \(x\not\equiv y\pmod A\) 且 \(x\not\equiv y\pmod B\)。

求这张图的生成树数量 \(\bmod 998244353\)。

\(1\leq A,B\leq 10^{18}\)。

利用特征多项式进行转化：

• 记 \(A\) 表示度数矩阵 - 邻接矩阵，\(ans=\text{det}(A\ 删去一行一列)\)。

  而 \([\lambda^1]p_A(\lambda)=(-1)^{n-1}\sum_{i=1}^n\text{det}(A\ 删去第\ i\ 行第\ i\ 列)\)，且根据矩阵树定理，删去哪一行哪一列都平等，故 \(ans=\frac{1}{n}\sum_{i=1}^n\text{det}(A\ 删去第\ i\ 行第\ i\ 列)=\frac{1}{n}\frac{[\lambda^1]p_A(\lambda)}{(-1)^{n-1}}\)。转化为了求 \([\lambda^1]p_A(\lambda)\)。

  又因为 \([\lambda^1]p_A(\lambda)=(-1)^{n-1}\sum_{i=1}^n(\prod_{j\neq i}\lambda_j)\)，我们实际上要算 \(\sum_{i=1}^n(\prod_{j\neq i}\lambda_j)\)。

打表发现 \(A\) 是个循环矩阵。循环矩阵的特征多项式有特殊性质：

• 可以想想 \(0\) 的出边与 \(1\) 的出边有什么区别。相当于邻接矩阵 \(0\) 那行向右平移了一格，就得到了邻接矩阵 \(1\) 那行。并且 \(0,1\) 的度数是一样的，度数矩阵也是整行向右平移一格。

  设 \(A\) 中 \(0\) 那行的生成函数为 \(f(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1}\)，则 \(\lambda=f(\omega)\)。

然后找 \(A\) 自己的性质：

• 显然 \(A\) 中一行的和为 \(0\)，所以 \(\lambda_1=f(\omega^0)=0\)，\(\sum_{i=1}^n(\prod_{j\neq i}\lambda_j)=\prod_{j\neq 1}\lambda_j=\prod_{k=1}^{n-1} f(\omega^k)\)。

• 容斥：

  \[f(x)=AB-A-B+\gcd(A,B)-(\sum_{i=0}^{AB-1}x^i-\sum_{i=0}^{B-1}x^{iA}-\sum_{i=0}^{A-1}x^{iB}+\sum_{i=0}^{\gcd(A,B)-1}x^{i\,\text{lcm}(A,B)}) \]

  其中，前半部分是度数矩阵第 \(0\) 行的生成函数，后半部分是邻接矩阵 \(0\) 那行的生成函数。

• 根据单位根反演，\([n\mid x]=\frac{1}{n}\sum_{i=0}^{n-1}(\omega_n^k)^i\)，所以 \(\sum_{i=0}^{AB-1}(\omega_{AB}^k)^i=[AB\mid k]AB\)，\(\sum_{i=0}^{B-1}(\omega_{AB}^k)^{iA}=\sum_{i=0}^{B-1}(\omega_{AB}^{kA})^i=\sum_{i=0}^{B-1}(\omega_B^k)^i=[B\mid k]B\)。同理：

  \[\begin{aligned} f(\omega^k) &=AB-A-B+\gcd(A,B)-([AB\mid k]AB-[B\mid k]B-[A\mid k]A+[\gcd(A,B)\mid k]\gcd(A,B))\\ &=[AB\not\mid k]AB-[A\not\mid k]A-[B\not\mid k]B+[\gcd(A,B)\not\mid k]\gcd(A,B) \end{aligned} \]

怎么求 \(\sum_{k=1}^{n-1}f(\omega^k)\)？

• 首先由于 \(1\leq k<AB\)，不可能 \(AB\mid k\)。

  然后分类讨论考虑 \(A\mid k,B\mid k,\gcd(A,B)\mid k\)。

  1. \(A\mid k\land B\mid k\Rightarrow \gcd(A,B)\mid k\)：每个 \(k\) 的贡献为 \(AB\)，这样的 \(k\) 的个数为 \(\frac{AB}{\text{lcm}(A,B)}-1\)（去掉 \(AB\)）。

  2. \(A\mid k\land B\not\mid k\Rightarrow \gcd(A,B)\mid k\)：\(AB-B\)，\(\frac{AB}{A}-\frac{AB}{\text{lcm}(A,B)}\)。

  3. \(A\not\mid k\land B\mid k\Rightarrow \gcd(A,B)\mid k\)：\(AB-A\)，\(\frac{AB}{B}-\frac{AB}{\text{lcm}(A,B)}\)。

  4. \(A\not\mid k\land B\not\mid k\)：

     若 \(\gcd(A,B)\mid k\)，\(AB-A-B\)，\(\frac{AB}{\gcd(A,B)}-\frac{AB}{A}-\frac{AB}{B}+\frac{AB}{\text{lcm}(A,B)}\)；

     若 \(\gcd(A,B)\not\mid k\)，\(AB-A-B+\gcd(A,B)\)，\(AB-\frac{AB}{\gcd(A,B)}\)。

总结一下 \(ans\) 的转化：

\[ans=\frac{1}{n}\frac{[\lambda^1]p_A(\lambda)}{(-1)^{n-1}}=\frac{1}{n}\frac{(-1)^{n-1}\sum_{i=1}^n(\prod_{j\neq i}\lambda_j)}{(-1)^{n-1}}=\frac{1}{n}\prod_{j\neq 1}\lambda_j=\frac{1}{n}\prod_{k=1}^{n-1} f(\omega^k) \]

#include<bits/stdc++.h>
#define ll __int128
using namespace std;
const int mod=998244353;
long long a,b,d;
ll n;
int qpow_(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod)
		if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
int qpow(ll x,ll n){return qpow_(x%mod,n%(mod-1));}
signed main(){
	scanf("%lld%lld",&a,&b),n=(ll)a*b,d=__gcd(a,b);
	int ans=1ll*qpow(n,d-1)*qpow(n-b,b-d)%mod*qpow(n-a,a-d)%mod*qpow(n-a-b,n/d-a-b+d)%mod*qpow(n-a-b+d,n-n/d)%mod;
	printf("%lld\n",1ll*ans*qpow(n,mod-2)%mod);
	return 0;
}