线性代数相关
0. 行列式
0.1. 符号与定义
\(\mathrm{det}(A)\)，又记作 \(|A|\)，等于 \(\sum_{p}(-1)^{\tau(p)}\prod_{i=1}^nA_{i,p_i}\)，其中 \(p\) 为 \(1\sim n\) 的排列，\(\tau(p)\) 表示排列 \(p\) 的逆序对数。
0.2. 基本性质
-
Lemma 0.2.1. 对于一个上三角矩阵，它的行列式等于主对角线所有值的乘积。
Proof. 根据行列式定义，要使 \(a_{n,p_n}\neq 0\)，\(p_n\) 只能取 \(n\)；要使 \(a_{n-1,p_{n-1}}\neq 0\)，\(p_{n-1}\) 只能取 \(n-1\) 或 \(n\) ，而 \(p_n\) 已经等于 \(n\) 所以 \(p_{n-1}\) 只能等于 \(n-1\) …… 以此类推，得到唯一的排列 \(p_i=i\)，此时 \(\tau(p)=0\)。\(\square\)
-
Lemma 0.2.2. 单位矩阵的行列式为 \(1\)。即 \(|I|=1\)。
根据 \(I_{i,j}=[i=j]\) 以及 Lemma 0.2.1. 易证。
-
Lemma 0.2.3. 交换矩阵两行，行列式变号。
Proof. 交换排列的两个数会让排列逆序对数量改变奇偶性，因此相当于为每个 \((-1)^{\tau(p)}\prod_{i=1}^na_{i,p_i}\) 乘上了 \(-1\)。根据乘法分配律将 \(-1\) 提出。\(\square\)
-
Lemma 0.2.4. 将矩阵某一行乘上常数，行列式乘上相同常数。
Proof. 不妨设 \(A_i\gets A_i\times k\)，记变换得到的矩阵为 \(A'\)。因为 \(A'_{i,j}=kA_{i,j}\) 且 \(A'_{c,j}=A_{c,j}\ (c\neq i)\) 则 \(\prod_{i=1}^nA'_{i,p_i}=k\prod_{i=1}^nA_{i,p_i}\)。根据乘法分配律将 \(k\) 提出。\(\square\)
-
Lemma 0.2.5. 若矩阵有相同的两行，则行列式为 \(0\)。
Proof. 交换相同的两行会让行列式变号，但矩阵保持不变。即 \(|A|=-|A|\)，得到 \(|A|=0\)。\(\square\)
-
Lemma 0.2.6. 若矩阵有两行存在倍数关系，则行列式为 \(0\)。
结合 Lamma 0.2.4. 与 0.2.5. 易证。
-
Lemma 0.2.7. 若两个矩阵至多有一行不等，则将这不等的一行相加得到的新矩阵的行列式等于原矩阵行列式之和。
Proof. 不妨设不等的行编号为 \(c\)，则 \(A_i=B_i=S_i,A_c+B_c=S_c\ (i\neq c)\)。根据乘法分配律得 \(\prod_{i=1}^nS_{i,p_i}=\left(\prod_{i=1}^nA_{i,p_i}\right)+\left(\prod_{i=1}^nB_{i,p_i}\right)\)，即 \(|S|=|A|+|B|\)。\(\square\)
-
Lemma 0.2.8. 将矩阵的某一行加上另一行的倍数，行列式不变。
Proof. 不妨设 \(A_d\gets A_d+kA_c\ (c\neq d)\)，令操作后的 \(A\) 矩阵为 \(C\)。令 \(B_i=A_i,B_d=kA_c\ (i\neq d)\)，则 \(B_c=B_d\)。根据 Lemma 0.2.6. 可知 \(|B|=0\)。根据 Lemma 0.2.7 可知 \(|A|+|B|=|C|\)，即 \(|C|=|A|\)。\(\square\)
1. 高斯消元
1.1. 算法介绍
高斯消元用来求解形如
的 \(m\) 元线性方程组。不妨记 \(a_{i,m}=b_i\)，则 \(a\) 被称为该线性方程组的增广矩阵。
高斯消元的过程如下：
- 记已经处理好的方程个数为 \(r\)。枚举每一个未知数 \(x_c\ (0\leq c<m)\) 作为主元。
- 找到 \(j\in[r,n)\) 使得 \(a_{j,c}\) 最大。交换 \(a_r\) 与 \(a_j\) 两行。
- 若 \(a_{r,c}=0\)，则 \(x_c\) 要么无解，要么有无数解。回到第一步枚举下一个主元，即 \(c\gets c+1\)。
- 将 \(a_r\) 中的每个数除以 \(a_{r,c}\)。注意需要预先存储 \(a_{r,c}\) 或倒序枚举第二维。该操作后 \(a_{r,c}=1\)。
- 将 \(a_{j,c}\ (j\in(r,n))\) 消为 \(0\)。具体来说，将 \(a_j\) 减去 \(\dfrac{a_{j,c}}{a_{r,c}}a_i\)。因为 \(a_{r,c}=1\)，故枚举 \(k\in[c,m)\) 并将 \(a_{j,k}\) 减去 \(a_{j,c}\times a_{i,k}\)。
最终我们得到一个上三角矩阵，其中第 \(r\) 行至第 \(n-1\) 行的前 \(m\) 个数都为 \(0\)。
如果最终 \(r<n\) 则表明方程组无解或有无数解：若存在 \(i\in [r,n)\) 满足 \(a_{i,n}\neq 0\)，即出现 \(0\) 等于一个 \(\neq 0\) 的数，则方程组无解；否则 \(n-r\) 即为自由元个数。注意如果 \(n>m\) 且方程组有解，则需要 \(n\gets m\)，再根据 \(r\) 与 \(n\) 的大小判断是否有无数解。
此时显然有 \(n=m\)。
接下来只需要将每个解回带即可。
- 显然 \(a_{n-1,m}\) 已经是 \(x_{n-1}\) 的解。
- 枚举从大到小枚举行数 \(i\in [0,n-1)\)，只需要将所有 \(x_j\ (j\in(i,n))\) 带入即可。注意此时 \(a_{j,n}\) 已经是 \(x_j\) 的解，所以只需要用 \(a_{i,n}\) 减去 \(\sum_{j=i+1}^{n-1}a_{i,j}\times a_{j,n}\) 即为 \(x_i\) 的解。
时间复杂度 \(n^3\)。
下面给出 \(n=m\) 时的代码。\(n\neq m\) 时只需稍作修改。
int gauss(){ // 返回方程自由元的个数
int r=0,c=0;
for(;c<n;c++){
int cur=r;
for(int i=r;i<n;i++)if(fabs(a[i][c])>fabs(a[cur][c]))cur=i;
if(fabs(a[cur][c])<eps)continue;
swap(a[cur],a[r]);
for(int i=n;i>=c;i--)a[r][i]/=a[r][c];
for(int i=r+1;i<n;i++)if(fabs(a[i][c])>eps)
for(int j=n;j>=c;j--)a[i][j]-=a[r][j]*a[i][c];
r++;
}
if(r<n){
for(int i=r;i<n;i++)if(fabs(a[i][n])>eps)return -1;
return n-r;
}
for(int i=n-2;~i;i--)for(int j=i+1;j<n;j++)
a[i][n]-=a[j][n]*a[i][j];
return 0;
}
1.2. 应用：矩阵求逆
给出 \(A\)，求它在模数 \(p\) 下的逆。设单位矩阵为 \(E\)，\(A^{-1}=P\)，则显然有 \(PA=E\)。
将 \(A\) 进行线性变换相当于左乘一个矩阵，而将 \(A\) 线性变换到 \(E\) 的过程就等于左乘它的逆，即 \(P\)。因为相同的线性变换对应相同的左乘矩阵，因此令 \(B=E\)，则 \(PB=P\)，所以在高斯消元 \(A\) 时，对 \(B\) 进行相同的变换即可得到 \(P\)。
注意点：此时枚举列下标时要从 \(1\) 开始而不是 \(i\)，因为虽然 \(A_{i,j}=0\ (i>j)\)，但其对应的 \(B_{i,j}\) 并不一定为 \(0\)。
时间复杂度 \(n^3\)。
1.3. 应用：求行列式
根据定义直接求的时间复杂度为 \(n!n\)，显然无法承受。但如果我们能将 \(A\) 转化为上三角矩阵，那么则当且仅当排列 \(p\) 满足 \(p_i=i\) 时，\(\prod_{i=1}^na_{i,p_i}\neq 0\)。
根据 Lemma 0.2.3（\(\mathrm{swap}(A_i,A_j)\ (i\neq j)\)，则 \(|A|\) 反号），0.2.4（\(A_i\gets k A_i\)，则 \(|A|\gets k|A|\)）与 0.2.8（\(A_i\gets A_i + k A_j\ (i\neq j)\)，则 \(|A|\) 不变）可以高斯消元 \(n^3\) 求行列式。代码见下方例题 IV.
1.4. 例题
I. P3389 【模板】高斯消元法 - 洛谷
板子题。
#include <bits/stdc++.h>
using namespace std;
#define gc getchar()
#define pb push_back
inline int read(){
int x=0; char s=gc;
while(!isdigit(s))s=gc;
while(isdigit(s))x=(x<<1)+(x<<3)+s-'0',s=gc;
return x;
}
const int N=105;
const double eps=1e-8;
int n;
double a[N][N];
int gauss(){
int r=0,c=0;
for(;c<n;c++){
int cur=r;
for(int i=r;i<n;i++)if(fabs(a[i][c])<fabs(a[cur][c]))cur=i;
if(fabs(a[cur][c])<eps)continue;
swap(a[cur],a[r]);
for(int i=n;i>=c;i--)a[r][i]/=a[r][r];
for(int i=r+1;i<n;i++)if(fabs(a[i][c])>eps)
for(int j=n;j>=c;j--)a[i][j]-=a[r][j]*a[i][c];
r++;
}
if(r<n){
for(int i=r;i<n;i++)if(fabs(a[i][n])>eps)return -1;
return n-r;
}
for(int i=n-2;~i;i--)for(int j=i+1;j<n;j++)
a[i][n]-=a[j][n]*a[i][j];
return 0;
}
int main(){
cin>>n;
for(int i=0;i<n;i++)for(int j=0;j<=n;j++)cin>>a[i][j];
if(gauss()!=0)puts("No Solution"),exit(0);
for(int i=0;i<n;i++)printf("%.2lf\n",a[i][n]);
return 0;
}
II. P2455 [SDOI2006]线性方程组
仍然是板子题。
#include <bits/stdc++.h>
using namespace std;
#define gc getchar()
#define pb push_back
inline int read(){
int x=0; char s=gc;
while(!isdigit(s))s=gc;
while(isdigit(s))x=(x<<1)+(x<<3)+s-'0',s=gc;
return x;
}
const int N=105;
const double eps=1e-8;
int n;
double a[N][N];
int gauss(){
int r=0,c=0;
for(;c<n;c++){
int cur=r;
for(int i=r;i<n;i++)if(fabs(a[i][c])>fabs(a[cur][c]))cur=i;
if(fabs(a[cur][c])<eps)continue;
swap(a[cur],a[r]);
for(int i=n;i>=c;i--)a[r][i]/=a[r][c];
for(int i=r+1;i<n;i++)if(fabs(a[i][c])>eps)
for(int j=n;j>=c;j--)a[i][j]-=a[r][j]*a[i][c];
r++;
}
if(r<n){
for(int i=r;i<n;i++)if(fabs(a[i][n])>eps)return -1;
return 0;
}
for(int i=n-2;~i;i--)for(int j=i+1;j<n;j++)
a[i][n]-=a[j][n]*a[i][j];
return 1;
}
int main(){
cin>>n;
for(int i=0;i<n;i++)for(int j=0;j<=n;j++)cin>>a[i][j];
int c=gauss();
if(c!=1)cout<<c<<endl;
else for(int i=0;i<n;i++)printf("x%d=%.2lf\n",i+1,a[i][n]);
return 0;
}
III. P4783 【模板】矩阵求逆
板子题。具体算法见 6.2.
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N=4e2+5;
const int mod=1e9+7;
ll ksm(ll a,ll b){
ll s=1;
while(b){
if(b&1)s=s*a%mod;
a=a*a%mod,b>>=1;
}
return s;
}
ll inv(ll x){return ksm(x,mod-2);}
ll n,a[N][N],b[N][N];
void gauss(){
for(int i=1;i<=n;i++){
int p=-1;
for(int j=i;j<=n;j++)if(a[j][i]){p=j; break;}
if(p==-1)puts("No Solution"),exit(0);
swap(a[i],a[p]),swap(b[i],b[p]);
ll iv=inv(a[i][i]);
for(int j=1;j<=n;j++)a[i][j]=a[i][j]*iv%mod,b[i][j]=b[i][j]*iv%mod;
for(int j=i+1;j<=n;j++){
ll c=a[j][i];
for(int k=n;k;k--)
b[j][k]=(b[j][k]-c*b[i][k]%mod+mod)%mod,
a[j][k]=(a[j][k]-c*a[i][k]%mod+mod)%mod;
}
}
for(int i=n-1;i;i--)
for(int j=i+1;j<=n;j++)
for(int k=1;k<=n;k++)
b[i][k]=(b[i][k]-a[i][j]*b[j][k]%mod+mod)%mod;
}
int main(){
cin>>n;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
cin>>a[i][j],b[i][j]=i==j;
gauss();
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
cout<<b[i][j]<<(j==n?'\n':' ');
return 0;
}
IV. P7112 【模板】行列式求值
注意到给出的模数 \(p\) 不一定为质数，所以无法直接求逆元。
如果直接让两行辗转相除，就不需要用到精确的除法也可以进行消元。时间复杂度看似 \(n^3\log p\) 实际上可以证明是 \(n^2(n+\log p)\) 的。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N=600+5;
ll n,p,a[N][N],sym,prod=1;
int main(){
cin>>n>>p;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
cin>>a[i][j],a[i][j]%=p;
for(int i=1;i<n;i++){
for(int j=i+1;j<=n;j++){
while(a[i][i]){
ll q=a[j][i]/a[i][i];
for(int k=i;k<=n;k++)a[j][k]=(a[j][k]-q*a[i][k]%p+p)%p;
swap(a[i],a[j]),sym^=1;
} swap(a[i],a[j]),sym^=1;
}
}
for(int i=1;i<=n;i++)prod=(prod*a[i][i])%p;
cout<<(sym?p-prod:prod)%p<<endl;
return 0;
}
2. Matrix-Tree 定理
2.1. 算法介绍
前置知识：高斯消元。
记 \(D\) 为无向图的度数矩阵，即 \(D_{i,i}=deg_i,D_{i,j}=0\ (i\neq j)\)，记 \(E\) 为无向图的邻接矩阵。记 \(K=D-E\)，\(K’\) 为 \(K\) 的 \(n-1\) 阶主子式，则 \(|K'|\) 为该无向图的生成树个数。
可以有重边和自环。可以发现自环对 \(D\) 和 \(E\) 的影响可以相互抵消，因此自环不会影响 \(K\) 每个位置的值。
证明略。
2.2. 扩展：有向图生成树个数
既然是有向图那么生成树就有根。不妨设根为 \(r\)。
记 \(D_{in}\) 为无向图的入度矩阵，\(D_{out}\) 为无向图的出度矩阵。
\(K_{in}=D_{in}-E\)，\(K_{in}'\) 为 \(K_{in}\) 去掉第 \(r\) 行和第 \(r\) 列得到的矩阵，则外向生成树个数为 \(|K'_{in}|\)。同理，内向生成树个数为 \(|K'_{out}|\)。
2.3. 扩展：边权积的和
对于一条 \((u,v,w)\) 的边，相当于有 \(w\) 条 \((u,v)\) 这样的重边。因此可以归约到普通生成树计数问题。
2.4. 扩展：边权和的和
考虑将加法转化为乘法。不难想到将边权看做形如 \(1+wx\) 的多项式，则答案为该一次多项式矩阵行列式的一次项系数。因此我们只需维护形如 \(a+bx\) 的一次多项式即可。
乘法：\((a+bx)(c+dx)\equiv ac+(ad+bc)x\pmod {x^2}\)。
除法：\(\dfrac{a+bx}{c+dx}\equiv \dfrac{(a+bx)(c-dx)}{c^2-d^2x^2}\equiv \dfrac{a}{c}+\dfrac{bc-ad}{c^2}x\pmod{x^2}\)。
2.5. 例题
I. P6178 【模板】Matrix-Tree 定理
板子题。时间复杂度 \(n^3\)。
#include <bits/stdc++.h>
using namespace std;
#define gc getchar()
#define ll long long
inline int read(){
int x=0,sign=0; char s=gc;
while(!isdigit(s))sign|=s=='-',s=gc;
while(isdigit(s))x=(x<<1)+(x<<3)+(s-'0'),s=gc;
return sign?-x:x;
}
const int N=300+5;
const int mod=1e9+7;
ll ksm(ll a,ll b){
ll s=1;
while(b){
if(b&1)s=s*a%mod;
a=a*a%mod,b>>=1;
} return s;
}
ll n,m,t,a[N][N];
void add(ll &x,ll y){x=(x+y+mod)%mod;}
ll det(){
ll sgn=0,prod=1;
for(int i=2;i<n;i++){
for(int j=i+1;j<=n;j++)if(a[j][i]&&!a[i][i])swap(a[i],a[j]),sgn^1;
ll inv=ksm(a[i][i],mod-2);
for(int j=i+1;j<=n;j++){
ll tmp=inv*a[j][i]%mod;
for(int k=i;k<=n;k++)add(a[j][k],-tmp*a[i][k]%mod);
}
}
for(int i=2;i<=n;i++)prod=(prod*a[i][i])%mod;
return (sgn?mod-prod:prod)%mod;
}
int main(){
cin>>n>>m>>t;
for(int i=1;i<=m;i++){
int u=read(),v=read(),w=read();
add(a[u][v],-w),add(a[v][v],w);
if(!t)add(a[v][u],-w),add(a[u][u],w);
}
cout<<det()<<endl;
return 0;
}
*II. P3317 [SDOI2014]重建
不那么套路的 Matrix-Tree。
首先对于一个生成树它的权值应该是 \(\prod_{e\notin T}(1-p_e)\prod_{e\in T}p_e\) 而不仅仅只有后面一部分，看起来不太好搞的样子。
一步比较巧妙的转化：\(\prod_{e\notin T}(1-p_e)=\dfrac{\prod_{e}(1-p_e)}{\prod_{e\in T}(1-p_e)}\)，这样答案即为 \(\prod_e(1-p_e)\sum_T\prod_{e\in T}\dfrac{p_e}{1-p_e}\)，矩阵数定理解决。
如果出现 \(p_e=1\) 可以将其变为 \(1-eps\) 以避免乘或除以 \(0\) 的情况。时间复杂度 \(n^3\)。
#include <bits/stdc++.h>
using namespace std;
const int N=50+5;
const double eps=1e-8;
int n;
double ans=1,a[N][N];
double det(){
for(int i=2;i<n;i++){
int pos=i;
for(int j=i+1;j<=n;j++)if(fabs(a[j][i])>eps)pos=j;
swap(a[i],a[pos]);
if(fabs(a[i][i])<eps)return 0;
for(int j=i+1;j<=n;j++)
for(int k=n;k>=i;k--)
a[j][k]-=a[j][i]/a[i][i]*a[i][k];
}
for(int i=2;i<=n;i++)ans*=a[i][i];
return fabs(ans);
}
int main(){
cin>>n;
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
cin>>a[i][j];
if(i!=j){
if(1-a[i][j]<eps)a[i][j]-=eps;
if(i<j)ans*=(1-a[i][j]);
a[i][j]=-a[i][j]/(1-a[i][j]);
}
}
}
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
if(i!=j)a[i][i]-=a[i][j];
printf("%.5lf\n",det());
return 0;
}
III. P4111 [HEOI2015]小 Z 的房间
对相邻的两个 . 连一条边，则答案为最终连出来的无向图的生成树个数。时间复杂度 \((nm)^3\)。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N=85;
const int mod=1e9;
ll n,m,cnt,a[N][N],buc[N][N];
char s[N][N];
void add(int u,int v){a[u][u]++,a[v][v]++,a[u][v]--,a[v][u]--;}
int gauss(){
ll sgn=0,ans=1;
for(int i=2;i<=cnt;i++)
for(int j=i+1;j<=cnt;j++){
while(a[i][i]){
ll q=a[j][i]/a[i][i];
for(int k=i;k<=n*m;k++)a[j][k]=(a[j][k]-q*a[i][k]%mod+mod)%mod;
swap(a[i],a[j]),sgn^=1;
}
swap(a[i],a[j]),sgn^=1;
}
for(int i=2;i<=cnt;i++)ans=ans*a[i][i]%mod;
return (sgn?mod-ans:ans)%mod;
}
int main(){
cin>>n>>m;
for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)cin>>s[i][j],buc[i][j]=cnt+=s[i][j]=='.';
for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)if(s[i][j]=='.'){
if(s[i+1][j]=='.')add(buc[i][j],buc[i+1][j]);
if(s[i][j+1]=='.')add(buc[i][j],buc[i][j+1]);
}
cout<<gauss()<<endl;
return 0;
}
*IV. P6624 [省选联考 2020 A 卷] 作业题
从大到小枚举 \(\gcd d\) ，如果边权能整除 \(d\) 则保留该边，否则丢弃。然后做矩阵树定理求边权和的和。求得边权 \(\gcd\) 不小于 \(d\) 的生成树个数后再容斥，即减去 \(d\) 的倍数的答案，该部分是 \(w\ln w\) 的。
总时间复杂度 \(n^3w\)，显然过不去。
注意到连接某个点的所有边的不同约数个数至多为 \(n\max d(w)\)，而 \(\max d(w)=144\)，因此只需要对这最多 \(144n\) 个约数分别做矩阵树定理即可。时间复杂度 \(n^4 \max(d_w)\)。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define pll pair <ll,ll>
#define fi first
#define se second
const int N=30+5;
const int mod=998244353;
const int W=152502;
ll ksm(ll a,ll b){
ll s=1;
while(b){
if(b&1)s=s*a%mod;
b>>=1,a=a*a%mod;
}
return s;
}
pll operator + (pll a,pll b){return {(a.fi+b.fi)%mod,(a.se+b.se)%mod};}
pll operator - (pll a,pll b){return {(a.fi-b.fi%mod+mod)%mod,(a.se-b.se%mod+mod)%mod};}
pll operator * (pll a,pll b){return {a.fi*b.fi%mod,(a.fi*b.se+a.se*b.fi)%mod};}
pll operator / (pll a,pll b){
ll inv=ksm(b.fi,mod-2);
return {a.fi*inv%mod,(a.se*b.fi-a.fi*b.se%mod+mod)%mod*inv%mod*inv%mod};
}
ll n,m,cnt,d[W],sum[W],ans[W],e[N][N];
pll a[N][N];
ll solve(int x){
bool sgn=0; pll ans={1,0};
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++)
if(e[i][j]&&e[i][j]%x==0)
a[i][i]=a[i][i]-(a[i][j]={mod-1,mod-e[i][j]});
for(int i=2;i<n;i++){
for(int j=i+1;j<=n;j++)
if(a[i][j].fi){
swap(a[i],a[j]),sgn^=1;
break;
}
for(int j=i+1;j<=n;j++){
pll inv=a[j][i]/a[i][i];
for(int k=i;k<=n;k++)a[j][k]=a[j][k]-inv*a[i][k];
}
}
for(int i=2;i<=n;i++)ans=ans*a[i][i],assert(a[i][i].fi>=0&&a[i][i].se>=0);
return (sgn?mod-ans.se:ans.se)%mod;
}
int main(){
for(int i=1;i<W;i++)for(int j=i;j<W;j+=i)d[j]++;
cin>>n>>m;
for(int i=1;i<=m;i++){
int u,v,w; cin>>u>>v>>w;
e[u][v]=e[v][u]=w,sum[u]+=d[w],sum[v]+=d[w];
}
int id=1;
for(int i=2;i<=n;i++)if(sum[i]<sum[id])id=i;
for(int i=W-1;i;i--){
bool ok=0;
for(int j=1;j<=n;j++)ok|=e[id][j]&&e[id][j]%i==0;
if(ok){
ans[i]=solve(i);
for(int j=i+i;j<W;j+=i)ans[i]=(ans[i]-ans[j]+mod)%mod;
cnt=(cnt+i*ans[i])%mod;
}
}
cout<<cnt<<endl;
return 0;
}
V. P4336 [SHOI2016]黑暗前的幻想乡
感动，终于有一道自己想出来的 Matrix-Tree 定理题目了。
\(2^{n-1}\) 枚举所有承包公司是否被选上的情况，并对被选上的承包公司包含的所有边所形成的图进行矩阵树定理，得到仅使用被选上的承包公司承包的边所形成的生成树个数。显然是个容斥，即 \(\sum_{i=1}^{2^{n-1}-1}(-1)^{n-1-\mathrm{bitcount}(i)}ans_i\)。
时间复杂度 \(2^{n-1}n^3\)，可以进行适当剪枝，如若有节点度数为 \(0\) 则直接返回 \(0\) 等。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define mem(x,v) memset(x,v,sizeof(x))
#define pii pair <int,int>
#define fi first
#define se second
const int N=17;
const int mod=1e9+7;
ll ksm(ll a,ll b){
ll s=1;
while(b){
if(b&1)s=s*a%mod;
b>>=1,a=a*a%mod;
} return s;
}
int n;
vector <pii> e[N];
ll a[N][N],ans;
ll solve(int x){
mem(a,0);
for(int i=0;i<n-1;i++)
if((x>>i)&1)
for(pii it:e[i])
a[it.fi][it.se]=--a[it.se][it.fi];
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++)
if(i!=j)
a[i][i]-=a[i][j];
if(a[i][i]==0)return 0;
}
ll sgn=0,ans=1;
for(int i=2;i<=n;i++){
for(int j=i+1;j<=n;j++)
if(a[j][i]){
swap(a[i],a[j]),sgn^=1;
break;
}
ll inv=ksm(a[i][i],mod-2);
for(int j=i+1;j<=n;j++)
for(int k=n;k>=i;k--)
a[j][k]=(a[j][k]-a[j][i]*inv%mod*a[i][k]%mod+mod)%mod;
}
for(int i=2;i<=n;i++)ans=ans*(mod+a[i][i])%mod;
return (sgn?mod-ans:ans)%mod;
}
int main(){
cin>>n;
for(int i=0;i<n-1;i++){
int m,u,v; cin>>m;
while(m--)cin>>u>>v,e[i].push_back({u,v});
}
for(int i=1;i<(1<<n-1);i++){
int cnt=0;
for(int j=0;j<n-1;j++)cnt+=(i>>j)&1;
if((n-1-cnt)&1)ans=(ans-solve(i)+mod)%mod;
else ans=(ans+solve(i))%mod;
}
cout<<ans<<endl;
return 0;
}
VI. P4455 [CQOI2018]社交网络
板子题。边的节点颠倒一下就是求以 \(1\) 为根的外向生成树个数。
#include <bits/stdc++.h>
using namespace std;
const int N=255;
const int mod=1e4+7;
int ksm(int a,int b){
int s=1;
while(b){
if(b&1)s=s*a%mod;
b>>=1,a=a*a%mod;
} return s;
}
int n,m,sgn,ans=1,a[N][N];
int main(){
cin>>n>>m;
for(int i=1,u,v;i<=m;i++)cin>>u>>v,a[v][u]=-1;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
a[j][j]-=(i!=j)*a[i][j];
for(int i=2;i<=n;i++){
for(int j=i+1;j<=n;j++)
if(a[j][i]){
swap(a[i],a[j]),sgn^=1;
break;
}
int inv=ksm(a[i][i],mod-2);
for(int j=i+1;j<=n;j++)
for(int k=n;k>=i;k--)
a[j][k]=(a[j][k]-a[j][i]*inv%mod*a[i][k]%mod+mod)%mod;
}
for(int i=2;i<=n;i++)ans=ans*(a[i][i]+mod)%mod;
cout<<(sgn?mod-ans:ans)%mod<<endl;
return 0;
}
*VII. CF917D Stranger Trees
nb tea！
类似求边权和的和，如果我们将给出树中的边权看做 \(x^1\)，其余边看做 \(x^0\)，那么恰好取到 \(k\) 条边的生陈树总数即为 \(x^k\) 前的系数。
显然对这个多项式矩阵做矩阵树定理是不现实的。但是假设最终得到的多项式为 \(a_0x^0+a_1x^1+\cdots+a_{n-1}x^{n-1}\)，我们可以插值求出每一项的系数，即将 \(x\) 带入 \(n\) 个取值，用矩阵树定理得到右边的常数项，高斯消元解 \(n\) 个 \(n\) 元一次方程即可。
一共要插 \(n\) 个值，每次插值 \(n^3\)。故时间复杂度 \(n^4\)。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N=100+5;
const ll mod=1e9+7;
ll ksm(ll a,ll b){
ll s=1; while(b){
if(b&1)s=s*a%mod;
b>>=1,a=a*a%mod;
} return s;
}
ll n,e[N][N],a[N][N],c[N][N];
void solve(ll x){
memset(a,0,sizeof(a));
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
if(i!=j){
if(e[i][j])a[i][j]=-x;
else a[i][j]=-1;
a[i][i]-=a[i][j];
}
ll sgn=0,ans=1;
for(int i=2;i<=n;i++){
for(int j=i+1;j<=n;j++)
if(a[j][i]){
swap(a[i],a[j]),sgn^=1;
break;
}
ll inv=ksm(a[i][i],mod-2);
for(int j=i+1;j<=n;j++)
for(int k=n;k>=i;k--)
a[j][k]=(a[j][k]-a[j][i]*inv%mod*a[i][k])%mod;
}
for(int i=2;i<=n;i++)ans=ans*(a[i][i]+mod)%mod;
c[x][n+1]=(sgn?mod-ans:ans)%mod;
for(ll i=1,j=1;i<=n;i++)c[x][i]=j,j=j*x%mod;
}
void gauss(){
memcpy(a,c,sizeof(c));
for(int i=1;i<n;i++){
for(int j=i+1;j<=n;j++)
if(a[j][i]){
swap(a[i],a[j]);
break;
}
ll inv=ksm(a[i][i],mod-2);
for(int j=i+1;j<=n;j++)
for(int k=n+1;k>=i;k--)
a[j][k]=(a[j][k]-a[j][i]*inv%mod*a[i][k])%mod;
}
for(int i=n-1;i;i--)
for(int j=i+1;j<=n;j++)
a[i][n+1]=(a[i][n+1]-a[i][j]*a[j][n+1]%mod*ksm(a[j][j],mod-2))%mod;
}
int main(){
cin>>n;
for(int i=1,u,v;i<n;i++)cin>>u>>v,e[u][v]=e[v][u]=1;
for(int i=1;i<=n;i++)solve(i);
gauss();
for(int i=1;i<=n;i++)cout<<(a[i][n+1]*ksm(a[i][i],mod-2)%mod+mod)%mod<<" ";
return 0;
}
