ybtAu「图论」第8章 生成树计数
A. 【例题1】重建
题目要求：
\[\large \sum_T\prod_{(u,v)\in T}G_{u,v}\prod_{(u,v)\not\in T}(1-G_{u,v})
\]
即：
\[\large \sum_T\prod_{(u,v)\in T}G_{u,v}\frac{\prod_{(u,v)}(1-G_{u,v})}{\prod_{(u,v)\in T}(1-G_{u,v})}\\
\large=\prod_{(u,v)}(1-G_{u,v})\sum_T\prod_{(u,v)\in T}\frac{G_{u,v}}{1-G_{u,v}}
\]
令 \(W_{u,v}=\frac{G_{u,v}}{1-G_{u,v}}\)，使用矩阵树定理求解即可。
注意精度。
#include <iostream>
#define N 55
const double eps=1e-7;
int n;
double a[N][N];
double gauss()
{
double ret=1;
for(int i=1;i<n;i++)
{
int mx=i;
for(int j=i+1;j<n;j++) if(abs(a[j][i])>abs(a[mx][i])) mx=j;
if(a[mx][i]<eps) return 0;
if(mx!=i)
{
ret*=-1;
for(int j=1;j<n;j++) std::swap(a[i][j],a[mx][j]);
}
ret*=a[i][i];
for(int j=i+1;j<n;j++)
{
double dt=a[j][i]/a[i][i];
for(int k=i;k<n;k++) a[j][k]-=dt*a[i][k];
}
}
return ret;
}
signed main()
{
std::ios::sync_with_stdio(0);
std::cin.tie(0),std::cout.tie(0);
std::cin>>n;
double prd=1;
for(int i=1;i<=n;i++) for(int j=1;j<=n;j++)
{
double x;
std::cin>>x;
if(x==1) x-=eps;
if(i<j) prd*=1-x;
a[i][j]=-x/(1-x);
}
for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) if(i!=j) a[i][i]-=a[i][j];
std::cout<<gauss()*prd;
}
B. 【例题2】树的计数
这里写一下矩阵树定理的解法。
列出该二分图去掉最后一行与最后一列的 Laplace 矩阵（\(I\) 表示单位矩阵，\(J\) 表示全是 \(1\) 的矩阵）：
\[\large
\begin{vmatrix}
mI_n & -J_{m-1,n}\\
-J_{n,m-1} & nI_{m-1}
\end{vmatrix}
\]
由行列式的性质，用左面 \(n\) 列 \(\times \frac{1}{m}\) 去消右上方矩阵，得：
\[\large
\begin{vmatrix}
mI_n & 0\\
-J_{n,m-1} & nI_{m-1}-\frac{n}{m}J_{m-1,m-1}
\end{vmatrix}
\]
用上面 \(n\) 行 \(\times \frac{1}{m}\) 去消左下方矩阵，得：
\[\large
\begin{vmatrix}
mI_n & 0\\
0 & nI_{m-1}-\frac{n}{m}J_{m-1,m-1}
\end{vmatrix}
\]
左上角的矩阵相对好处理，对于右下角的矩阵：
\[\large
\begin{vmatrix}
n-\frac{n}{m} & -\frac{n}{m} & -\frac{n}{m}\\
-\frac{n}{m} & n-\frac{n}{m} & -\frac{n}{m}\\
-\frac{n}{m} & -\frac{n}{m} & n-\frac{n}{m}
\end{vmatrix}
\]
把右面第 \(2\) 到 \(m-1\) 列都加到左边，第 \(1\) 列全都变成 \(\frac{n}{m}\)；再把第 \(1\) 列加到右面所有列，有：
\[\large
\begin{vmatrix}
\frac{n}{m} & 0 & 0\\
\frac{n}{m} & n & 0\\
\frac{n}{m} & 0 & n
\end{vmatrix}
\]
变成了一个下三角矩阵。
对整个矩阵（包括左上角）求行列式，为 \(m^nn^{m-2}\frac{n}{m}=m^{n-1}n^{m-1}\)。
注意要开 __int128。
#include <iostream>
#define int long long
int n,m,p;
__int128 qpow(__int128 x,int y) {__int128 r=1;for(;y;y>>=1,x=x*x%p) if(y&1) r=r*x%p;return r;}
signed main()
{
std::ios::sync_with_stdio(0);
std::cin.tie(0),std::cout.tie(0);
std::cin>>n>>m>>p;
std::cout<<(int)(qpow(n,m-1)*qpow(m,n-1)%p);
}
C. 作业题
题目要求：
\[\large \sum_T(\sum_{e_i\in T}w_{e_i}\gcd_{e_i\in T}(w_{e_i}))\\
\large =\sum_T(\sum_{e_i\in T}w_{e_i}\sum_{d\mid w_{e_i}}\varphi(d))\\
\large =\sum_d\varphi(d)\sum_{T,\forall e_i\in T,d\mid w_{e_i}}\sum_{e_i\in T}w_{e_i}\\
\]
要求边权和，而矩阵树定理只能用于求解边权积，因此需要对边权进行变换。
令 \(W_{e_i}=w_{e_i}x+1\)，此时 \(\prod_{e_i\in T}W_{e_i}\) 的一次项系数就是边权和。于是可以使用矩阵树定理求解。
关于二项式的运算细节见代码。
#include <iostream>
#define R 152501
#define mod 998244353
#define int long long
#define N 40
int n,m,U[N*N],V[N*N],wt[N*N],pri[R],phi[R+5],ex[R+5],cnt;
int qpow(int x,int y) {int r=1;for(;y;y>>=1,x=x*x%mod) if(y&1) r=r*x%mod;return r;}
struct Line
{
int x,y;
Line(int _x=0,int _y=0):x(_x),y(_y) {}
Line operator +(const Line &g) const {return {(x+g.x)%mod,(y+g.y)%mod};}
Line operator -(const Line &g) const {return {(x-g.x+mod)%mod,(y-g.y+mod)%mod};}
Line operator *(const Line &g) const {return {(y*g.x%mod+x*g.y%mod)%mod,y*g.y%mod};}
Line operator /(const Line &g) const
{
int inv=qpow(g.y,mod-2);
return {(x*g.y%mod-y*g.x%mod+mod)%mod*inv%mod*inv%mod,y*inv%mod};
}
} a[N][N];
void sieve()
{
phi[1]=1;
for(int i=2;i<=R;i++)
{
if(!ex[i]) pri[++cnt]=i,phi[i]=i-1;
for(int j=1;j<=cnt&&i*pri[j]<=R;j++)
{
ex[i*pri[j]]=1;
if(i%pri[j]) phi[i*pri[j]]=phi[i]*(pri[j]-1);
else {phi[i*pri[j]]=phi[i]*pri[j];break;}
}
}
}
Line gauss()
{
Line ret={0,1};
bool fg=0;
for(int i=1;i<n;i++)
{
if(!a[i][i].y) for(int j=i+1;j<n;j++) if(a[j][i].y)
{
fg^=1;
for(int k=1;k<n;k++) std::swap(a[i][k],a[j][k]);
break;
}
ret=ret*a[i][i];
Line inv=Line(0,1)/a[i][i];
for(int j=i+1;j<n;j++)
{
Line dt=a[j][i]*inv;
for(int k=i;k<n;k++) a[j][k]=a[j][k]-dt*a[i][k];
}
}
if(fg) ret={mod-ret.x,mod-ret.y};
return ret;
}
signed main()
{
std::ios::sync_with_stdio(0);
std::cin.tie(0),std::cout.tie(0);
std::cin>>n>>m;
sieve();
for(int i=1;i<=m;i++) std::cin>>U[i]>>V[i]>>wt[i];
int ans=0;
for(int d=1;d<=R;d++)
{
int cnt=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<=m;i++) if(wt[i]%d==0)
{
a[U[i]][U[i]]=a[U[i]][U[i]]+(Line){wt[i],1};
a[V[i]][V[i]]=a[V[i]][V[i]]+(Line){wt[i],1};
a[U[i]][V[i]]=a[U[i]][V[i]]-(Line){wt[i],1};
a[V[i]][U[i]]=a[V[i]][U[i]]-(Line){wt[i],1};
cnt++;
}
if(cnt<n-1) continue;
//for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) printf("%dx+%d%c",a[i][j].x,a[i][j].y," \n"[j==n]);
//printf("Check %d phi %d det %d\n",d,phi[d],gauss().x);
(ans+=phi[d]*gauss().x%mod)%=mod;
}
std::cout<<ans;
}
浙公网安备 33010602011771号