一些定义如下：
度数矩阵：一个
n
n
n 阶无向图
G
G
G 的度数矩阵
D
D
D 的大小为
n
×
n
n \times n
n×n，且
D
D
D 仅在
D
i
,
i
D_{i,i}
Di,i(
i
=
1
,
2
,
⋯
n
i=1,2,\cdots n
i=1,2,⋯n) 处有值，其中
D
i
,
i
D_{i,i}
Di,i 的值为点
i
i
i 的度数。(出 / 入) 度矩阵：一个
n
n
n 阶有向图
G
G
G 的 (出 / 入) 度矩阵
D
D
D 的大小为
n
×
n
n\times n
n×n，且
D
D
D 仅在
D
i
,
i
D_{i,i}
Di,i(
i
=
1
,
2
,
⋯
n
i=1,2,\cdots n
i=1,2,⋯n) 处有值，其中
D
i
,
i
D_{i,i}
Di,i 的值为点
i
i
i 的 (出 / 入) 度。邻接矩阵 —— 无向图：一个
n
n
n 阶无向图
G
G
G 的邻接矩阵
A
A
A 的大小为
n
×
n
n \times n
n×n，
A
A
A 在
A
i
,
j
A_{i,j}
Ai,j(
i
≠
j
i\neq j
i=j) 处有值，其中
A
i
,
j
A_{i,j}
Ai,j 的值为点
i
,
j
i,j
i,j 间的连边个数。邻接矩阵 —— 有向图：一个
n
n
n 阶有向图
G
G
G 的邻接矩阵
A
A
A 的大小为
n
×
n
n \times n
n×n，
A
A
A 在
A
i
,
j
A_{i,j}
Ai,j(
i
≠
j
i\neq j
i=j) 处有值，其中
A
i
,
j
A_{i,j}
Ai,j 的值为从点
i
i
i 连向点
j
j
j 的有向边个数。
矩阵树定理
由于此定理的证明极其复杂，需要大量高等数学知识，这里只给出结论。
矩阵树定理是用于求解图上生成树计数问题的重要定理，内容如下：
设
G
G
G 为一
n
n
n 阶无向图，定义
G
G
G 的基尔霍夫(Kirchhoff) 矩阵
K
K
K 为其度数矩阵与其邻接矩阵之差，则
G
G
G 的无根生成树的个数为
K
K
K 的任意一个
n
−
1
n - 1
n−1 阶主子式对应行列式的绝对值。
定理内容十分简洁，没什么好讲的。代码实现上基尔霍夫矩阵可以直接按照定义计算，对行列式求值用高斯消元即可，时间复杂度
O
(
n
3
)
\Omicron(n^3)
O(n3)。
下面给出高斯消元求行列式的代码：
double det(int n
){
double d
= 1;
for(int i
= 1; i
<= n
; i
++){
if(fabs(a
[i
][i
]) < eps
){
bool flag
= 0;
for(int j
= i
+ 1; j
<= n
; j
++){
if(fabs(a
[j
][i
]) > eps
){
d
*= -1,flag
= 1;
for(int k
= i
; k
<= n
; k
++) swap(a
[i
][k
],a
[j
][k
]);
break;
}
}
if(!flag
) return 0;
}
for(int j
= i
+ 1; j
<= n
; j
++){
double t
= a
[j
][i
] / a
[i
][i
];
for(int k
= i
; k
<= n
; k
++) a
[j
][k
] -= t
* a
[i
][k
];
}
d
*= a
[i
][i
];
}
return d
;
}
有向图上的推广
设
G
G
G 为一
n
n
n 阶有向图，定义
G
G
G 的基尔霍夫(Kirchhoff) 矩阵
K
K
K 为其 (出 / 入) 度矩阵与其邻接矩阵之差，则
G
G
G 中以任一点
i
i
i 为根的 (内 / 外) 向生成树的个数为
K
i
,
i
K_{i,i}
Ki,i 的余子式对应行列式的绝对值。
变元矩阵树定理
将上述两个定理中的度数、边数改为相应边边权和，则得到的值从生成树个数变为 对应的生成树所有边权之积 的和。
容易发现矩阵树定理是此定理每条边边权取一时的特殊情况。
例题
[SDOI2014]重建
分析
下面我们用
e
e
e 表示一条边，
p
e
p_e
pe 表示
e
e
e 连通的概率。 容易发现题目让我们求的是这个式子（其中
T
T
T 表示生成树，下同）
∑
T
∏
e
∈
T
p
e
∏
e
∉
T
(
1
−
p
e
)
\sum_T\prod_{e\in T}p_e\prod_{e\not\in T}(1-p_e)
T∑e∈T∏pee∈T∏(1−pe)
对比一下变元矩阵树定理的式子
∑
T
∏
e
∈
T
w
e
\sum_T\prod_{e\in T}w_e
T∑e∈T∏we
发现我们要求的式子的边权似乎不那么统一，考虑能不能转化一下。 设
D
=
∏
e
(
1
−
p
e
)
D=\prod_e(1-p_e)
D=∏e(1−pe)，那么可以将
D
D
D 从式子中提出来，式子变为
D
∑
T
∏
e
∈
T
p
e
1
−
p
e
D\sum_{T}\prod_{e\in T}\frac{p_e}{1-p_e}
DT∑e∈T∏1−pepe
于是只要设
e
e
e 的边权为
p
e
1
−
p
e
\frac{p_e}{1-p_e}
1−pepe 就是模板题了，这题可以说相当裸。 要注意当
p
e
p_e
pe 小于
e
p
s
\rm eps
eps 时直接取
e
p
s
\rm eps
eps，否则计算时会出现无穷，大于
1
−
e
p
s
1 - \rm eps
1−eps 同理。
代码
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std
;
const int maxn
= 55;
const double eps
= 1e-8;
int n
;
double t
,D
= 1,a
[maxn
][maxn
];
int read(){
int x
= 0;
char c
= getchar();
while(c
< '0' || c
> '9') c
= getchar();
while(c
>= '0' && c
<= '9') x
= x
* 10 + (c
^ 48),c
= getchar();
return x
;
}
double det(int n
){
double d
= 1;
for(int i
= 1; i
<= n
; i
++){
if(fabs(a
[i
][i
]) < eps
){
for(int j
= i
+ 1; j
<= n
; j
++){
if(fabs(a
[j
][i
]) > eps
){
d
*= -1;
for(int k
= i
; k
<= n
; k
++) swap(a
[i
][k
],a
[j
][k
]);
break;
}
}
}
for(int j
= i
+ 1; j
<= n
; j
++){
double t
= a
[j
][i
] / a
[i
][i
];
for(int k
= i
; k
<= n
; k
++) a
[j
][k
] -= t
* a
[i
][k
];
}
d
*= a
[i
][i
];
}
return d
;
}
int main(){
n
= read();
for(int i
= 1; i
<= n
; i
++)
for(int j
= 1; j
<= n
; j
++){
scanf("%lf",&t
);
if(t
< eps
) t
= eps
;
if(t
> 1 - eps
) t
= 1 - eps
;
if(i
> j
) D
*= 1 - t
;
if(i
!= j
) a
[i
][j
] = t
/ (t
- 1);
}
for(int i
= 1; i
<= n
; i
++)
for(int j
= 1; j
<= n
; j
++)
if(i
!= j
) a
[i
][i
] -= a
[i
][j
];
printf("%f\n",det(n
- 1) * D
);
return 0;
}
[SHOI2016]黑暗前的幻想乡
分析
直接容斥，然后对每一种情况重新建立矩阵，时间复杂度为
O
(
2
n
n
3
)
\Omicron(2^nn^3)
O(2nn3)
代码
#include <iostream>
#include <cstdio>
#include <vector>
using namespace std
;
const int maxn
= 20,mod
= 1e9 + 7;
int n
,m
,ans
,a
[maxn
][maxn
];
vector
<pair
<int,int> > e
[maxn
];
int read(){
int x
= 0;
char c
= getchar();
while(c
< '0' || c
> '9') c
= getchar();
while(c
>= '0' && c
<= '9') x
= x
* 10 + (c
^ 48),c
= getchar();
return x
;
}
inline int add(int x
,int y
){
if(x
+ y
< mod
) return x
+ y
;
else return x
+ y
- mod
;
}
inline int dec(int x
,int y
){
if(x
- y
>= 0) return x
- y
;
else return x
- y
+ mod
;
}
int qpow(int x
,int k
){
int d
= 1,t
= x
;
while(k
){
if(k
& 1) d
= 1ll * d
* t
% mod
;
t
= 1ll * t
* t
% mod
,k
>>= 1;
}
return d
;
}
int det(int n
){
int f
= 1,d
= 1;
for(int i
= 1; i
<= n
; i
++){
if(a
[i
][i
] == 0){
for(int j
= i
+ 1; j
<= n
; j
++){
if(a
[j
][i
]){
f
^= 1;
for(int k
= i
; k
<= n
; k
++) swap(a
[i
][k
],a
[j
][k
]);
break;
}
}
}
int inv
= qpow(a
[i
][i
],mod
- 2);
for(int j
= i
+ 1; j
<= n
; j
++){
int t
= 1ll * a
[j
][i
] * inv
% mod
;
for(int k
= i
; k
<= n
; k
++) a
[j
][k
] = dec(a
[j
][k
],1ll * t
* a
[i
][k
] % mod
);
}
d
= 1ll * d
* a
[i
][i
] % mod
;
}
return f
? d
: (mod
- d
);
}
int main(){
n
= read(),m
= 1 << (n
- 1);
for(int i
= 0; i
< n
- 1; i
++){
int t
= read();
while(t
--) e
[i
].push_back({read(),read()});
}
for(int i
= 0; i
< m
; i
++){
for(int j
= 1; j
<= n
; j
++)
for(int k
= 1; k
<= n
; k
++)
a
[j
][k
] = 0;
int cnt
= 0;
for(int j
= 0; j
< n
- 1; j
++){
if((i
>> j
& 1) == 0) continue;
cnt
++;
for(int k
= e
[j
].size() - 1; k
>= 0; k
--){
int u
= e
[j
][k
].first
,v
= e
[j
][k
].second
;
a
[u
][u
] ++,a
[v
][v
] ++;
a
[u
][v
] --,a
[v
][u
] --;
}
}
for(int j
= 1; j
<= n
; j
++)
for(int k
= 1; k
<= n
; k
++)
if(a
[j
][k
] < 0) a
[j
][k
] += mod
;
if((n
- cnt
) % 2) ans
= add(ans
,det(n
- 1));
else ans
= dec(ans
,det(n
- 1));
}
printf("%d\n",ans
);
return 0;
}
作业题
分析
今年省选的 Day2T3，现在回去看怎么感觉像是 Day2 最简单的题 ¿
先推一下这个式子
A
n
s
=
∑
T
∑
e
∈
T
w
e
gcd
e
∈
T
(
w
e
)
=
∑
T
∑
e
∈
T
w
e
∑
d
∣
gcd
e
∈
T
(
w
e
)
φ
(
d
)
=
∑
T
∑
e
∈
T
w
e
∑
∀
e
∈
T
,
d
∣
w
e
φ
(
d
)
=
∑
d
=
1
V
φ
(
d
)
(
∑
T
,
∀
e
∈
T
,
d
∣
w
e
∑
e
∈
T
w
e
)
\begin{aligned} Ans&=\sum_{T}\sum_{e\in T}w_e\gcd_{e\in T}(w_e)\\ &=\sum_{T}\sum_{e\in T}w_e\sum_{d\,|\gcd_{e\in T}(w_e)}\varphi(d)\\ &=\sum_{T}\sum_{e\in T}w_e\sum_{\forall e\in T,d|w_e}\varphi(d)\\ &=\sum_{d=1}^V\varphi(d)\left(\sum_{T,\atop \forall e\in T,d|w_e}\sum_{e\in T} w_e\right) \end{aligned}
Ans=T∑e∈T∑wee∈Tgcd(we)=T∑e∈T∑wed∣gcde∈T(we)∑φ(d)=T∑e∈T∑we∀e∈T,d∣we∑φ(d)=d=1∑Vφ(d)⎝⎜⎛∀e∈T,d∣weT,∑e∈T∑we⎠⎟⎞
其中，
V
V
V 表示值域最大值，下同。
先考虑右边的部分，如果没有
d
∣
w
e
d \,| \,w_e
d∣we 这个限制要怎么快速求？ 这是一个比较经典的套路，我们可以把每条边的权值设为一个一次多项式
1
+
w
x
1+wx
1+wx （其中
w
w
w 表示这条边的边权），那么每个生成树边权的和就转化成了生成树边权的积所得多项式的一次项系数，于是就可以直接上变元矩阵树定理了。
再看回原式子，加上限制以后怎么搞？ 如果我们对每个
d
d
d 重新加边，那么只加满足
d
∣
w
e
d\, | \,w_e
d∣we 的边就一定可以满足限制条件，就可以用上面的方法求了。但这样的时间复杂度为
O
(
V
n
3
)
\Omicron(Vn^3)
O(Vn3)，不可能通过，考虑优化。
容易发现当图不连通时对答案是一定没有贡献的，而如果我们加入的边不足
n
−
1
n-1
n−1 条，图一定不连通，于是我们可以只对边数大于
n
−
1
n-1
n−1 的
d
d
d 求行列式。因为对每条边来说，它最多只会被它边权的每个因数算一遍，而每次跑行列式都需要至少
n
−
1
n-1
n−1 条边，故复杂度上界为
O
(
t
m
n
−
1
n
3
)
=
O
(
t
n
4
)
\Omicron(\frac{tm}{n-1}n^3)=\Omicron(tn^4)
O(n−1tmn3)=O(tn4)（
m
m
m 是
n
2
n^2
n2 级别的），其中
t
t
t 表示值域中最大的因数个数。 打表可知
t
=
144
t=144
t=144，可以通过，且实际上远跑不满。
代码
#include <iostream>
#include <cstdio>
using namespace std
;
const int maxn
= 35,maxm
= 500,maxv
= 1.6e5,mod
= 998244353;
int n
,m
,mx
,ans
,tot
,u
[maxm
],v
[maxm
],w
[maxm
],vis
[maxv
],p
[maxv
],phi
[maxv
];
inline int add(int x
,int y
){
if(x
+ y
< mod
) return x
+ y
;
else return x
+ y
- mod
;
}
inline int dec(int x
,int y
){
if(x
- y
>= 0) return x
- y
;
else return x
- y
+ mod
;
}
int qpow(int x
,int k
){
int d
= 1,t
= x
;
while(k
){
if(k
& 1) d
= 1ll * d
* t
% mod
;
t
= 1ll * t
* t
% mod
,k
>>= 1;
}
return d
;
}
struct node
{
int x
,y
;
node
operator + (node p
){
return {add(x
,p
.x
),add(y
,p
.y
)};
}
node
operator - (node p
){
return {dec(x
,p
.x
),dec(y
,p
.y
)};
}
node
operator * (node p
){
return {1ll * x
* p
.x
% mod
,add(1ll * x
* p
.y
% mod
,1ll * y
* p
.x
% mod
)};
}
}a
[maxn
][maxn
];
int read(){
int x
= 0;
char c
= getchar();
while(c
< '0' || c
> '9') c
= getchar();
while(c
>= '0' && c
<= '9') x
= x
* 10 + (c
^ 48),c
= getchar();
return x
;
}
inline node
Inv(node p
){
int t
= qpow(p
.x
,mod
- 2);
t
= 1ll * t
* t
% mod
;
p
.y
= mod
- 1ll * p
.y
* t
% mod
,p
.x
= 1ll * p
.x
* t
% mod
;
return p
;
}
int det(int n
){
int f
= 1;
node d
= {1,0};
for(int i
= 1; i
<= n
; i
++){
if(!a
[i
][i
].x
){
bool flag
= 0;
for(int j
= i
+ 1; j
<= n
; j
++){
if(a
[j
][i
].x
){
f
^= 1,flag
= 1;
for(int k
= i
; k
<= n
; k
++) swap(a
[i
][k
],a
[j
][k
]);
break;
}
}
if(!flag
) return 0;
}
node inv
= Inv(a
[i
][i
]);
for(int j
= i
+ 1; j
<= n
; j
++){
node t
= a
[j
][i
] * inv
;
for(int k
= i
; k
<= n
; k
++) a
[j
][k
] = a
[j
][k
] - t
* a
[i
][k
];
}
d
= d
* a
[i
][i
];
}
return f
? d
.y
: mod
- d
.y
;
}
void Sieve(int n
){
phi
[1] = 1;
for(int d
= 2; d
<= n
; d
++){
if(!vis
[d
]) p
[++tot
] = d
,phi
[d
] = d
- 1;
for(int i
= 1; i
<= tot
&& p
[i
] * d
<= mx
; i
++){
int v
= p
[i
] * d
;
vis
[v
] = 1;
if(d
% p
[i
] == 0){
phi
[v
] = p
[i
] * phi
[d
];
break;
}
phi
[v
] = phi
[p
[i
]] * phi
[d
];
}
}
}
int main(){
n
= read(),m
= read();
for(int i
= 1; i
<= m
; i
++) u
[i
] = read(),v
[i
] = read(),w
[i
] = read(),mx
= max(mx
,w
[i
]);
Sieve(mx
);
for(int d
= 1; d
<= mx
; d
++){
int cnt
= 0;
for(int i
= 1; i
<= m
; i
++) if(w
[i
] % d
== 0) cnt
++;
if(cnt
< n
- 1) continue;
for(int i
= 1; i
<= n
; i
++)
for(int j
= 1; j
<= n
; j
++)
a
[i
][j
] = {0,0};
for(int i
= 1; i
<= m
; i
++)
if(w
[i
] % d
== 0){
a
[u
[i
]][v
[i
]] = {mod
- 1,mod
- w
[i
]};
a
[v
[i
]][u
[i
]] = {mod
- 1,mod
- w
[i
]};
a
[u
[i
]][u
[i
]] = a
[u
[i
]][u
[i
]] + (node
){1,w
[i
]};
a
[v
[i
]][v
[i
]] = a
[v
[i
]][v
[i
]] + (node
){1,w
[i
]};
}
ans
= add(ans
,1ll * phi
[d
] * det(n
- 1) % mod
);
}
printf("%d\n",ans
);
return 0;
}
