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;
}