图论学习笔记II
xtx1092515503
2020-07-19 14:14:24
这里是图论学习笔记II。
将有关图上的树或是相近的算法归类到了本文中。
目前（可能）包括最小生成树、Kruskal重构树、Matrix-Tree、BEST定理（虽然这个大概不算图上建出的树，但是与Matrix-Tree类似那么就放过来罢）、虚树、圆方树等知识点。
# I.最小生成树
暂无。
# II.Kruskal重构树
原先的Kruskal重构树笔记已经被整合进图论学习笔记了。
Kruskal重构树，是用于求出有关一张图中，某点仅经过边权 $\leq$ 某个值 $v$ 的边所得到的子图的有关信息的工具。
但事实上，其应用还有更多。
我们先讲述其构造方法：
1. 将所有边按照边权递增排序。
2. 依次枚举每一条边。假如此时边的两个端点处于两个不同集合中，按照常规Kruskal算法，此时应该将该边加入MST，并在dsu上合并这两个点；然而，在Kruskal重构树算法中，我们将先新建一个点，然后将 $x,y$ 所在子树的根当作该新建点的两个儿子。同时，将该新建点的权值设作该边的边权。
具体而言，假设我们有三个点，以及三条边 $(1,2,1),(2,3,2),(1,3,4)$。下面我们构造重构树：
1. 第一条边 $(1,2,1)$，两个端点 $(1,2)$ 不位于同一子树中，故新建点 $4$，点权为 $1$，两个儿子为 $1$ 和 $2$。
2. 第二条边 $(2,3,2)$，两个端点 $(2,3)$ 不位于同一子树中，故新建点 $5$。此时，$2$ 所在子树的根是 $4$，$3$ 所在子树的根是 $3$，故连边 $5\rightarrow4,5\rightarrow3$，并将 $5$ 的点权设作 $2$。
3. 第三条边 $(1,3,4)$，两个端点 $(1,3)$ 已都位于 $5$ 的子树中，故不作任何操作。
就此我们构造出了Kruskal重构树。
它有什么性质呢？
1. 其成一个大根堆
因为我们加入的边权递增，故一个点的权值一定大于其两个儿子的权值，假如我们把原图中点的权值设作 $\infty$ 的话。
2. 原图中两点间所有路径中，路径上所有边的最大值最小的那一条路径的上述最大值为重构树上两点LCA的权值。
首先，路径上所有边的最大值最小的路径，一定是MST上的路径；而重构树上两点间路径唯一等价于MST上路径；又重构树有大根堆性质，故最大值一定出现在LCA处。
这也就给我们提供了截出到 $x$ 只经过权值不超过 $v$ 的点集的方法：在重构树上树上倍增找到最后一个权值不大于 $v$ 的 $x$ 的祖先 $y$，然后 $y$ 子树中所有节点即为该点集。因为 $x$ 到 $y$ 子树内任何点的LCA一定也在 $y$ 的子树内，$x$ 到 $y$ 子树外任何点的LCA一定也在 $y$ 的子树外，按照大根堆性质，前者一定在点集内，后者一定在点集外。
Kruskal重构树的复杂度为 $O(n\log n)$，来自于排序和**不能按秩合并，只能路径压缩**的冰茶姬。
## I.[Peaks](https://www.luogu.com.cn/problem/P4197)
可以被直接当作模板了。
树上倍增到需要的节点，这样所有合法的点集就变成了重构树上某节点的子树中所有节点，也即dfs序上一段区间。于是问题变为经典的区间 $k$ 大数问题，直接主席树带走。
时间复杂度 $O(n\log n)$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
int n,m,q,a[100100],dsu[200100],num,val[200100],anc[200100][20],L[200100],R[200100],rev[100100],cnt,rt[100100],tot;
struct node{
int x,y,z;
friend bool operator<(const node&u,const node&v){return u.z<v.z;}
void read(){scanf("%d%d%d",&x,&y,&z);}
}e[500100];
vector<int>v[200100],u;
int find(int x){return dsu[x]==x?x:dsu[x]=find(dsu[x]);}
void dfs(int x,int fa){
anc[x][0]=fa;
if(v[x].empty())L[x]=R[x]=++tot,rev[tot]=x;
else L[x]=tot+1,dfs(v[x][0],x),dfs(v[x][1],x),R[x]=tot;
}
int doubling(int x,int k){
for(int i=19;i>=0;i--)if(anc[x][i]&&val[anc[x][i]]<=k)x=anc[x][i];
return x;
}
struct SegTree{int lson,rson,sum;}seg[3201000];
#define mid ((l+r)>>1)
void modify(int &x,int y,int l,int r,int P){
if(l>P||r<P)return;
x=++cnt,seg[x]=seg[y],seg[x].sum++;
if(l!=r)modify(seg[x].lson,seg[y].lson,l,mid,P),modify(seg[x].rson,seg[y].rson,mid+1,r,P);
}
int query(int x,int y,int l,int r,int k){
if(l==r)return u[l-1];
if(seg[seg[x].lson].sum-seg[seg[y].lson].sum>=k)return query(seg[x].lson,seg[y].lson,l,mid,k);
else return query(seg[x].rson,seg[y].rson,mid+1,r,k-(seg[seg[x].lson].sum-seg[seg[y].lson].sum));
}
int main(){
scanf("%d%d%d",&n,&m,&q),num=n;
for(int i=1;i<=n;i++)scanf("%d",&a[i]),dsu[i]=i,u.push_back(a[i]);
for(int i=1;i<=m;i++)e[i].read();
sort(e+1,e+m+1);
for(int i=1;i<=m;i++){
int x=find(e[i].x),y=find(e[i].y);
if(x==y)continue;
num++,dsu[x]=dsu[y]=dsu[num]=num,val[num]=e[i].z;
v[num].push_back(x),v[num].push_back(y);
}
for(int i=1;i<=num;i++)if(dsu[i]==i)dfs(i,0);
sort(u.begin(),u.end()),u.resize(m=unique(u.begin(),u.end())-u.begin());
for(int i=1;i<=n;i++)a[i]=lower_bound(u.begin(),u.end(),a[i])-u.begin()+1;
for(int j=1;j<=19;j++)for(int i=1;i<=num;i++)anc[i][j]=anc[anc[i][j-1]][j-1];
for(int i=1;i<=n;i++)modify(rt[i],rt[i-1],1,m,a[rev[i]]);
for(int i=1,x,y,z;i<=q;i++){
scanf("%d%d%d",&x,&y,&z);
x=doubling(x,y);
// for(int j=L[x];j<=R[x];j++)printf("%d ",rev[j]);puts("");
if(R[x]-L[x]+1<z){puts("-1");continue;}
printf("%d\n",query(rt[R[x]],rt[L[x]-1],1,m,R[x]-L[x]+2-z));
}
return 0;
}
```
## II.[[NOI2018] 归程](https://www.luogu.com.cn/problem/P4768)
建出重构树，则问题转变为子树中所有点到 $1$ 的距离的最小值。于是预先 Dijkstra 跑一遍求出所有点到 $1$ 的最短路径，然后再遍历重构树求出子树最短路径，然后回答询问时就直接树上倍增即可。
时间复杂度 $O(n\log n)$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
int T,n,m,ch[400100][2],dis[400100],num,anc[400100][20],dsu[400100],val[400100],Q,K,S,lasans;
vector<pair<int,int> >v[200100];
struct node{
int x,y,z;
friend bool operator<(const node&u,const node&v){return u.z>v.z;}
node(){}
node(int X,int Y,int Z){x=X,y=Y,z=Z;}
}e[400100];
priority_queue<pair<int,int> >q;
void Dijkstra(){
q.push(make_pair(0,1)),dis[1]=0;
while(!q.empty()){
int x=q.top().second;q.pop();
for(auto y:v[x])if(dis[y.first]>dis[x]+y.second)dis[y.first]=dis[x]+y.second,q.push(make_pair(-dis[y.first],y.first));
}
}
int find(int x){return dsu[x]==x?x:dsu[x]=find(dsu[x]);}
void dfs(int x,int fa){anc[x][0]=fa;if(x>n)dfs(ch[x][0],x),dfs(ch[x][1],x),dis[x]=min(dis[ch[x][0]],dis[ch[x][1]]);}
int doubling(int x,int k){for(int i=19;i>=0;i--)if(anc[x][i]&&val[anc[x][i]]>k)x=anc[x][i];return x;}
int main(){
scanf("%d",&T);
while(T--){
scanf("%d%d",&n,&m),memset(dis,0x7f,sizeof(dis)),num=n;
for(int i=1,x,y,z,w;i<=m;i++)scanf("%d%d%d%d",&x,&y,&z,&w),v[x].push_back(make_pair(y,z)),v[y].push_back(make_pair(x,z)),e[i]=node(x,y,w);
Dijkstra();
sort(e+1,e+m+1);
for(int i=1;i<=n;i++)dsu[i]=i;
for(int i=1;i<=m;i++){
int x=find(e[i].x),y=find(e[i].y);
if(x==y)continue;
num++,dsu[x]=dsu[y]=dsu[num]=num,val[num]=e[i].z;
ch[num][0]=x,ch[num][1]=y;
}
for(int i=1;i<=n;i++)v[i].clear();
dfs(num,0);
for(int j=1;j<=19;j++)for(int i=1;i<=num;i++)anc[i][j]=anc[anc[i][j-1]][j-1];
scanf("%d%d%d",&Q,&K,&S),lasans=0;
for(int i=1,x,y;i<=Q;i++){
scanf("%d%d",&x,&y),x=(x+K*lasans-1)%n+1,y=(y+K*lasans)%(S+1);
printf("%d\n",lasans=dis[doubling(x,y)]);
}
}
return 0;
}
```
## III.[[ARC098D] Donation](https://www.luogu.com.cn/problem/AT4144)
Kruskal重构树还可以用来优化DP哦~~
结论1.最优情形下，当我们向一个节点”捐赠“后，我们不会再次访问该节点。
明显，如果再次访问，则在再次访问时再捐献一定不更劣。
我们定义 $c_x=\max(0,a_x-b_x)$，则约束等价于离开 $x$ 点时剩余钱数不小于 $c_x$。
明显我们肯定是希望从大到小地依次给每个 $c_x$ 捐赠，因为反正离开 $x$ 时无论如何都会剩下 $c_x$，如果先给大的捐赠，这剩下的钱说不定就在之后捐出去了；但是如果后给大的 $c_x$ 捐赠，这剩下的钱就有可能烂在手里捐不出去了。
但是，这只是希望，并不代表我们真的做得到，因为在从一个 $c_x$ 到另一个 $c_y$ 的过程中，我们说不定会路过一个拥有更大 $c_z$ 的点 $z$ 而导致需要更多的入场费。
但是，我们发现，从 $c_x$ 到 $c_y$ 的本质是让路径最大值最小，同Kruskal重构树的本质一致。
为了让重构树也可以应用于点权，我们将连接两点 $(u,v)$ 的边权设作 $\max(c_u,c_v)$。这样，我们便可建出重构树来。
设对于重构树上一个非原图节点，其 $c_x$ 为建出重构树时的边权。
设 $f_x$ 表示若我们给重构树上节点 $x$ 的子树内所有节点全都被捐赠过，所需要的最小钱数。
明显，根据结论1，我们捐赠的路径一定是 $\text{一个儿子}\rightarrow x\rightarrow\text{另一个儿子}$，不可能出现两个儿子子树内的节点全捐完后再来捐 $x$ 的情形，因为从一个儿子中节点走到另一个儿子中节点的最优路径必经过 $x$，故该路径需要的最大剩余钱数也是 $c_x$，而在一个儿子内部乱跑的最大钱数一定不大于 $c_x$。
因为我们离开 $x$ 时剩余钱数一定不小于 $c_x$，而凭借着这么多钱在第一个儿子里一定可以随便乱跑而不必多带钱，因此设 $s_x$ 表示 $x$ 节点子树中所有节点的 $b_x$ 之和，则从第一个儿子 $l$ 出来时所需费用一定仅需 $s_l$ 即可。然后，处理 $x$ 本身以及另一个儿子 $r$ 所需费用，为 $\max(f_r,c_x)$，因为要保证经过 $x$ 的过路限制以及付给 $r$ 子树中点的费用。于是总代价为 $s_l+\max(f_r,c_x)$。同理，$r\rightarrow x\rightarrow l$ 也是一条合法路径，也可以贡献为备选答案。
对于叶子节点，其 $f$ 直接设作 $b_x+c_x$ 即可。
时间复杂度 $O(n\log n)$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int n,m,b[100100],c[100100],dsu[200100],num,val[200100],ch[400100][2];
ll f[200100],s[200100];
struct node{
int x,y,z;
friend bool operator<(const node&u,const node&v){return u.z<v.z;}
void read(){scanf("%d%d",&x,&y),z=max(c[x],c[y]);}
}e[100100];
int find(int x){return dsu[x]==x?x:dsu[x]=find(dsu[x]);}
void dfs(int x){
if(x<=n)f[x]=b[x]+c[x],s[x]=b[x];
else dfs(ch[x][0]),dfs(ch[x][1]),f[x]=min(s[ch[x][0]]+max(1ll*val[x],f[ch[x][1]]),s[ch[x][1]]+max(1ll*val[x],f[ch[x][0]])),s[x]=s[ch[x][0]]+s[ch[x][1]];
}
int main(){
scanf("%d%d",&n,&m),num=n;
for(int i=1,A,B;i<=n;i++)scanf("%d%d",&A,&B),b[i]=B,c[i]=max(0,A-B);
for(int i=1;i<=m;i++)e[i].read();
sort(e+1,e+m+1);
for(int i=1;i<=n;i++)dsu[i]=i;
for(int i=1;i<=m;i++){
int x=find(e[i].x),y=find(e[i].y);
if(x==y)continue;
num++,dsu[x]=dsu[y]=dsu[num]=num,val[num]=e[i].z;
ch[num][0]=x,ch[num][1]=y;
}
dfs(num);
printf("%lld\n",f[num]);
return 0;
}
```
## IV.[[IOI2018] werewolf 狼人](https://www.luogu.com.cn/problem/P4899)
明显在重构树上，所有人能到的点集是dfs序上一段区间，狼能到的点集也是一段区间（注意这里是两棵不同的树，人树的边权是两边的 $\min$，且求最小重构树；而狼树的边权是两边 $\min$，求的是最大重构树）。而我们要判断两段区间里是否有相同的数。
于是对于任意一棵树的dfs序上每个位置，求出其元素在另一棵树上对应的dfs序，然后问题就转换为二维数点问题，判断矩形内是否有元素。简单离线后BIT即可。
但是IOI赛场上是交互式，因此强制在线，故要用主席树。
但是LG上不需要。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
int n,m,q,a[200100];
struct EDGE{
int u,v,w;
EDGE(){}
EDGE(int U,int V,int W){u=U,v=V,w=W;}
};
struct MST{
EDGE e[400100];
vector<int>v[400100];
int dsu[400100],tot,val[400100],anc[400100][20],L[400100],R[400100],cnt,dfn[200100];
int find(int x){return dsu[x]==x?x:dsu[x]=find(dsu[x]);}
void dfs(int x){
L[x]=cnt+1;if(x<=n)dfn[++cnt]=x;
for(auto y:v[x])anc[y][0]=x,dfs(y);
R[x]=cnt;
}
int doubling(int x,int lim){for(int i=19;i>=0;i--)if(anc[x][i]&&val[anc[x][i]]<=lim)x=anc[x][i];return x;}
void Kruskal(){
tot=n;
sort(e+1,e+m+1,[](EDGE x,EDGE y){return x.w<y.w;});
for(int i=1;i<=n;i++)dsu[i]=i;
for(int i=1;i<=m;i++){
int X=find(e[i].u),Y=find(e[i].v);
if(X==Y)continue;
val[++tot]=e[i].w;
dsu[tot]=dsu[X]=dsu[Y]=tot;
v[tot].push_back(X),v[tot].push_back(Y);
}
// for(int i=1;i<=tot;i++)printf("%d ",val[i]);puts("");
dfs(tot);
for(int j=1;j<=19;j++)for(int i=1;i<=tot;i++)anc[i][j]=anc[anc[i][j-1]][j-1];
// for(int i=1;i<=n;i++)printf("%d ",dfn[i]);puts("");
// for(int i=1;i<=tot;i++)printf("[%d,%d]",L[i],R[i]);puts("");
}
}mnst,mxst;
int pos[200100],val[200100],res[200100],t[200100];
void ADD(int x){while(x<=n)t[x]++,x+=x&-x;}
int SUM(int x){int ret=0;while(x)ret+=t[x],x-=x&-x;return ret;}
vector<pair<int,int> >v[200100];
int main(){
scanf("%d%d%d",&n,&m,&q);
for(int i=1,x,y;i<=m;i++)scanf("%d%d",&x,&y),x++,y++,mnst.e[i]=EDGE(x,y,max(x,y)),mxst.e[i]=EDGE(x,y,-min(x,y));
mnst.Kruskal(),mxst.Kruskal();
for(int i=1;i<=n;i++)pos[mnst.dfn[i]]=i;
for(int i=1;i<=n;i++)val[i]=pos[mxst.dfn[i]];
for(int i=1,s,t,l,r,x;i<=q;i++){
scanf("%d%d%d%d",&s,&t,&l,&r),s++,t++,l++,r++;
x=mxst.doubling(s,-l);int L1=mxst.L[x]-1,R1=mxst.R[x];//printf("%d:[%d,%d]\n",x,L1+1,R1);
x=mnst.doubling(t,r);int L2=mnst.L[x]-1,R2=mnst.R[x];//printf("%d:[%d,%d]\n",x,L2+1,R2);
v[L1].push_back(make_pair(L2,i));
v[L1].push_back(make_pair(R2,-i));
v[R1].push_back(make_pair(L2,-i));
v[R1].push_back(make_pair(R2,i));
}
for(int i=1;i<=n;i++){
ADD(val[i]);
for(auto j:v[i])if(j.second>0)res[j.second]+=SUM(j.first);else res[-j.second]-=SUM(j.first);
}
for(int i=1;i<=q;i++)printf("%d\n",!!res[i]);
return 0;
}
```
# III.Matrix Tree定理
线性代数基础详见[线性代数学习笔记](https://www.luogu.com.cn/blog/Troverld/linear-algebra)。
矩阵树定理是用来作**生成树计数**的好东西。其定理主要表述如下：
设 $\mathbb{G}$ 为一无向图。则其无向无根生成树的数量，可以被如下算法求得：
设矩阵 $S_1$，其中第 $i$ 行第 $j$ 列的数是 $(i,j)$ 间无向边数。
设矩阵 $S_2$，其是一对角矩阵（即只有对角线有数的矩阵），其中第 $i$ 行第 $i$ 列的数是 $i$ 的度数。
设矩阵 $S_3=S_2-S_1$
则其生成树数量，为 $S_3$ 删去任意一行一列所得到的矩阵的行列式大小。
~~证明？OI需要证明吗？~~
------
其有二变种。一个是求有根外向生成树数，另一个是求有根内向生成树数。
这两种的 $S_1$ 都是从 $i$ 到 $j$ 的有向边数。但是，外向树的 $S_2$ 是入度数，内向树的 $S_2$ 是出度数。
同时，这时就不能随便删行删列了——必须删去根所在的那一行一列。
------
## I. [【模板】Matrix-Tree 定理](https://www.luogu.com.cn/problem/P6178)
就是模板啦。但是，这题有边权，怎样处理？
我们考虑如果将一条带权边拆作（权值）条无权边，则任何一棵带权边的生成树，其中的一条边，就可以换成任何一条无权边，一共（权值）条，是等效的。依据乘法原理，我们将每条带权边看作无权边，总权值仍然是符合实际意义的。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
int n,m,T,a[310][310];
int ksm(int x,int y=mod-2){
int z=1;
for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
return z;
}
int Determinant(){
int lam=(n&1?1:-1);
for(int i=2;i<=n;i++){
int pos=i;
for(;pos<n;pos++)if(a[pos][i])break;
if(pos>n)return 0;
if(pos!=i)swap(a[pos],a[i]),lam=-lam;
for(int j=i+1;j<=n;j++){
int del=1ll*a[j][i]*ksm(a[i][i])%mod;
for(int k=i;k<=n;k++)(a[j][k]+=mod-1ll*a[i][k]*del%mod)%=mod;
}
}
(lam+=mod)%=mod;
for(int i=2;i<=n;i++)lam=1ll*lam*a[i][i]%mod;
return lam;
}
int main(){
scanf("%d%d%d",&n,&m,&T);
for(int i=1,x,y,z;i<=m;i++){
scanf("%d%d%d",&x,&y,&z);
if(x==y)continue;
if(T==0)(a[x][y]+=z)%=mod,(a[y][x]+=z)%=mod,(a[x][x]+=mod-z)%=mod,(a[y][y]+=mod-z)%=mod;
else (a[x][y]+=z)%=mod,(a[y][y]+=mod-z)%=mod;
}
printf("%d\n",Determinant());
return 0;
}
```
## II.[[SHOI2016]黑暗前的幻想乡](https://www.luogu.com.cn/problem/P4336)
考虑容斥套Matrix-Tree。我们考虑全部 $n-1$ 个公司全部可选的生成树数量，减去选 $n-2$ 个公司的数量，加上 $n-3$ 个……于是只需枚举哪些公司然后容斥一下即可。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
int n,res;
vector<pair<int,int> >v[20];
int a[20][20];
int ksm(int x,int y=mod-2){
int z=1;
for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
return z;
}
int Determinant(){
int lam=-1;
// for(int i=1;i<n;i++){for(int j=1;j<n;j++)printf("%d ",a[j][i]);puts("");}
for(int i=1;i<n;i++){
int j=i;
for(;j<n;j++)if(a[j][i])break;
if(j>=n)return 0;
if(i!=j)lam=-lam,swap(a[i],a[j]);
for(j=i+1;j<n;j++){
int del=1ll*a[j][i]*ksm(a[i][i])%mod;
for(int k=i;k<n;k++)(a[j][k]+=mod-1ll*a[i][k]*del%mod)%=mod;
}
}
// for(int i=1;i<n;i++){for(int j=1;j<n;j++)printf("%d ",a[j][i]);puts("");}
(lam+=mod)%=mod;
for(int i=1;i<n;i++)lam=1ll*lam*a[i][i]%mod;
return lam;
}
int main(){
scanf("%d",&n);
for(int i=0,z,x,y;i<n-1;i++){
scanf("%d",&z);
while(z--)scanf("%d%d",&x,&y),v[i].push_back(make_pair(x,y));
}
for(int i=0;i<(1<<(n-1));i++){
memset(a,0,sizeof(a));
for(int j=0;j<n-1;j++){
if(i&(1<<j))continue;
for(auto p:v[j]){
(++a[p.first][p.first])%=mod;
(++a[p.second][p.second])%=mod;
(a[p.first][p.second]+=mod-1)%=mod;
(a[p.second][p.first]+=mod-1)%=mod;
}
}
(res+=1ll*(__builtin_popcount(i)&1?1:mod-1)*Determinant()%mod)%=mod;
}
printf("%d\n",res);
return 0;
}
```
## III.[TopCoder13369-TreeDistance](https://vjudge.net/problem/TopCoder-13369)
这里是本题的矩阵树定理解法。DP解法详见重题——[DP刷题笔记III](https://www.luogu.com.cn/blog/Troverld/dp-shua-ti-bi-ji-iii)，CXX.[CF917D Stranger Trees](https://www.luogu.com.cn/problem/CF917D)。
我们考虑为所有树上的边赋权为 $1$，完全图中不在树上的边赋权为 $x$。对其求矩阵树定理。明显，其形式是一 $n-1$ 阶矩阵的行列式，也即一关于 $x$ 的 $n-1$ 次多项式。则，多项式上常数项系数即为出现 $0$ 条新边的生成树数，一次项系数即为出现 $1$ 条新边的生成树数……然后对系数求一个前缀和，即是最多出现 $k$ 条新边的生成树数。
但问题来了，没人教过我们应该如何对一个有未知项的矩阵求行列式呀（具体是因为消元的时候不知道应该怎样判定“$0$”）？
没关系。按照我们之前的结论，其行列式是一 $n-1$ 次多项式。那我们直接挑 $n$ 个固定的 $x$ 各跑一遍矩阵树，然后拉格朗日插一下就插出多项式了。当然，如果你嫌拉格朗日还要手工进行多项式乘除法太麻烦，也可以直接上高斯消元插多项式。
时间复杂度 $O(n^4)$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
class TreeDistance{
private:
int n,m,a[52][52],g[52][52],f[52];
int ksm(int x,int y=mod-2){
int z=1;
for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
return z;
}
int Determinant(){
int lam=(n&1?1:-1);
// for(int i=1;i<n;i++){for(int j=1;j<n;j++)printf("%d ",a[j][i]);puts("");}
for(int i=1;i<n;i++){
int j=i;
for(;j<n;j++)if(a[j][i])break;
if(j>=n)return 0;
if(i!=j)lam=-lam,swap(a[i],a[j]);
for(j=i+1;j<n;j++){
int del=1ll*a[j][i]*ksm(a[i][i])%mod;
for(int k=i;k<n;k++)(a[j][k]+=mod-1ll*a[i][k]*del%mod)%=mod;
}
}
// for(int i=1;i<n;i++){for(int j=1;j<n;j++)printf("%d ",a[j][i]);puts("");}
(lam+=mod)%=mod;
for(int i=1;i<n;i++)lam=1ll*lam*a[i][i]%mod;
return lam;
}
void Gauss(){
memset(a,0,sizeof(a));
for(int i=1;i<=n;i++){
for(int j=1,k=1;j<=n;j++,k=1ll*k*i%mod)a[i][j]=k;
a[i][n+1]=f[i];
}
// for(int i=1;i<=n;i++){for(int j=1;j<=n+1;j++)printf("%d ",a[i][j]);puts("");}
for(int i=1;i<=n;i++){
int j=i;
for(;j<=n;j++)if(a[j][i])break;
if(i!=j)swap(a[i],a[j]);
for(j=1;j<=n;j++){
if(j==i)continue;
int del=1ll*a[j][i]*ksm(a[i][i])%mod;
for(int k=i;k<=n+1;k++)(a[j][k]+=mod-1ll*a[i][k]*del%mod)%=mod;
}
}
}
public:
int countTrees(vector<int>p,int K){
m=K,n=p.size()+1;
memset(g,0,sizeof(g));
for(int i=2;i<=n;i++)g[p[i-2]+1][i]=g[i][p[i-2]+1]=1;
// for(int i=1;i<=n;i++){for(int j=1;j<=n;j++)printf("%d ",g[i][j]);puts("");}
for(int i=1;i<=n;i++){
memset(a,0,sizeof(a));
for(int j=1;j<=n;j++)for(int k=1;k<=n;k++){
if(j==k)continue;
if(g[j][k])(++a[j][k])%=mod,(a[k][k]+=mod-1)%=mod;
else (a[j][k]+=i)%=mod,(a[k][k]+=mod-i)%=mod;
}
f[i]=Determinant();
// printf("%d\n",f[i]);
}
Gauss();
int res=0;
for(int i=1;i<=min(n,m+1);i++)(res+=1ll*a[i][n+1]*ksm(a[i][i])%mod)%=mod;
return res;
}
}my;
/*
int main(){
printf("%d\n",my.countTrees({0,0},1));
printf("%d\n",my.countTrees({0,1,2,2,0},50));
}
*/
```
## IV.[[HEOI2015]小 Z 的房间](https://www.luogu.com.cn/problem/P4111)
首先，这题实际上很板子。但是，目的是通过本题介绍一下模数非质数时的高斯消元。
模数非质数时，无法求逆元；但是要说如何把一个东西消成 $0$ 的话，可以想到**辗转相除法**，即$\gcd$时要用的那东西。直接上辗转相除，可以被证明复杂度仍是 $O(n^3)$ 的。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int mod=1e9;
int n,m,id[10][10],p,a[100][100];
char s[10][10];
int Determinant(){
int lam=(p&1?1:-1);
for(int i=1;i<p;i++){
int j=i;
for(;j<p;j++)if(a[j][i])break;
if(j>=p)return 0;
swap(a[i],a[j]),lam=-lam;
for(j=i+1;j<p;j++)while(a[j][i]){
int del=a[i][i]/a[j][i];
for(int k=i;k<p;k++)(a[i][k]+=mod-1ll*a[j][k]*del%mod)%=mod;
swap(a[i],a[j]),lam=-lam;
}
}
(lam+=mod)%=mod;
for(int i=1;i<p;i++)lam=1ll*lam*a[i][i]%mod;
return lam;
}
int main(){
scanf("%d%d",&n,&m);
for(int i=0;i<n;i++)scanf("%s",s[i]);
for(int i=0;i<n;i++)for(int j=0;j<m;j++)if(s[i][j]=='.')id[i][j]=++p;
for(int i=0;i<n;i++)for(int j=0;j<m;j++){
if(i&&id[i][j]&&id[i-1][j])a[id[i][j]][id[i-1][j]]=a[id[i-1][j]][id[i][j]]=mod-1,a[id[i][j]][id[i][j]]++,a[id[i-1][j]][id[i-1][j]]++;
if(j&&id[i][j]&&id[i][j-1])a[id[i][j]][id[i][j-1]]=a[id[i][j-1]][id[i][j]]=mod-1,a[id[i][j]][id[i][j]]++,a[id[i][j-1]][id[i][j-1]]++;
}
// for(int i=1;i<=p;i++){for(int j=1;j<=p;j++)printf("%d ",a[i][j]);puts("");}
printf("%d\n",Determinant());
return 0;
}
```
## V.[[SDOI2014]重建](https://www.luogu.com.cn/problem/P3317)
这里就涉及到矩阵树定理的本质了——你费尽心思，又是建矩阵又是求行列式的，究竟是为了什么呢？
矩阵树定理的本质，是求
$$\sum\limits_{\mathbb{T}}\prod\limits_{e_i\in\mathbb{T}}e_i$$
而本题要求的，是
$$\sum\limits_{\mathbb{T}}\prod\limits_{e_i\in\mathbb{T}}e_i\prod\limits_{e_i\notin\mathbb{T}}(1-e_i)$$
考虑提出去一个 $\prod(1-e_i)$，得到
$$\prod(1-e_i)\sum\limits_{\mathbb{T}}\prod\limits_{e_i\in\mathbb{T}}\dfrac{e_i}{1-e_i}$$
这里就可以直接上矩阵树了。
但是，如果我们发现有边权恰为 $1$，应该咋办呢？
发现，恰为 $1$ 的实际意义是两点间必有边，故使用冰茶姬将两者并一块即可。同时，因为还要保证是一棵树，冰茶姬维护的同时还要保证无环。这里的无环包括两种意义的无环，一是 $1$ 边自身不能成环，二是 $1$ 边构成的连通块中不能出现其它边（即，其它边的选择都必须是 $1-e_i$）。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const double eps=1e-8;
int n,dsu[60],val[60],m;
double a[60][60],g[60][60],tmp,all=1;
double Determinant(){
double ret=1;
for(int i=1;i<m;i++){
int j=i,k=i;
for(;j<m;j++)if(fabs(a[j][i])>fabs(a[k][i]))k=j;
if(i!=k)ret=-ret,swap(a[i],a[k]);
if(fabs(a[i][i])<eps)return 0;
for(j=i+1;j<m;j++){
double del=a[j][i]/a[i][i];
for(k=i;k<m;k++)a[j][k]-=del*a[i][k];
}
ret*=a[i][i];
}
return ret;
}
int find(int x){return dsu[x]==x?x:dsu[x]=find(dsu[x]);}
bool merge(int x,int y){
x=find(x),y=find(y);
if(x==y)return false;
dsu[y]=x;
return true;
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)dsu[i]=i;
for(int i=1;i<=n;i++)for(int j=1;j<=n;j++){
scanf("%lf",&g[i][j]);
if(i<=j)continue;
if(fabs(1-g[i][j])>eps){all*=(1-g[i][j]);continue;}
if(!merge(i,j)){puts("0");return 0;}//100% have circle,can't be a tree
}
for(int i=1;i<=n;i++)if(dsu[i]==i)val[i]=++m;
for(int i=1;i<=n;i++)for(int j=1;j<=n;j++){
if(find(i)==find(j))continue;
int I=val[find(i)],J=val[find(j)];
a[I][J]-=g[i][j]/(1-g[i][j]);
a[I][I]+=g[i][j]/(1-g[i][j]);
}
printf("%.10lf\n",all*Determinant());
return 0;
}
```
## VI.[[JSOI2008]最小生成树计数](https://www.luogu.com.cn/problem/P4208)
我觉得是至今为止刷过的所有矩阵树中最棒的一个。
首先，有一个性质是很明显的：在一张图的所有最小生成树中，其每个权值的边的数量都是固定的。不然，假如有两棵MST，其中存在一个权值，使得该权值的边在两棵树中数量不同，则一定存在一边权和更小的生成树。
有了这个性质，我们就可以再推出一个更强的性质——其中任意一棵MST中，某一权值的边所联通的点集都是相同的，都是连完MST上其他权值的边后尚未联通的点集。
很明显，第一条性质是第二条性质的子集（固定的点集中边的数量也是固定的），故我们只需简单证明一下第二条性质即可。
考虑Kruskal算法的过程。考虑生成树中最大权值的所有边，依据Kruskal算法，在我们未添加该权值的任何边之前，其他边所联通的集合是固定的（假如并非固定，即存在两种不同的方案使得一对点在一种方案上联通而另一种不连通，则我们完全可以取两种方案上所有边的某种并集使得总权值和更小）——这也意味着该权值的所有边联通的集合是固定的。则我们现在将该集合看做全部用最大权值的边联通，然后按照此种方法再去分析次大权值的所有边，即可简单得到第二条性质。
于是我们现在就直接套矩阵树即可——因为每一种边权所联通的集合是固定的，故我们对该集合（当然，用dsu缩去不由该边权的点所联通的点集后）以及该边权的所有边直接上一个矩阵树计算联通该集合的方案数。最后，因为每种边权的集合间独立，上乘法原理将所有边权的方案乘一块即可。
时间复杂度 $O(n^3+nm+m\log m)$，因为最小生成树上的边数固定是 $n-1$，故大部分操作都没有求行列式的 $O(n^3)$ 大。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int mod=31011;
int n,m,dsu[110],a[110][110],id[110],p,ret=1;
int find(int x){return dsu[x]==x?x:dsu[x]=find(dsu[x]);}
bool merge(int x,int y){x=find(x),y=find(y);if(y==x)return false;dsu[y]=x;return true;}
struct node{
int x,y,z;
friend bool operator<(const node&u,const node&v){return u.z<v.z;}
}e[1010];
vector<int>v;
int Determinant(){
for(int i=1;i<=p;i++)for(int j=1;j<=p;j++)a[i][j]%=mod,(a[i][j]+=mod)%=mod;
int lam=1;
for(int i=1;i<p;i++){
int j=i;
for(;j<p;j++)if(a[j][i])break;
if(j>=p)return 0;
if(i!=j)swap(a[i],a[j]),lam=-lam;
for(j=i+1;j<p;j++)while(a[j][i]){
int del=a[i][i]/a[j][i];
for(int k=i;k<p;k++)(a[i][k]+=mod-1ll*del*a[j][k]%mod)%=mod;
swap(a[i],a[j]),lam=-lam;
}
}
(lam+=mod)%=mod;
for(int i=1;i<p;i++)lam=1ll*lam*a[i][i]%mod;
for(int i=1;i<=p;i++)for(int j=1;j<=p;j++)a[i][j]=0;
return lam;
}
int main(){
scanf("%d%d",&n,&m);
for(int i=1;i<=m;i++)scanf("%d%d%d",&e[i].x,&e[i].y,&e[i].z);
for(int i=1;i<=n;i++)dsu[i]=i;
sort(e+1,e+m+1);
for(int i=1;i<=m;i++)if(merge(e[i].x,e[i].y))v.push_back(i);
if(v.size()<n-1){puts("0");return 0;}
for(int i=0,j=0,I=1,J=1;i<v.size();i=j,I=J){
while(j<v.size()&&e[v[i]].z==e[v[j]].z)j++;
for(int k=1;k<=n;k++)dsu[k]=k;
for(int k=0;k<v.size();k++)if(k<i||k>=j)merge(e[v[k]].x,e[v[k]].y);
p=0;for(int k=1;k<=n;k++)if(dsu[k]==k)id[k]=++p;
while(e[I].z!=e[v[i]].z)I++;
for(J=I;J<=m&&e[J].z==e[I].z;J++){
int X=id[find(e[J].x)],Y=id[find(e[J].y)];
a[X][Y]--,a[Y][X]--;
a[X][X]++,a[Y][Y]++;
}
ret=1ll*ret*Determinant()%mod;
}
printf("%d\n",ret);
return 0;
}
```
## VII.[[省选联考 2020 A 卷] 作业题](https://www.luogu.com.cn/problem/P6624)
考虑式子
$$\sum\limits_{\mathbb{T}}\sum\limits_{e_i\in\mathbb{T}}e_i\gcd\limits_{e_i\in\mathbb{T}}\{e_i\}$$
对 $\gcd$ 一类东西，常规思路是[莫比乌斯反演](https://www.luogu.com.cn/blog/Troverld/mu-bi-wu-si-fan-yan-gan-xing-xia-che)掉。此处我们考虑欧拉反演。
$$\sum\limits_{\mathbb{T}}\sum\limits_{e_i\in\mathbb{T}}e_i\sum\limits_{d|e_i\in\mathbb{T}}\varphi(d)$$
把 $d$ 扔到前面去
$$\sum\limits_d\varphi(d)\sum\limits_{\mathbb{T}}[d|e_i\in\mathbb{T}]\sum\limits_{e_i\in\mathbb{T}}e_i$$
我们考虑就枚举这个 $d$，这样就只需考虑所有 $d|e_i$ 的 $e_i$ 构成的子图即可。于是问题转换为求 $\sum\limits_{\mathbb{T}}\sum\limits_{e_i\in\mathbb{T}}e_i$。
矩阵树解决的是后面为 $\prod$ 类型的问题，但此处是 $\sum$。一种方法是通过 $\exp$ 操作，让 $e_i$ 变成 $3^{e_i}$（因为 $3$ 是 $998244353$ 的原根），这样直接上矩阵树就得到了 $3^{\sum e_i}$。通过BSGS即可求出 $\sum e_i$。但这样复杂度多一个 $\sqrt{998244353}$，大概三万多，顶的上求那个行列式的复杂度了——况且2020省选还不给用 `unordered_map`，还得手写哈希表。
所以我们考虑其它方法。一种naive的方法是枚举某条边处在多少棵生成树中（可以通过将这条边所连接的两个点缩点后，将其他边的权值全部赋作 $1$ 然后跑矩阵树得到）然后用次数乘上边权即可。但这样复杂度是 $O(mn^3)$，不可能过。
但我们有一种神奇的做法——将矩阵中所有东西全都变成**一次函数**，若一条边的边权为 $e$，其将被变成函数 $ex+1$。
这样，计算此矩阵的矩阵树（按照行列式的定义，其结果是一个 $n-1$ 阶多项式），会发现其第 $k$ 次项上系数，意味着有 $k$ 条边贡献了 $ex$，$n-k-1$ 条边贡献了 $1$——换句话说，就是所有贡献了 $x$ 的项的边权之积。则我们考虑其一次项系数，其意义是只有一条边贡献了 $ex$，其它边的贡献都是 $1$——我们会发现此描述的实际意义就是边 $e$ 所在的生成树数量。故我们只需要求出此矩阵的Matrix-Tree的一次项即可。
求矩阵树就是求行列式。求行列式就是消元。消元只需实现一次项矩阵里函数的加减乘除即可——自然，结果只需保留常数项和一次项。加减乘都好办，唯独除难办。
考虑我们求出 $\dfrac{k_1x+b_1}{k_2x+b_2}$。手动模拟多项式求逆可得结果为 $\dfrac{k_1b_2-k_2b_1}{(b_2)^2}x+\dfrac{b_1}{b_2}$。
则直接对元素为一次函数的矩阵求行列式即可。
复杂度最多为 $m\sqrt{e}n^3$——因为一条边的不同因数数量是 $O(\sqrt{e})$ 级别的，故所有边的不同因数的并集是 $O(m\sqrt{e})$ 的。或许也能卡过。但是，观察到必须有至少 $n-1$ 条合法的边才有可能出现合法生成树，故我们可以只判断能出现生成树的权值（当然，用并查集判断），复杂度优化为上界 $\dfrac{n^3m\sqrt{e}}{n}=n^2m\sqrt{e}=n^4\sqrt{e}$。但这只是理论上界，实际应用起来跑得飞快。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
int n,m,lim,u[910],v[910],w[910],dsu[910],phi[200100],pri[200100],res;
int find(int x){return dsu[x]==x?x:dsu[x]=find(dsu[x]);}
void merge(int x,int y){x=find(x),y=find(y);if(x!=y)dsu[y]=x;}
int ksm(int x,int y=mod-2){
int z=1;
for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
return z;
}
struct Poly{
int k,b;//y=kx+b
Poly(){k=b=0;}
Poly(int x){k=x,b=1;}
Poly(int x,int y){k=x,b=y;}
friend Poly operator+(const Poly&x,const Poly&y){return Poly((x.k+y.k)%mod,(x.b+y.b)%mod);}
void operator+=(const Poly&y){*this=*this+y;}
friend Poly operator-(const Poly&x,const Poly&y){return Poly((x.k-y.k+mod)%mod,(x.b-y.b+mod)%mod);}
void operator-=(const Poly&y){*this=*this-y;}
friend Poly operator*(const Poly&x,const Poly&y){return Poly((1ll*x.k*y.b+1ll*x.b*y.k)%mod,1ll*x.b*y.b%mod);}
void operator*=(const Poly&y){*this=*this*y;}
friend Poly operator/(const Poly&x,const Poly&y){return Poly((1ll*x.k*y.b%mod-1ll*x.b*y.k%mod+mod)%mod*ksm(1ll*y.b*y.b%mod)%mod,1ll*x.b*ksm(y.b)%mod);}
void operator/=(const Poly&y){*this=*this/y;}
bool operator~()const{return b!=0;}
Poly operator-()const{return Poly((mod-k)%mod,(mod-b)%mod);}
void print(){printf("(%dx+%d)",k,b);}
}a[32][32];
void Eular(){
phi[1]=1;
for(int i=2;i<=lim;i++){
if(!pri[i])pri[++pri[0]]=i,phi[i]=i-1;
for(int j=1;j<=pri[0]&&i*pri[j]<=lim;j++){
pri[i*pri[j]]=true;
if(!(i%pri[j])){phi[i*pri[j]]=phi[i]*pri[j];break;}
else phi[i*pri[j]]=phi[i]*phi[pri[j]];
}
}
}
int Determinant(){
Poly ret(0,1);
for(int i=1;i<n;i++){
int j=i;
for(;j<n;j++)if(~a[j][i])break;
if(j>=n)return 0;
if(j!=i)swap(a[j],a[i]),ret=-ret;
Poly inv=Poly(0,1)/a[i][i];
for(j=i+1;j<n;j++){
Poly lam=a[j][i]*inv;
for(int k=i;k<n;k++)a[j][k]-=lam*a[i][k];
}
ret*=a[i][i];
}
for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)a[i][j]=Poly();
return ret.k;
}
int main(){
scanf("%d%d",&n,&m);
for(int i=1;i<=m;i++)scanf("%d%d%d",&u[i],&v[i],&w[i]),lim=max(lim,w[i]);
Eular();
for(int i=1;i<=lim;i++){
for(int j=1;j<=n;j++)dsu[j]=j;
for(int j=1;j<=m;j++)if(!(w[j]%i))merge(u[j],v[j]);
int cnt=0;
for(int j=1;j<=n;j++)cnt+=(dsu[j]==j);
if(cnt!=1)continue;
for(int j=1;j<=m;j++)if(!(w[j]%i))a[u[j]][u[j]]+=Poly(w[j]),a[v[j]][v[j]]+=Poly(w[j]),a[u[j]][v[j]]-=Poly(w[j]),a[v[j]][u[j]]-=Poly(w[j]);
(res+=1ll*phi[i]*Determinant()%mod)%=mod;
}
printf("%d\n",res);
return 0;
}
```
# IV.BEST定理
BEST定理是用于有向图欧拉回路计数的一个定理。内容如下：
$$\boxed{ans=T(x)\prod\limits_{v\in\mathbb V}(deg_v-1)!}$$
其中 $T(x)$ 是以 $x$ 为根的内向树数，$deg_v$ 为节点 $v$ 的度数（显然，一张有欧拉回路的图应保证入度等于出度，且图联通）。
（实际上，因为保证出入度相等，因此这张图的内外向树数应该相等）
若其指定了起点为 $s$，则方案数就变成了
$$\boxed{ans=T(s)deg_s\prod\limits_{v\in\mathbb V}(deg_v-1)!}$$
因为其可以被看作是整个大回路在 $s$ 的某条出边处被截开了。
有亿些细节需要注意：
1. 若存在自环，就把它同普通边一样处理即可。
2. 孤立点要舍去。
3. 记得判连通性。（可以用冰茶姬）
4. 记得特判无边的情形。
5. 记得特判出入边是否相等。
6. 在有起点时别忘了乘上 $deg_s$。
7. 若重边间不计算顺序记得除去重边的阶乘。
时间复杂度 $O(n^3)$。
## I.[Which Dreamed It /【模板】BEST 定理](https://www.luogu.com.cn/problem/P5807)
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
const int mod=1e6+3;
const int N=2e5;
int T,n,g[110][110],h[110][110],deg[110],ged[110],m,id[110],fac[200100],caf[200100],res,dsu[110];
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
int Determinant(){
int ret=1;
for(int i=2;i<=m;i++){
int j=i;
while(j<=m&&!h[j][i])j++;
if(j>m)return 0;
if(j!=i)swap(h[j],h[i]),ret=-ret;
int inv=mod-ksm(h[i][i]);
for(j=i+1;j<=m;j++){
int lam=1ll*h[j][i]*inv%mod;
for(int k=i;k<=m;k++)(h[j][k]+=1ll*lam*h[i][k]%mod)%=mod;
}
}
(ret+=mod)%=mod;
for(int i=2;i<=m;i++)ret=1ll*ret*h[i][i]%mod;
return ret;
}
int find(int x){return dsu[x]==x?x:dsu[x]=find(dsu[x]);}
void merge(int x,int y){x=find(x),y=find(y);if(x!=y)dsu[y]=x;}
int main(){
fac[0]=1;for(int i=1;i<=N;i++)fac[i]=1ll*fac[i-1]*i%mod;
caf[n]=ksm(fac[n]);for(int i=N;i;i--)caf[i-1]=1ll*caf[i]*i%mod;
scanf("%d",&T);
while(T--){
scanf("%d",&n),m=0,memset(g,0,sizeof(g)),memset(h,0,sizeof(h)),memset(id,0,sizeof(id)),memset(deg,0,sizeof(deg)),memset(ged,0,sizeof(ged)),res=1;
for(int i=1;i<=n;i++)dsu[i]=i;
for(int i=1,x,y;i<=n;i++){
scanf("%d",&x),ged[i]=x;
while(x--)scanf("%d",&y),g[i][y]++,deg[y]++,merge(i,y);
}
// for(int i=1;i<=n;i++){for(int j=1;j<=n;j++)printf("%d ",g[i][j]);puts("");}
bool ok=true;
for(int i=1;i<=n;i++)if(deg[i]!=ged[i])ok=false;
for(int i=2;i<=n;i++)if(deg[i]&&find(i)!=find(1))ok=false;
if(!ok){puts("0");continue;}
if(deg[1]==g[1][1]){printf("%d\n",fac[deg[1]]);continue;}
for(int i=1;i<=n;i++)if(deg[i])id[i]=++m;
for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)if(id[i]&&id[j])h[id[i]][id[j]]=(mod-g[i][j])%mod;
for(int i=1;i<=m;i++)if(deg[i])(h[id[i]][id[i]]+=deg[i])%=mod;
// for(int i=1;i<=n;i++){for(int j=1;j<=n;j++)printf("%d ",h[i][j]);puts("");}
for(int i=1;i<=n;i++)if(deg[i])res=1ll*res*fac[deg[i]-1]%mod;
res=1ll*res*Determinant()%mod*deg[1]%mod;
printf("%d\n",res);
}
return 0;
}
```
# V.虚树
虚树，是一种针对**树上点集**的强力运算。它可以在$O(k\log k)$（其中$k$是点集大小）的时间内，建出一棵包含点集中所有节点，以及其中某些点的lca的树出来。这棵树就被称作**虚树**。之后就可以在虚树上进行操作了，例如树形DP等。
建出虚树的操作主要是这样的：
1. 我们维护一个栈，从栈顶到栈底构成原树上一条**深度递减**的链。注意这条链不是完整的链——它只记录了链上的关键节点。
2. 我们首先将原树的根节点加入虚树，方便维护。接着，将点集中所有点**按照dfs序排序**，然后依次加入虚树。
3. 当一个点被加入虚树时，如果栈为空，直接将该点加入栈，并直接结束函数；否则，求出它与栈顶节点的LCA。
4. 如果栈中第二个元素的深度大于等于LCA的深度，在虚树上连一条边$(\text{栈中第二个元素},\text{栈顶元素})$，并弹出栈顶节点。重复此操作，直到第二个元素的深度小于LCA的深度。
5. 如果栈顶元素的深度小于LCA的深度，连一条边$(LCA,\text{栈顶元素})$，并弹出栈顶元素。
6. 如果栈顶元素不同于LCA，将LCA入栈。
7. 将（3）中那个点入栈，并结束该函数。
8. 在所有东西全部加入虚树后，清空栈，在清空的过程中同时连边。
代码：
```cpp
void ins(int x){
if(!tp){stk[++tp]=x;return;}
int lca=LCA(x,stk[tp]);
while(tp>=2&&dep[lca]<dep[stk[tp-1]])v[stk[tp-1]].push_back(stk[tp]),tp--;
if(tp&&dep[lca]<dep[stk[tp]])v[lca].push_back(stk[tp--]);
if(!tp||stk[tp]!=lca)stk[++tp]=lca;
stk[++tp]=x;
}
```
理解：
因为我们事先将节点按照dfs序排序，所以当新节点的深度大于栈顶节点深度的时候，就意味着**原树中栈顶节点子树内所有节点都已经被加入虚树**，可以弹出了。
然后，这个lca为了维护虚树的形态，也应该被加入虚树。
总复杂度$O(k\log k)$，瓶颈在于按照dfs序排序。
------------
一些定义：
虚树：即上文代码所建出的树
原树/实树：原本的树
实点：给出点集中的点
虚点：在建出虚树时为了维护形态而添加的点（即实点的LCA之类）
树外点：原树中不在虚树内的点
树点：虚树内的点（大部分没有歧义的时候，“节点”就特指树点）
## I.[CF613D Kingdom and its Cities](https://www.luogu.com.cn/problem/CF613D)
建虚树时有几个事是一定不能忘的：
1. 记得将节点按照dfs序排序；
2. 记得将根节点加入虚树；
3. 该清空的一定都得清空。
这题我们就可以使用虚树解决。
首先，我们建出虚树，并给所有节点一个$sz$，其中只有点集中的节点的$sz$是$1$（因为虚树中还有非点集的点，即它们的LCA）。
我们可以特判掉无解的情况，即点集中有一对点**直接相连**。其它情况都有解。
我们考虑对于一个有$sz$的点，如果它的儿子也有$sz$，则原树中（从这个点到儿子的路径中）至少要删一个点，令答案加一；
否则，对于一个无$sz$的点，如果它有$sz$的儿子数量大于$1$，则删去这个点，令答案加一；否则，即它只有一个有$sz$的儿子或是一个都没有，直接继承那个儿子的$sz$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
int n,m,dfn[200100][20],DFN[100100],TOT,LG[200100],dep[100100],in[100100],tot,sz[100100];
namespace real{
vector<int>v[100100];
void dfs(int x,int fa){
dfn[in[x]=++tot][0]=x,dep[x]=dep[fa]+1,DFN[x]=++TOT;
for(auto y:v[x])if(y!=fa)dfs(y,x),dfn[++tot][0]=x;
}
}
int MIN(int x,int y){
return dep[x]<dep[y]?x:y;
}
int LCA(int x,int y){
x=in[x],y=in[y];
if(x>y)swap(x,y);
int k=LG[y-x+1];
return MIN(dfn[x][k],dfn[y-(1<<k)+1][k]);
}
namespace imag{
int stk[100100],tp,sz[100100],res,a[100100];
vector<int>v[100100];
bool cmp(int x,int y){
return DFN[x]<DFN[y];
}
void ins(int x){
sz[x]=1;
if(!tp){stk[++tp]=x;return;}
int lca=LCA(x,stk[tp]);
while(tp>=2&&dep[lca]<dep[stk[tp-1]])v[stk[tp-1]].push_back(stk[tp]),tp--;
if(tp&&dep[lca]<dep[stk[tp]])v[lca].push_back(stk[tp--]);
if(!tp||stk[tp]!=lca)stk[++tp]=lca;
stk[++tp]=x;
}
void fin(){
while(tp>=2)v[stk[tp-1]].push_back(stk[tp]),tp--;
tp--;
}
void dfs1(int x){
for(auto y:v[x]){
if(sz[x]&&sz[y]&&dep[y]==dep[x]+1)res=-1;
dfs1(y);
}
}
void dfs2(int x){
int tmp=sz[x];
for(auto y:v[x]){
dfs2(y);
if(tmp&&sz[y])res++;
if(!tmp)sz[x]+=sz[y];
}
if(!tmp&&sz[x]>1)sz[x]=0,res++;
}
void dfs3(int x){
sz[x]=0;
for(auto y:v[x])dfs3(y);
v[x].clear();
}
void work(){
int q;
res=0,scanf("%d",&q);
for(int i=1;i<=q;i++)scanf("%d",&a[i]);
sort(a+1,a+q+1,cmp);
if(a[1]!=1)stk[++tp]=1;
for(int i=1;i<=q;i++)ins(a[i]);
fin();
dfs1(1);
if(res==-1)puts("-1");
else dfs2(1),printf("%d\n",res);
dfs3(1);
}
}
int main(){
scanf("%d",&n);
for(int i=1,x,y;i<n;i++)scanf("%d%d",&x,&y),real::v[x].push_back(y),real::v[y].push_back(x);
real::dfs(1,0);
for(int i=2;i<=tot;i++)LG[i]=LG[i>>1]+1;
for(int j=1;j<=LG[tot];j++)for(int i=1;i+(1<<j)<=tot;i++)dfn[i][j]=MIN(dfn[i][j-1],dfn[i+(1<<(j-1))][j-1]);
scanf("%d",&m);
while(m--)imag::work();
return 0;
}
```
## II.[[SDOI2011]消耗战](https://www.luogu.com.cn/problem/P2495)
老套路，我们建出虚树。
这题虚树中的边是**带边权**的，边权为原树中两点路径中权值的$\min$。这个权值的$\min$可以通过倍增求出。
在建出虚树后，我们就可以考虑DP了。设当前点为$x$，$f_x$为$x$同$x$子树中所有“资源丰富”节点切断的最小代价。
如果$x$本身是“资源丰富”的点，那直接```return```，不需要访问子树；
否则，我们需要访问子树。设$y$为$x$的一个儿子，$(x,y)$的边权为$z$。
如果$y$是资源节点，$f_x$应该增加$z$；
否则，因为$y$节点可以在$(x,y)$边处或者$y$子树内切断，$f_x$应该增加$\min(z,f_y)$。
则最终答案即为$f_1$。
复杂度$O(n\log n)$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int n,m,dfn[250010],dep[250010],tot,anc[250010][20],mn[250010][20];
namespace real{
int head[250010],cnt;
struct node{
int to,next,val;
}edge[500010];
void ae(int u,int v,int w){
edge[cnt].next=head[u],edge[cnt].to=v,edge[cnt].val=w,head[u]=cnt++;
edge[cnt].next=head[v],edge[cnt].to=u,edge[cnt].val=w,head[v]=cnt++;
}
void dfs(int x){
dfn[x]=++tot;
for(int i=head[x],y;i!=-1;i=edge[i].next)if((y=edge[i].to)!=anc[x][0])dep[y]=dep[x]+1,anc[y][0]=x,mn[y][0]=edge[i].val,dfs(y);
}
}
int LCA(int x,int y){
if(dep[x]<dep[y])swap(x,y);
for(int i=19;i>=0;i--)if(dep[x]-(1<<i)>=dep[y])x=anc[x][i];
if(x==y)return x;
for(int i=19;i>=0;i--)if(anc[x][i]!=anc[y][i])x=anc[x][i],y=anc[y][i];
return anc[x][0];
}
int route(int x,int y){
int res=0x3f3f3f3f;
for(int i=19;i>=0;i--)if(dep[x]-(1<<i)>=dep[y])res=min(res,mn[x][i]),x=anc[x][i];
return res;
}
namespace imag{
int head[250010],cnt,stk[250010],tp,a[250010];
ll f[250010];
bool mark[250010];
struct node{
int to,next,val;
}edge[250010];
void ae(int u,int v){
edge[cnt].next=head[u],edge[cnt].to=v,edge[cnt].val=route(v,u),head[u]=cnt++;
}
void ins(int x){
mark[x]=true;
if(!tp){stk[++tp]=x;return;}
int lca=LCA(x,stk[tp]);
while(tp>=2&&dep[lca]<dep[stk[tp-1]])ae(stk[tp-1],stk[tp]),tp--;
if(tp&&dep[lca]<dep[stk[tp]])ae(lca,stk[tp--]);
if(!tp||stk[tp]!=lca)stk[++tp]=lca;
stk[++tp]=x;
}
void fin(){
while(tp>=2)ae(stk[tp-1],stk[tp]),tp--;
tp--;
}
bool cmp(int x,int y){
return dfn[x]<dfn[y];
}
void dfs1(int x){
if(mark[x])return;
for(int i=head[x];i!=-1;i=edge[i].next){
dfs1(edge[i].to);
if(mark[edge[i].to])f[x]+=edge[i].val;
else f[x]+=min(1ll*edge[i].val,f[edge[i].to]);
}
}
void dfs2(int x){
mark[x]=false,f[x]=0;
for(int i=head[x];i!=-1;i=edge[i].next)dfs2(edge[i].to);
head[x]=-1;
}
void work(){
int q;
scanf("%d",&q);
for(int i=1;i<=q;i++)scanf("%d",&a[i]);
sort(a+1,a+q+1,cmp);
stk[++tp]=1;
for(int i=1;i<=q;i++)ins(a[i]);
fin();
dfs1(1);
printf("%lld\n",f[1]);
dfs2(1),cnt=0;
}
}
int main(){
scanf("%d",&n),memset(real::head,-1,sizeof(real::head)),memset(imag::head,-1,sizeof(imag::head));
for(int i=1,x,y,z;i<n;i++)scanf("%d%d%d",&x,&y,&z),real::ae(x,y,z);
real::dfs(1);
for(int j=1;j<20;j++)for(int i=1;i<=n;i++)anc[i][j]=anc[anc[i][j-1]][j-1],mn[i][j]=min(mn[i][j-1],mn[anc[i][j-1]][j-1]);
scanf("%d",&m);
while(m--)imag::work();
return 0;
}
```
## III.[[HEOI2014]大工程](https://www.luogu.com.cn/problem/P4103)
仍然建出虚树。
我们考虑设$sz_x$表示$x$子树中实点（即原本点集中的点）的数量，再设$f_x$表示$x$到$x$子树中某个实点的最长路径，$g_x$则表示最短路径。
我们先考虑求$\min$和$\max$的部分。
对于一个实点，它初始值$f_x=g_x=0$（有一条到自己的路径）；对于一条虚点，它初始$f_x=-\infty,g_x=\infty$。
我们枚举它的儿子$y$。设$z$表示这条边的边权（实际上就是$y$的深度减去$x$的深度），则我们可以用$f_x+z+f_y$拼成一条最长路径，$g_x+z+g_y$拼成一条最短路径。然后更新$f$与$g$即可，就跟正常树中求直径的做法没啥区别。
我们再考虑求和的地方。
显然，对于一条边$(x,y)$（这里令$x$为$y$在虚树上的父亲），它总共会被经过$sz_y\times(q-sz_y)$次，其中$q$是点集大小。这也很好理解——对于$y$侧树中每个实点，它都会与$x$侧树中每个实点有一条边。则只需要用次数乘上边权，再求和，便是答案。
这题我担心$n\log n$空间的倍增或者ST表可能会MLE，就使用了树剖。但是总复杂度仍然是$O(n\log n)$的。
代码：
```
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int n,m,fa[1001000],dep[1001000],dfn[1001000],top[1001000];
namespace real{
int sz[1001000],son[1001000],tot;
vector<int>v[1001000];
void dfs1(int x,int Fa){
fa[x]=Fa,dep[x]=dep[Fa]+1,dfn[x]=++tot,sz[x]=1;
for(auto y:v[x]){
if(y==fa[x])continue;
dfs1(y,x),sz[x]+=sz[y];
if(sz[son[x]]<sz[y])son[x]=y;
}
}
void dfs2(int x){
if(son[x])top[son[x]]=top[x],dfs2(son[x]);
for(auto y:v[x])if(y!=fa[x]&&y!=son[x])top[y]=y,dfs2(y);
}
}
int LCA(int x,int y){
while(top[x]!=top[y])dep[top[x]]>dep[top[y]]?x=fa[top[x]]:y=fa[top[y]];
return dep[x]<dep[y]?x:y;
}
namespace imag{
int q,stk[1001000],tp,a[1001000],sz[1001000],mn,mx,f[1001000],g[1001000];
ll res;
vector<int>v[1001000];
bool cmp(int x,int y){
return dfn[x]<dfn[y];
}
void ins(int x){
sz[x]=1;
if(!tp){stk[++tp]=x;return;}
int lca=LCA(x,stk[tp]);
while(tp>=2&&dep[lca]<dep[stk[tp-1]])v[stk[tp-1]].push_back(stk[tp]),tp--;
if(tp&&dep[lca]<dep[stk[tp]])v[lca].push_back(stk[tp--]);
if(!tp||stk[tp]!=lca)stk[++tp]=lca;
stk[++tp]=x;
}
void fin(){
while(tp>=2)v[stk[tp-1]].push_back(stk[tp]),tp--;
tp--;
}
void dfs0(int x){
if(sz[x])f[x]=g[x]=0;
else f[x]=-0x3f3f3f3f,g[x]=0x3f3f3f3f;
for(auto y:v[x]){
dfs0(y);
mx=max(mx,f[x]+(dep[y]-dep[x])+f[y]),f[x]=max(f[x],(dep[y]-dep[x])+f[y]);
mn=min(mn,g[x]+(dep[y]-dep[x])+g[y]),g[x]=min(g[x],(dep[y]-dep[x])+g[y]);
}
}
void dfs1(int x){
for(auto y:v[x])dfs1(y),res+=1ll*(dep[y]-dep[x])*sz[y]*(q-sz[y]),sz[x]+=sz[y];
}
void dfs2(int x){
sz[x]=0;
for(auto y:v[x])dfs2(y);
v[x].clear();
}
void work(){
res=0,mn=0x3f3f3f3f,mx=-0x3f3f3f3f,scanf("%d",&q);
for(int i=1;i<=q;i++)scanf("%d",&a[i]);
sort(a+1,a+q+1,cmp);
if(a[1]!=1)stk[++tp]=1;
for(int i=1;i<=q;i++)ins(a[i]);
fin();
dfs0(1);
dfs1(1);
printf("%lld %lld %lld\n",res,mn,mx);
dfs2(1);
}
}
int main(){
scanf("%d",&n);
for(int i=1,x,y;i<n;i++)scanf("%d%d",&x,&y),real::v[x].push_back(y),real::v[y].push_back(x);
real::dfs1(1,0),top[1]=1,real::dfs2(1);
scanf("%d",&m);
while(m--)imag::work();
return 0;
}
```
## IV.[[HNOI2014]世界树](https://www.luogu.com.cn/problem/P3233)
~~人傻常数大没错了，$n\log n$还会TLE~~
首先当然是建出虚树来。
然后，对于虚树中每个节点（不管是否是实点），我们可以DP出管辖它的那个节点，设为$f_x$。这个可以通过**二次扫描与换根法**在$O(k)$的时间内通过两次dfs求出，假如你使用ST表求LCA的话。这里我使用了倍增，因此复杂度为$O(k\log n)$。
我们设$g_x$表示$x$节点（必须是个实点）管辖了多少节点，即询问的答案。
首先，我们先求出**不在虚树中的子树**的贡献。设$sz_x$为原树中$x$节点的子树大小。则对于虚树中每个点$x$，$sz_x-\sum\limits_{y\in son_x}sz_y$即是不在虚树中子树的贡献（$x$自身也包含在内）。这个可以通过一次dfs求出。
但是，这个$sz_y$可不能直接用$y$的$sz$——它应该是**原树上$x$在$y$方向的那个儿子（设为$Q$）的**$sz$。这个可以通过倍增求出，当然如果你不嫌麻烦，可以通过**长链剖分**来在$O(1)$时间内求出**树上$k$级祖先**。因为前面LCA我就使用了倍增，这里仍然使用倍增。
然后就是**虚树上路径的贡献**了。对于虚树上一条边$(x,y)$且$x$是$y$的父亲，如果$f_x=f_y$，则整条边都归$f_x$管辖，贡献是$sz_Q-sz_y$，其中$Q$是$x$在$y$方向的儿子。
否则，即$f_x\neq f_y$，我们可以直接找出从$y$点往上有$P$个节点是归$f_y$管辖的。则一共有$sz_P-sz_y$个节点归$y$管理，$sz_Q-sz_P$个节点归$x$管理。
注意到这里的$sz_Q$在子树贡献和路径贡献中被算了两次，一次加一次减。故两次可以直接抵消。
复杂度$O(n\log n)$（就算你丧心病狂用了ST表和长链剖分，复杂度瓶颈的```sort```你仍无法解决）
代码：
```cpp
#pragma GCC optimize(3)
#include<bits/stdc++.h>
using namespace std;
int n,m,anc[300100][20],dfn[300100],tot,dep[300100],sz[300100];
namespace real{
vector<int>v[300100];
void dfs(int x,int fa){
anc[x][0]=fa,dep[x]=dep[fa]+1,dfn[x]=++tot,sz[x]=1;
for(auto y:v[x])if(y!=fa)dfs(y,x),sz[x]+=sz[y];
}
}
int MIN(int x,int y){
return dep[x]<dep[y]?x:y;
}
int LCA(int x,int y){
if(dep[x]<dep[y])swap(x,y);
for(int i=19;i>=0;i--)if(dep[x]-(1<<i)>=dep[y])x=anc[x][i];
if(x==y)return x;
for(int i=19;i>=0;i--)if(anc[x][i]!=anc[y][i])x=anc[x][i],y=anc[y][i];
return anc[x][0];
}
int DIS(int x,int y){
return dep[x]+dep[y]-2*dep[LCA(x,y)];
}
int ANC(int x,int z){
for(int i=19;i>=0;i--)if((1<<i)<=z)x=anc[x][i],z-=(1<<i);
return x;
}
namespace imag{
int stk[300100],tp,res,a[300100],r[300100],f[300100],g[300100];
bool mark[300100];
vector<int>v[300100];
bool cmp(int x,int y){
return dfn[x]<dfn[y];
}
void ins(int x){
mark[x]=true;
if(!tp){stk[++tp]=x;return;}
int lca=LCA(x,stk[tp]);
while(tp>=2&&dep[lca]<dep[stk[tp-1]])v[stk[tp-1]].push_back(stk[tp]),tp--;
if(tp&&dep[lca]<dep[stk[tp]])v[lca].push_back(stk[tp--]);
if(!tp||stk[tp]!=lca)stk[++tp]=lca;
stk[++tp]=x;
}
void fin(){
while(tp>=2)v[stk[tp-1]].push_back(stk[tp]),tp--;
tp--;
}
void dfs1(int x){
if(mark[x])f[x]=x;
else f[x]=0;
for(auto y:v[x]){
dfs1(y);
if(!f[y])continue;
if(!f[x])f[x]=f[y];
else{
int X=DIS(x,f[x]),Y=DIS(x,f[y]);
if(X>Y)f[x]=f[y];
if(X==Y&&f[x]>f[y])f[x]=f[y];
}
}
}
void dfs2(int x){
for(auto y:v[x]){
if(!f[x])continue;
if(!f[y])f[y]=f[x];
else{
int X=DIS(y,f[x]),Y=DIS(y,f[y]);
if(X<Y)f[y]=f[x];
if(X==Y&&f[x]<f[y])f[y]=f[x];
}
dfs2(y);
}
}
void dfs3(int x){
g[f[x]]+=sz[x];
for(auto y:v[x])dfs3(y);
}
void dfs4(int x){
for(auto y:v[x]){
dfs4(y);
if(f[x]==f[y])g[f[x]]-=sz[y];
else{
int X=DIS(x,f[x]),Y=DIS(y,f[y]),Z=dep[y]-dep[x]-1,P=0;
int D=min(abs(X-Y),Z);
if(Y<X)P+=D;
Z-=D;
P+=(Z+(f[y]<f[x]))/2;
P=ANC(y,P);
// printf("(%d,%d):%d\n",x,y,P);
g[f[y]]+=sz[P]-sz[y];
g[f[x]]-=sz[P];
}
}
}
void dfs5(int x){
g[x]=mark[x]=0;
for(auto y:v[x])dfs5(y);
v[x].clear();
}
void work(){
int q;
res=0,scanf("%d",&q);
for(int i=1;i<=q;i++)scanf("%d",&a[i]),r[i]=a[i];
sort(a+1,a+q+1,cmp);
if(a[1]!=1)stk[++tp]=1;
for(int i=1;i<=q;i++)ins(a[i]);
fin();
dfs1(1);
dfs2(1);
dfs3(1);
// for(int i=1;i<=q;i++)printf("%d ",g[r[i]]);puts("");
dfs4(1);
for(int i=1;i<=q;i++)printf("%d ",g[r[i]]);puts("");
dfs5(1);
}
}
int main(){
scanf("%d",&n);
for(int i=1,x,y;i<n;i++)scanf("%d%d",&x,&y),real::v[x].push_back(y),real::v[y].push_back(x);
real::dfs(1,0);
for(int j=1;j<20;j++)for(int i=1;i<=n;i++)anc[i][j]=anc[anc[i][j-1]][j-1];
scanf("%d",&m);
while(m--)imag::work();
return 0;
}
```
## V.[CF639F Bear and Chemistry](https://www.luogu.com.cn/problem/CF639F)
大毒瘤题一道。
先边双缩点缩成森林，再对每组询问建出虚树，再连上边跑Tarjan求边双即可。
~~是不是很simple？但是相信我，码起来会让你发疯的~~
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
int n,m,q,id[300100],lim,tree[300100],totans;
int read(){
int x=0;
char c=getchar();
while(c>'9'||c<'0')c=getchar();
while(c>='0'&&c<='9')x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x;
}
int deco(){
int x=(read()+totans)%n;
if(!x)x=n;
// printf("DECO:%d\n",x);
return id[x];
}
namespace ECC{
int dfn[300100],low[300100],tot,stk[300100],tp,head[300100],cnt,col[300100],c;
struct node{
int to,next;
}edge[1200100];
void ae(int u,int v){
// printf("%d %d\n",u,v);
edge[cnt].next=head[u],edge[cnt].to=v,head[u]=cnt++;
edge[cnt].next=head[v],edge[cnt].to=u,head[v]=cnt++;
}
void Tarjan(int x,int fa){
low[x]=dfn[x]=++tot,stk[++tp]=x;
for(int i=head[x],y;i!=-1;i=edge[i].next){
if(!dfn[y=edge[i].to])Tarjan(y,i),low[x]=min(low[x],low[y]);
else if((i^fa)!=1)low[x]=min(low[x],dfn[y]);
}
if(low[x]!=dfn[x])return;
c++;
int y;
do col[y=stk[tp--]]=c;while(y!=x);
}
}
namespace Tree{
vector<int>v[300100];
int dfn[300100],tot,anc[300100][20],dep[300100];
int LCA(int x,int y){
if(dep[x]>dep[y])swap(x,y);
for(int i=19;i>=0;i--)if(dep[x]<=dep[y]-(1<<i))y=anc[y][i];
if(x==y)return x;
for(int i=19;i>=0;i--)if(anc[x][i]!=anc[y][i])x=anc[x][i],y=anc[y][i];
return anc[x][0];
}
void dfs(int x,int fa,int ID){
tree[x]=ID,dfn[x]=++tot,dep[x]=dep[fa]+1,anc[x][0]=fa;
for(int i=1;anc[x][i-1];i++)anc[x][i]=anc[anc[x][i-1]][i-1];
for(auto y:v[x])if(y!=fa)dfs(y,x,ID);
}
namespace imag{
int stk[300100],tp,num[300100],len,ed,inq[300100],all,U[300100],V[300100];
vector<int>sp;
void ins(int x){
sp.push_back(x);
if(!tp){stk[++tp]=x;return;}
int lca=LCA(x,stk[tp]);
while(tp>=2&&dep[lca]<dep[stk[tp-1]])ECC::ae(stk[tp-1],stk[tp]),tp--;
if(tp&&dep[lca]<dep[stk[tp]])ECC::ae(lca,stk[tp--]);
if(!tp||stk[tp]!=lca)stk[++tp]=lca,sp.push_back(lca);
stk[++tp]=x;
}
void fin(){
while(tp>=2)ECC::ae(stk[tp-1],stk[tp]),tp--;
tp=0;
}
bool solve(int ip){
all=len=read(),ed=read();
for(int i=1;i<=len;i++)inq[i]=num[i]=deco();
for(int i=1;i<=ed;i++)U[i]=deco(),V[i]=deco(),num[++all]=U[i],num[++all]=V[i];
sort(num+1,num+all+1,[](int x,int y){return tree[x]==tree[y]?dfn[x]<dfn[y]:tree[x]<tree[y];});
all=unique(num+1,num+all+1)-num-1;
// printf("QWQ:");for(int i=1;i<=len;i++)printf("%d ",inq[i]);puts("");
// for(int i=1;i<=ed;i++)printf("(%d,%d)",U[i],V[i]);puts("");
// printf("%d::",all);for(int i=1;i<=all;i++)printf("%d ",num[i]);puts("");
ECC::tot=ECC::cnt=ECC::c=0;
for(int i=1;i<=ed;i++)ECC::ae(U[i],V[i]);
for(int i=1;i<=all;i++){
if(tree[num[i]]!=tree[num[i-1]]){
fin();
if(num[i]!=tree[num[i]])sp.push_back(tree[num[i]]),stk[++tp]=tree[num[i]];
}
ins(num[i]);
}
fin();
for(int i=1;i<=all;i++)if(!ECC::dfn[num[i]])ECC::Tarjan(num[i],-1);
bool ok=true;
for(int i=2;i<=len;i++)if(ECC::col[inq[i]]!=ECC::col[inq[i-1]])ok=false;
for(auto i:sp)ECC::head[i]=-1,ECC::col[i]=ECC::dfn[i]=ECC::low[i]=0;
sp.clear();
if(!ok)return false;
(totans+=ip)%=n;
return true;
}
}
}
int main(){
n=read(),m=read(),q=read(),memset(ECC::head,-1,sizeof(ECC::head));
for(int i=1,x,y;i<=m;i++)x=read(),y=read(),ECC::ae(x,y);
for(int i=1;i<=n;i++)if(!ECC::dfn[i])ECC::Tarjan(i,-1);
lim=ECC::c;for(int i=1;i<=n;i++)id[i]=ECC::col[i];
for(int i=1;i<=n;i++)for(int j=ECC::head[i];j!=-1;j=ECC::edge[j].next)if(id[i]!=id[ECC::edge[j].to])Tree::v[id[i]].push_back(id[ECC::edge[j].to]);
for(int i=1;i<=lim;i++)ECC::head[i]=-1,ECC::col[i]=ECC::dfn[i]=ECC::low[i]=0;
for(int i=1;i<=lim;i++)if(!tree[i])Tree::dfs(i,0,i);
// for(int i=1;i<=n;i++)printf("%d ",id[i]);puts("");
// for(int i=1;i<=lim;i++)printf("%d ",tree[i]);puts("");
for(int i=1;i<=q;i++)puts(Tree::imag::solve(i)?"YES":"NO");
return 0;
}
```
## VI.[LOJ#3077. 「2019 集训队互测 Day 4」绝目编诗](https://loj.ac/p/3077)
神题。
乍一看好像和虚树半毛钱关系都没有呀？没关系，过亿会就有了。
我们不妨先从暴力开始想起。
~~暴力怎么写？暴力怎么写？加边加边加边，搜就完事了。~~
没错，这里的暴力就是爆搜——搜出所有环来，然后判断是否有两个环长度相等即可。
但是这里的暴力也需要一些艺术——因为我们要不重不漏地统计所有的环。我们得加上这几条限制：
1. 爆搜时，只搜**从一个点出发再回到该点本身**的简单环，不包括 $\rho$ 型甚至更诡异形状的。
2. 既然是简单环，我们就不允许经过同一个点两次。
3. 为了避免将长度为 $2$ 的简单环也计入，我们还要记录环上前一个节点，避免走反边回到上一个节点。
这样，我们就保证可以搜出所有简单环。但是，一个简单环可能被重复搜出。更具体地说，对于一个长度为 $x$ 的简单环，其会旋转 $x$ 次，翻转 $1$ 次，总计被搜出 $2x$ 次。因此，当我们发现某一个长度的环出现了至少 $2x+1$ 次，才可以判出现了长度相同的环。
暴力的复杂度是指数级别的，期望得分 $60\%$。
代码：
```cpp
//Version: Completely brute solution
#include<bits/stdc++.h>
using namespace std;
int n,m,sz[10100],dep[10100];
vector<int>v[10100];
void dfs(int x,int fa){
dep[x]=dep[fa]+1;
for(auto y:v[x])if(y!=fa){
if(dep[y]){
if(dep[y]!=1)continue;
sz[dep[x]]++;
if(sz[dep[x]]>2*dep[x]){puts("Yes");exit(0);}
continue;
}
dfs(y,x);
}
dep[x]=0;
}
int main(){
scanf("%d%d",&n,&m);
if(m>=2*n){puts("Yes");return 0;}
for(int i=1,x,y;i<=m;i++)scanf("%d%d",&x,&y),v[x].push_back(y),v[y].push_back(x);
for(int i=1;i<=n;i++)dfs(i,0);
puts("No");
return 0;
}
```
考虑优化。
首先，众所周知，一个 $n$ 个点的简单图，其最大简单环的长度最多只能为 $n$。又因为不存在长度为 $1,2$ 的简单环，所以若一个简单图不存在长度相等的简单环，其所拥有的简单环数量最多是 $3\sim n$ 各一个，共 $n-2$ 个。
考虑求出该图的任意一极大生成森林，其边数最多是 $n-1$ 条。考虑往上加任意一条边，环数都必定增加至少 $1$。故添加至多 $n-2$ 条边后，$3\sim n$ 所有长度的环都出现了。再加任意一条边，都会出现长度相同的环。因此，不存在长度相同环的图，其 $m$ 最大只能是 $(n-1)+(n-2)=2n-3$，超过此界直接判 `Yes` 即可。
于是 $m$ 就被压缩到了和 $n$ 同级。
下面我们加入第二个优化。因为 `No` 的最多情况下，环的总长是 $3\sim n$ 各出现一次，共 $n^2$ 级别。因此我们考虑仅访问环边，复杂度就是 $n^3$ 的（因为每个环都会被访问 $2x$ 次）。如何做到呢？很简单，只需在每次到达一个新节点时，$O(n+m)$ 地 `dfs` 一下看能不能回到起点，不能就直接剪枝剪掉即可。这样子，复杂度就被优化到了严格 $O(n^4)$，再也不是什么玄学的指数复杂度了。
期望得分 $80\%$。
代码：
```cpp
//Version: improved brute solution (n^4)
#include<bits/stdc++.h>
using namespace std;
int n,m,sz[10100],dep[10100],vis[10100],id;
vector<int>v[10100];
bool reach(int x,int fa,int O){
vis[x]=O;
for(auto y:v[x])if(y!=fa){
if(dep[y]==1)return true;
if(dep[y]||vis[y]==O)continue;
if(reach(y,x,O))return true;
}
return false;
}
void dfs(int x,int fa){
dep[x]=dep[fa]+1;
if(!reach(x,fa,++id)){dep[x]=0;return;}
for(auto y:v[x])if(y!=fa){
if(dep[y]){
if(dep[y]!=1)continue;
sz[dep[x]]++;
if(sz[dep[x]]>2*dep[x]){puts("Yes");exit(0);}
continue;
}
dfs(y,x);
}
dep[x]=0;
}
int main(){
scanf("%d%d",&n,&m);
if(m>=2*n){puts("Yes");return 0;}
for(int i=1,x,y;i<=m;i++)scanf("%d%d",&x,&y),v[x].push_back(y),v[y].push_back(x);
for(int i=1;i<=n;i++)dfs(i,0);
puts("No");
return 0;
}
```
下面我们祭出最终优化。
明显，$m\leq2n-3$ 这个限制是极松的。能否找到更紧的限制？
我们来考虑任何一张无长度相等的环的图。考虑极端情况，此时 $3\sim n$ 的环各出现一次。
我们此时来删掉一些边。考虑对于每条边，都有 $\dfrac{1}{\sqrt{n}}$ 的概率将其删掉。则因为 $m,n$ 同级，所以期望删掉的边的条数是 $O\Big(\dfrac{m}{\sqrt{n}}\big)=O(\sqrt{n})$ 级别的。
现在，考虑一长度为 $l$ 的环。在删边结束后，该环幸存的概率是 $\Big(1-\dfrac{1}{\sqrt{n}}\Big)^l$。显然，每个环剩余的概率都是独立的，故全图剩余环的期望就是 $\sum\limits_{l=3}^n\Big(1-\dfrac{1}{\sqrt{n}}\Big)^l$，是等比数列求和的形式。又因为公比 $q=1-\dfrac{1}{\sqrt{n}}$，故 $q\in(0,1)$，所以该式 $<\dfrac{1}{1-q}=\sqrt{n}$。
因此，我们在删掉期望 $\sqrt{n}$ 条边后，期望剩余 $\sqrt{n}$ 个环。又因为期望涵盖了一切可能，所以定存在一种删边方案，使得删边后剩余环数量不超过 $\sqrt{n}$。此时，再删去至多 $\sqrt{n}$ 条边即可消去剩余的 $\sqrt{n}$ 个环。两者结合起来，即对于任一无长度相等的环的图，删去至多 $2\sqrt{n}$ 条边即可消去一切环，换句话说，得到一个边数至多为 $n-1$ 的森林。正因如此，我们就证明了无长度相等的环的图具有的边数至多是 $n+2\sqrt{n}$ 级别的。这样，对于 $m>n+2\sqrt{n}$，也可以直接判 `Yes` 了。
（事实上，这个界也不是实的，但对于本题来说够用了）
考虑求出原图的任意生成森林。这时，非树边至多就只有 $2\sqrt{n}$ 条，所牵涉到的点最多只有 $4\sqrt{n}$ 个。对这 $4\sqrt{n}$ 个点建出虚树，然后虚树上的边便带了权。再把非树边加入虚树，便得到了一张点数至多 $4\sqrt{n}$ 的图。同时，因为虚树是一棵树，边数是 $\sqrt{n}$ 级别的，再加入 $\sqrt{n}$ 级别的非树边，总边数仍是 $\sqrt{n}$ 级别的。这样，我们便得到了一张点、边均为 $\sqrt{n}$ 级别的带权图。明显我们上述 $n^4$ 的算法对带权图仍然适用，因此复杂度便被优化到了 $n^2$。
需要注意的是，在加入非树边后，该虚树便可能出现**重边**。这时不能简单记录环上上一个节点，而应记录环上上一条边。
同时，因为边带了权，长度为 $x$ 的环上不一定就有恰好 $x$ 个点，也就不一定出现恰 $2x$ 次了。为了避免这种情况，我们记录当前长度为 $x$ 的环上点数，若搜出的新环的长度与之不同，明显是两个不同的环，可以直接判 `Yes` 了；否则，因为知晓了点数，就和之前一样统计即可。
时间复杂度 $O(n^2)$。因为压根跑不满，再加上LOJ逆天的评测姬，轻松跑过。期望得分 $100\%$。
.代码：
```cpp
//Version: correct solution (n^2)
#include<bits/stdc++.h>
using namespace std;
int n,m;
namespace OriginalGraph{
int dsu[10100];
int find(int x){return dsu[x]==x?x:dsu[x]=find(dsu[x]);}
bool merge(int x,int y){
x=find(x),y=find(y);
if(x==y)return false;
dsu[x]=y;return true;
}
}
using namespace OriginalGraph;
namespace ImprovedGraph{
int tms[10100],sz[10100],dep[310],ped[310],vis[310],id,head[310],cnt;
struct node{int to,next,val;}edge[40100];
void ae(int u,int v,int w){
// printf("%d %d %d\n",u,v,w);
edge[cnt].next=head[u],edge[cnt].to=v,edge[cnt].val=w,head[u]=cnt++;
edge[cnt].next=head[v],edge[cnt].to=u,edge[cnt].val=w,head[v]=cnt++;
}
bool reach(int x,int fe,int O){
vis[x]=O;
for(int i=head[x],y;i!=-1;i=edge[i].next)if((i^fe)!=1){
y=edge[i].to;
if(!dep[y])return true;
if(dep[y]!=-1||vis[y]==O)continue;
if(reach(y,i,O))return true;
}
return false;
}
void circ(int x,int fe){
if(!reach(x,fe,++id))return;
for(int i=head[x],y;i!=-1;i=edge[i].next)if((i^fe)!=1){
y=edge[i].to;
if(dep[y]!=-1){
if(dep[y])continue;
if(sz[dep[x]+edge[i].val]){
if(ped[x]+1!=sz[dep[x]+edge[i].val]){puts("Yes");exit(0);}
tms[dep[x]+edge[i].val]++;
if(tms[dep[x]+edge[i].val]>2*sz[dep[x]+edge[i].val]){puts("Yes");exit(0);}
}else sz[dep[x]+edge[i].val]=ped[x]+1,tms[dep[x]+edge[i].val]++;
continue;
}
dep[y]=dep[x]+edge[i].val,ped[y]=ped[x]+1,circ(y,i),dep[y]=-1,ped[y]=0;
}
}
}
using namespace ImprovedGraph;
namespace SpanningTree{
vector<int>v[10100];
vector<pair<int,int> >e;
int dis[10100],num[10100],tot,las[10100];
void dfs(int x,int fa){
// printf("%d\n",x);
dis[x]=dis[fa]+1;
if(las[x])num[x]=++tot;
for(auto y:v[x]){
if(y==fa)continue;
dfs(y,x);
if(!las[y])continue;
if(las[x]){
if(las[x]==x)ae(num[x],num[las[y]],dis[las[y]]-dis[x]);
else num[x]=++tot,ae(num[x],num[las[x]],dis[las[x]]-dis[x]),ae(num[x],num[las[y]],dis[las[y]]-dis[x]),las[x]=x;
}else las[x]=las[y];
}
}
}
using namespace SpanningTree;
int main(){
scanf("%d%d",&n,&m),memset(head,-1,sizeof(head));
if(m>=2*sqrt(n)+n){puts("Yes");return 0;}
for(int i=1;i<=n;i++)dsu[i]=i;
for(int i=1,x,y;i<=m;i++){
scanf("%d%d",&x,&y);
if(merge(x,y))v[x].push_back(y),v[y].push_back(x);
else e.push_back(make_pair(x,y)),las[x]=x,las[y]=y;
}
for(int i=1;i<=n;i++)if(!dis[i])dfs(i,0);
for(auto i:e)ae(num[i.first],num[i.second],1);
memset(dep,-1,sizeof(dep));
for(int i=1;i<=tot;i++)dep[i]=0,circ(i,-1),dep[i]=-1;
puts("No");
return 0;
}
```
# VI.圆方树
圆方树是一种神奇的建图方式，它在处理**图上的一条路径**时有奇效。
它的构建方式是这样的：
对于一张无向图，搜出图中的所有**点双连通分量**。注意到割点可能同时存在于多个点双连通分量中。
然后，对于每个点双，我们新建一个点来代表这整个点双。这个新建的点，被称作**方点**，而原树中的点被称作**圆点**。在圆方树中，**方点会连向它代表的全部圆点**。
贴一张图即可。
树建好了，剩下的就是怎么应用了。
## I.[CF487E Tourists](https://www.luogu.com.cn/problem/CF487E)
用这题作圆方树的入门题还是很合适的。
首先，先建出圆方树出来。我们可以给方点赋一个权值，即为它连着的所有圆点的权值的$\min$。然后只需要在圆方树上查询路径$\min$即可。使用树剖即可。
但这个做法会被叉掉：当原图是一张菊花图时，花心的圆点将会连向$n-1$个方点，也就意味着修改该点权值时会影响$n-1$个方点，就可以把复杂度叉到$O(n^2)$。
有什么办法呢？
考虑到树的性质。我们可以在每个方节点开一个```std::multiset<int>```，储存它所有儿子的值。则方点的权值即为```multiset```中的最小值。这样，修改一个圆点时，就只需要改动它父亲的```multiset```即可。
然后，在询问时，如果路径的LCA是一个方点，将答案与该方点的父亲的权值取$\min$即可。
这个算法利用了**圆方树中的圆点会被很多方点统计**的性质，最终让每个圆点仅被一个方点计入，不重不漏，非常神仙。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
int n,m,q,c,val[200100];
namespace SCT{//square circle tree?
int dfn[200100],rev[200100],son[200100],sz[200100],fa[200100],top[200100],dep[200100],tot,mn[800100];
vector<int>v[200100];
void ae(int x,int y){
v[x].push_back(y),v[y].push_back(x);
}
multiset<int>s[100100];
void dfs1(int x,int Fa){
fa[x]=Fa,dep[x]=dep[Fa]+1,sz[x]=1;
for(auto y:v[x]){
if(y==Fa)continue;
dfs1(y,x);
if(x>n)s[x-n].insert(val[y]);
sz[x]+=sz[y];
if(sz[y]>sz[son[x]])son[x]=y;
}
if(x>n)val[x]=*s[x-n].begin();
}
void dfs2(int x){
dfn[x]=++tot,rev[tot]=x;
if(son[x])top[son[x]]=top[x],dfs2(son[x]);
for(auto y:v[x])if(y!=fa[x]&&y!=son[x])top[y]=y,dfs2(y);
}
#define lson x<<1
#define rson x<<1|1
#define mid ((l+r)>>1)
void pushup(int x){mn[x]=min(mn[lson],mn[rson]);}
void build(int x,int l,int r){
if(l==r){mn[x]=val[rev[l]];return;}
build(lson,l,mid),build(rson,mid+1,r),pushup(x);
}
int query(int x,int l,int r,int L,int R){
if(l>R||r<L)return 0x3f3f3f3f;
if(L<=l&&r<=R)return mn[x];
return min(query(lson,l,mid,L,R),query(rson,mid+1,r,L,R));
}
void update(int x,int l,int r,int P){
if(l>P||r<P)return;
if(l==r){mn[x]=val[rev[P]];return;}
update(lson,l,mid,P),update(rson,mid+1,r,P),pushup(x);
}
#undef lson
#undef rson
#undef mid
void prepare(){
dfs1(1,0),top[1]=1,dfs2(1);
// for(int x=1;x<=c;x++)printf("%d::FA:%d SN:%d SZ:%d DP:%d TP:%d DN:%d RV:%d\n",x,fa[x],son[x],sz[x],dep[x],top[x],dfn[x],rev[x]);
build(1,1,c);
}
int ask(int x,int y){
int res=0x3f3f3f3f;
while(top[x]!=top[y]){
if(dep[top[x]]<dep[top[y]])swap(x,y);
res=min(res,query(1,1,c,dfn[top[x]],dfn[x])),x=fa[top[x]];
}
if(dep[x]>dep[y])swap(x,y);
res=min(res,query(1,1,c,dfn[x],dfn[y]));
if(x>n)res=min(res,val[fa[x]]);
// printf("LCA:%d\n",x);
return res;
}
void modify(int x,int y){
if(fa[x])s[fa[x]-n].erase(s[fa[x]-n].find(val[x]));
val[x]=y;
if(fa[x])s[fa[x]-n].insert(val[x]),val[fa[x]]=*s[fa[x]-n].begin(),update(1,1,c,dfn[fa[x]]);
update(1,1,c,dfn[x]);
// for(int i=1;i<=c;i++)printf("%d:%d ",i,val[i]);puts("");
}
}
namespace Graph{
vector<int>v[100100];
int dfn[100100],low[100100],tot;
stack<int>s;
void Tarjan(int x){
dfn[x]=low[x]=++tot,s.push(x);
for(auto y:v[x]){
if(!dfn[y]){
Tarjan(y);
low[x]=min(low[x],low[y]);
if(low[y]>=dfn[x]){
c++;
while(s.top()!=y)SCT::ae(c,s.top()),s.pop();
SCT::ae(c,s.top()),s.pop();
SCT::ae(c,x);
}
}else low[x]=min(low[x],dfn[y]);
}
}
}
char s[10];
int main(){
scanf("%d%d%d",&n,&m,&q),c=n;
for(int i=1;i<=n;i++)scanf("%d",&val[i]);
for(int i=1,x,y;i<=m;i++)scanf("%d%d",&x,&y),Graph::v[x].push_back(y),Graph::v[y].push_back(x);
Graph::Tarjan(1);
SCT::prepare();
for(int i=1,x,y;i<=q;i++){
scanf("%s%d%d",s,&x,&y);
if(s[0]=='C')SCT::modify(x,y);
else printf("%d\n",SCT::ask(x,y));
}
return 0;
}
```
## II.[[APIO2018] Duathlon 铁人两项](https://www.luogu.com.cn/problem/P4630)
我们考虑对于这样一个三元组$\left<s,c,f\right>$，假如我们固定了$s$和$f$，$c$有多少种可能的取值呢？
显然，$c$的取值等于$s\rightarrow f$的简单路径的并集的大小减$2$，因为$s$和$f$不能作为$c$。
那这个并集的大小呢？
我们先假设$s$与$f$位于同一个点双里面，则它们之间必定存在一条路径经过任何一个点$c$。故实际上并集大小就是点双大小。
现在它们不在同一个点双中，则并集大小则为一路上经过的所有点双的大小之和，减去割点的大小（它们被包括在两个点双之中），再减去$2$。
则我们发现，如果将方点的点权赋为它代表的点双的大小，圆点的点权赋为$-1$，则圆方树上$s\rightarrow f$路径上所有点权之和，就等于（并集大小减二）。
因此我们实际上要求的就是所有不同的圆点点对间点权之和。
我们考虑对于路径上的一个点$c$，它究竟被包含到多少条路径中即可。这个简单DP一下就行了。
注意这里的路径是**有向路径**，故记得乘二。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll res;
int n,m,c,all;
namespace RST{
vector<int>v[200100];
int val[200100],sz[200100];
void dfs(int x,int fa){
sz[x]=(x<=n);
for(auto y:v[x])if(y!=fa)dfs(y,x),res+=2ll*sz[y]*sz[x]*val[x],sz[x]+=sz[y];
res+=2ll*sz[x]*(all-sz[x])*val[x];
}
void ae(int x,int y){
// printf("%d %d\n",x,y);
v[x].push_back(y),v[y].push_back(x),val[x]++;
}
}
namespace Graph{
vector<int>v[100100];
int dfn[100100],low[100100],tot;
stack<int>s;
void Tarjan(int x){
all++;
dfn[x]=low[x]=++tot,s.push(x);
for(auto y:v[x]){
if(!dfn[y]){
Tarjan(y),low[x]=min(low[x],low[y]);
if(low[y]>=dfn[x]){
c++;
while(s.top()!=y)RST::ae(c,s.top()),s.pop();
RST::ae(c,s.top()),s.pop();
RST::ae(c,x);
}
}else low[x]=min(low[x],dfn[y]);
}
}
}
int main(){
scanf("%d%d",&n,&m),c=n;
for(int i=1,x,y;i<=m;i++)scanf("%d%d",&x,&y),Graph::v[x].push_back(y),Graph::v[y].push_back(x);
for(int i=1;i<=n;i++)RST::val[i]=-1;
for(int i=1;i<=n;i++){
if(Graph::dfn[i])continue;
all=0;
Graph::Tarjan(i),Graph::s.pop();
RST::dfs(i,0);
}
printf("%lld\n",res);
return 0;
}
```
## III.[[SDOI2018]战略游戏](https://www.luogu.com.cn/problem/P4606)
~~这题我居然能1A，神奇，神奇~~
本题是老缝合怪了，强行把一个**圆方树**板子跟一个**虚树**板子缝到了一起。不会虚树的可以参见笔者的[虚树学习笔记](https://www.luogu.com.cn/blog/Troverld/xu-shu-xue-xi-bi-ji)。
具体来说，首先我们先建出圆方树出来；然后，再在圆方树上针对给定的点集跑出虚树出来；然后，对圆方树上的圆点数量做一个**树上前缀和**（本题只能删圆点，不能删方点），在虚树上每一条边用儿子的前缀和减去父亲的前缀和即可。注意到如果虚树中的虚点（即LCA节点）是圆点的话，则它也是可以被删去的，应该被计入答案。
代码（```namespace```套```namespace```真有意思）：
```cpp
#include<bits/stdc++.h>
using namespace std;
int T,n,m,q,c;
namespace Tree{
int dep[200100],sum[200100],anc[200100][20],dfn[200100],tot;
namespace real{
vector<int>v[200100];
void ae(int x,int y){v[x].push_back(y),v[y].push_back(x);}
void dfs(int x,int fa){
sum[x]=sum[fa]+(x<=n),anc[x][0]=fa,dep[x]=dep[fa]+1,dfn[x]=++tot;
for(int i=1;i<=19;i++)anc[x][i]=anc[anc[x][i-1]][i-1];
for(auto y:v[x])if(y!=fa)dfs(y,x);
}
}
int LCA(int x,int y){
if(dep[x]>dep[y])swap(x,y);
for(int i=19;i>=0;i--)if(dep[x]<=dep[y]-(1<<i))y=anc[y][i];
if(x==y)return x;
for(int i=19;i>=0;i--)if(anc[x][i]!=anc[y][i])x=anc[x][i],y=anc[y][i];
return anc[x][0];
}
namespace imag{
int a[100100],stk[200100],tp,sz[200100],all,res;
bool marked[200100];
vector<int>v[200100];
bool cmp(int x,int y){return dfn[x]<dfn[y];}
void ae(int x,int y){v[x].push_back(y);}
void ins(int x){
marked[x]=true,sz[x]=1;
if(!tp){stk[++tp]=x;return;}
int lca=LCA(x,stk[tp]);
while(tp>=2&&dep[stk[tp-1]]>dep[lca])ae(stk[tp-1],stk[tp]),tp--;
if(tp&&dep[stk[tp]]>dep[lca])ae(lca,stk[tp--]);
if(!tp||stk[tp]!=lca)stk[++tp]=lca;
stk[++tp]=x;
}
void fin(){
while(tp>=2)ae(stk[tp-1],stk[tp]),tp--;
tp--;
}
void dfs1(int x){
int msz=0;
for(auto y:v[x]){
dfs1(y);
if(sz[y]<all)res+=sum[anc[y][0]]-sum[x];
sz[x]+=sz[y],msz=max(msz,sz[y]);
}
if(!marked[x]&&msz<sz[x]&&x<=n)res++;
}
void dfs2(int x){
sz[x]=marked[x]=0;
for(auto y:v[x])dfs2(y);
v[x].clear();
}
int work(){
scanf("%d",&all),res=0;
for(int i=1;i<=all;i++)scanf("%d",&a[i]);
sort(a+1,a+all+1,cmp);
if(a[1]!=1)stk[tp=1]=1;
for(int i=1;i<=all;i++)ins(a[i]);
fin();
dfs1(1),dfs2(1);
return res;
}
}
}
namespace Graph{
vector<int>v[100100];
int dfn[100100],low[100100],tot;
stack<int>s;
void Tarjan(int x){
dfn[x]=low[x]=++tot,s.push(x);
for(auto y:v[x]){
if(!dfn[y]){
Tarjan(y),low[x]=min(low[x],low[y]);
if(low[y]>=dfn[x]){
Tree::real::v[++c].clear();
while(s.top()!=y)Tree::real::ae(c,s.top()),s.pop();
Tree::real::ae(c,s.top()),s.pop();
Tree::real::ae(c,x);
}
}else low[x]=min(low[x],dfn[y]);
}
}
}
int main(){
scanf("%d",&T);
while(T--){
scanf("%d%d",&n,&m),c=n,Graph::tot=Tree::tot=0;
for(int i=1;i<=n;i++)Graph::v[i].clear(),Tree::real::v[i].clear(),Graph::dfn[i]=Graph::low[i]=0;
for(int i=1,x,y;i<=m;i++)scanf("%d%d",&x,&y),Graph::v[x].push_back(y),Graph::v[y].push_back(x);
Graph::Tarjan(1),Graph::s.pop();
Tree::real::dfs(1,0);
scanf("%d",&q);
while(q--)printf("%d\n",Tree::imag::work());
}
return 0;
}
```
## IV.[[GYM102900K]Traveling Merchant](https://codeforces.ml/gym/102900/problem/K)
首先，观察到路径一定是一个 $\rho$ 形的东西，其中在 $\rho$ 的交点之前，一直都是黑白点交替，到了交点处是两个同色点。
于是我们就只保留异色边建一张图，则问题就转变为给你多对同色点，询问有无从 $1$ 经过其中一个点到达另一个点的简单路径。
考虑建出异色边图的圆方树。则我们只需以 $1$ 为根遍历圆方树，将询问的两个圆点升到它们的父亲——两个方点，然后判断它们是否成祖孙关系即可。
时间复杂度 $O(n)$。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
int T,n,m,c;
char s[200100];
namespace Tree{
vector<int>v[400100];
void ae(int x,int y){v[x].push_back(y),v[y].push_back(x);}
int dfn[400100],sz[400100],tot,fa[400100];
void dfs(int x,int P){
fa[x]=P,dfn[x]=++tot,sz[x]=1;
for(auto y:v[x])if(y!=P)dfs(y,x),sz[x]+=sz[y];
}
}
namespace Graph{
vector<int>v[200100];
int dfn[200100],low[200100],tot;
stack<int>s;
void Tarjan(int x){
dfn[x]=low[x]=++tot,s.push(x);
for(auto y:v[x]){
if(!dfn[y]){
Tarjan(y),low[x]=min(low[x],low[y]);
if(low[y]>=dfn[x]){
c++;
while(s.top()!=y)Tree::ae(c,s.top()),s.pop();
Tree::ae(c,s.top()),s.pop();
Tree::ae(c,x);
}
}else low[x]=min(low[x],dfn[y]);
}
}
}
vector<pair<int,int> >q;
int main(){
scanf("%d",&T);
while(T--){
scanf("%d%d%s",&n,&m,s+1),c=n;
for(int i=1,x,y;i<=m;i++){
scanf("%d%d",&x,&y),x++,y++;
if(s[x]!=s[y])Graph::v[x].push_back(y),Graph::v[y].push_back(x);
else q.push_back(make_pair(x,y));
}
Graph::Tarjan(1),Tree::dfs(1,0);
bool ok=false;
for(auto i:q){
if(!Graph::dfn[i.first]||!Graph::dfn[i.second])continue;
int x=i.first,y=i.second;
if(Tree::fa[x])x=Tree::fa[x];
if(Tree::fa[y])y=Tree::fa[y];
if(Tree::dfn[x]>Tree::dfn[y])swap(x,y);
if(Tree::dfn[x]+Tree::sz[x]-1>=Tree::dfn[y]){ok=true;break;}
}
puts(ok?"yes":"no");
for(int i=1;i<=n;i++)Graph::dfn[i]=Graph::low[i]=0,Graph::v[i].clear();Graph::tot=0;
for(int i=1;i<=c;i++)Tree::dfn[i]=Tree::sz[i]=Tree::fa[i]=0,Tree::v[i].clear();c=0;
while(!Graph::s.empty())Graph::s.pop();
q.clear();
}
return 0;
}
```
## V.[[USACO17DEC]Push a Box P](https://www.luogu.com.cn/problem/P4082)
思想很简单，发现任意推动箱子的时刻牛总在箱子旁，而这总共是 $4nm$ 种状态，可以建图储存，然后在上面搜索，搜出所有从起始状态可以到达的状态即可。我们需要连的边只有牛推了一格箱子的边（这个非常简单）以及牛不推箱子，从箱子的一方走到另一方的边。
于是我们要支持查询从一点到另一点不经过指定的某一点是否可行。一开始我的想法是建出圆方树后判定第三点是否在前两点的路径上。这个可以通过求两点间路径长度来判定。但是这样如果要保证复杂度的话，LCA就成为了瓶颈。$O(n)-O(1)$ LCA是存在的，但是我也不会。后来想想看每次询问的两个点应该不会隔太远，所以就暴力上了发树剖。当然是TLE了，开O2后又变成了MLE，于是又费尽心思用指针垃圾回收不用的数组，终于不MLE了，又回到了TLE。
后来想想看，实际上就等价于判定二者是否在一个点双里。判定点双，只需要建出圆方树后，判断其是否是兄弟（有着共同的方点父亲）或者爷孙（方点是一个的儿子，一个的父亲）关系即可。于是把树剖删了，还是TLE。
然后继续卡。在建圆方树的父子关系时，我是用 `vector` 存边然后搜父亲的；但是实际上只需要在搜出一个点双后，令割点为方点父亲，其余点为方点儿子即可。这样就省了一个 `vector` 和一次搜索。吸氧后终于不TLE了，然后WA11。
数据下来一看，发现自己之前代码中有个错误，在计算初始状态时，因为起点和终点不一定紧贴，所以上述思考不一定正确，有可能起点与终点相邻的格子并不在同一个点双内，但是它们间唯一的路径也并不需要经过终点。所以就又加了发爆搜，卡过去了。
时间复杂度 $O(nm)$。常数大得要死。
代码：
```cpp
#include<bits/stdc++.h>
using namespace std;
int n,m,q,dx[4]={1,0,-1,0},dy[4]={0,1,0,-1},cnt;
char s[1510][1510];
#define valid(x,y,z) (x+dx[z]>=0&&x+dx[z]<n&&y+dy[z]>=0&&y+dy[z]<m&&s[x+dx[z]][y+dy[z]]!='#')
int id(int x,int y,int z){return z*n*m+x*m+y;}
int id(int x,int y){return x*m+y;}
#define ID(x,y,z) (x+dx[z])*m+(y+dy[z])
vector<int>v[9000010];
int fa[4500010],dfn[2250010],low[2250010],stk[2250010],tp,tot;
bool VCC(int x,int y){return fa[x]==fa[y]||fa[x]!=-1&&fa[fa[x]]==y||fa[y]!=-1&&fa[fa[y]]==x;}
void Tarjan(int x){
dfn[x]=low[x]=++tot,stk[++tp]=x;
for(auto y:v[x]){
if(!dfn[y]){
Tarjan(y),low[x]=min(low[x],low[y]);
if(dfn[x]>low[y])continue;
while(stk[tp]!=y)fa[stk[tp--]]=cnt;
fa[stk[tp--]]=cnt;
fa[cnt++]=x;
}else low[x]=min(low[x],dfn[y]);
}
}
void init(){
for(int i=0;i<n;i++)for(int j=0;j<m;j++){
if(s[i][j]=='#')continue;
for(int k=0;k<4;k++)if(valid(i,j,k))v[id(i,j)].push_back(ID(i,j,k));
}
for(int i=0;i<n*m;i++)if(!dfn[i])Tarjan(i),tp--;
}
int ax,ay,bx,by;
bool vis[9000010];
void dfs(int x){
vis[x]=true;
for(auto y:v[x])if(!vis[y])dfs(y);
}
int col[1510][1510],ccc;
void dfs(int x,int y){
col[x][y]=ccc;
for(int z=0;z<4;z++)if(valid(x,y,z)&&!col[x+dx[z]][y+dy[z]])dfs(x+dx[z],y+dy[z]);
}
int main(){
// freopen("P4082_6.out","w",stdout);
scanf("%d%d%d",&n,&m,&q),cnt=n*m;
for(int i=0;i<n;i++)scanf("%s",s[i]);
memset(fa,-1,sizeof(fa)),init();
for(int i=0;i<=4*n*m;i++)v[i].clear();
for(int i=0;i<n;i++)for(int j=0;j<m;j++)if(s[i][j]!='#'&&!col[i][j])ccc++,dfs(i,j);
for(int i=0;i<n;i++)for(int j=0;j<m;j++){
if(s[i][j]=='#')continue;
for(int k=0;k<4;k++)if(valid(i,j,k)&&valid(i,j,k^2))v[id(i,j,k)].push_back(k*n*m+(i+dx[k^2])*m+(j+dy[k^2]));
for(int k=0;k<4;k++)for(int p=k+1;p<4;p++){
if(!valid(i,j,k)||!valid(i,j,p))continue;
if(col[i+dx[k]][j+dy[k]]!=col[i+dx[p]][j+dy[p]])continue;
if(!VCC(ID(i,j,k),ID(i,j,p)))continue;
v[id(i,j,k)].push_back(id(i,j,p));
v[id(i,j,p)].push_back(id(i,j,k));
}
if(s[i][j]=='A')ax=i,ay=j;
if(s[i][j]=='B')bx=i,by=j;
}
memset(col,0,sizeof(col)),ccc=1;
s[bx][by]='#';
dfs(ax,ay);
s[bx][by]='B';
for(int i=0;i<4;i++){
if(!valid(bx,by,i))continue;
if(!col[bx+dx[i]][by+dy[i]])continue;
v[4*n*m].push_back(id(bx,by,i));
}
dfs(4*n*m);
for(int i=1,x,y;i<=q;i++){
scanf("%d%d",&x,&y),x--,y--;
bool ok=false;
for(int k=0;k<4;k++)if(valid(x,y,k)&&vis[id(x,y,k)]){ok=true;break;}
if(x==bx&&y==by)ok=true;
puts(ok?"YES":"NO");
}
return 0;
}
```