这里是矩阵树定理学习笔记。线性代数基础详见线性代数学习笔记。
矩阵树定理是用来作生成树计数的好东西。其定理主要表述如下：
设 \(\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 定理
就是模板啦。但是，这题有边权，怎样处理？
我们考虑如果将一条带权边拆作（权值）条无权边，则任何一棵带权边的生成树，其中的一条边，就可以换成任何一条无权边，一共（权值）条，是等效的。依据乘法原理，我们将每条带权边看作无权边，总权值仍然是符合实际意义的。
代码：
#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;
}
考虑容斥套Matrix-Tree。我们考虑全部 \(n-1\) 个公司全部可选的生成树数量，减去选 \(n-2\) 个公司的数量，加上 \(n-3\) 个……于是只需枚举哪些公司然后容斥一下即可。
代码：
#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;
}
这里是本题的矩阵树定理解法。DP解法详见重题——DP刷题笔记III，CXX.CF917D Stranger Trees。
我们考虑为所有树上的边赋权为 \(1\)，完全图中不在树上的边赋权为 \(x\)。对其求矩阵树定理。明显，其形式是一 \(n-1\) 阶矩阵的行列式，也即一关于 \(x\) 的 \(n-1\) 次多项式。则，多项式上常数项系数即为出现 \(0\) 条新边的生成树数，一次项系数即为出现 \(1\) 条新边的生成树数……然后对系数求一个前缀和，即是最多出现 \(k\) 条新边的生成树数。
但问题来了，没人教过我们应该如何对一个有未知项的矩阵求行列式呀（具体是因为消元的时候不知道应该怎样判定“\(0\)”）？
没关系。按照我们之前的结论，其行列式是一 \(n-1\) 次多项式。那我们直接挑 \(n\) 个固定的 \(x\) 各跑一遍矩阵树，然后拉格朗日插一下就插出多项式了。当然，如果你嫌拉格朗日还要手工进行多项式乘除法太麻烦，也可以直接上高斯消元插多项式。
时间复杂度 \(O(n^4)\)。
代码：
#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));
}
*/
首先，这题实际上很板子。但是，目的是通过本题介绍一下模数非质数时的高斯消元。
模数非质数时，无法求逆元；但是要说如何把一个东西消成 \(0\) 的话，可以想到辗转相除法，即\(\gcd\)时要用的那东西。直接上辗转相除，可以被证明复杂度仍是 \(O(n^3)\) 的。
代码：
#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;
}
这里就涉及到矩阵树定理的本质了——你费尽心思，又是建矩阵又是求行列式的，究竟是为了什么呢？
矩阵树定理的本质，是求
而本题要求的，是
考虑提出去一个 \(\prod(1-e_i)\)，得到
这里就可以直接上矩阵树了。
但是，如果我们发现有边权恰为 \(1\)，应该咋办呢？
发现，恰为 \(1\) 的实际意义是两点间必有边，故使用冰茶姬将两者并一块即可。同时，因为还要保证是一棵树，冰茶姬维护的同时还要保证无环。这里的无环包括两种意义的无环，一是 \(1\) 边自身不能成环，二是 \(1\) 边构成的连通块中不能出现其它边（即，其它边的选择都必须是 \(1-e_i\)）。
代码：
#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;
}
我觉得是至今为止刷过的所有矩阵树中最棒的一个。
首先，有一个性质是很明显的：在一张图的所有最小生成树中，其每个权值的边的数量都是固定的。不然，假如有两棵MST，其中存在一个权值，使得该权值的边在两棵树中数量不同，则一定存在一边权和更小的生成树。
有了这个性质，我们就可以再推出一个更强的性质——其中任意一棵MST中，某一权值的边所联通的点集都是相同的，都是连完MST上其他权值的边后尚未联通的点集。
很明显，第一条性质是第二条性质的子集（固定的点集中边的数量也是固定的），故我们只需简单证明一下第二条性质即可。
考虑Kruskal算法的过程。考虑生成树中最大权值的所有边，依据Kruskal算法，在我们未添加该权值的任何边之前，其他边所联通的集合是固定的（假如并非固定，即存在两种不同的方案使得一对点在一种方案上联通而另一种不连通，则我们完全可以取两种方案上所有边的某种并集使得总权值和更小）——这也意味着该权值的所有边联通的集合是固定的。则我们现在将该集合看做全部用最大权值的边联通，然后按照此种方法再去分析次大权值的所有边，即可简单得到第二条性质。
于是我们现在就直接套矩阵树即可——因为每一种边权所联通的集合是固定的，故我们对该集合（当然，用dsu缩去不由该边权的点所联通的点集后）以及该边权的所有边直接上一个矩阵树计算联通该集合的方案数。最后，因为每种边权的集合间独立，上乘法原理将所有边权的方案乘一块即可。
时间复杂度 \(O(n^3+nm+m\log m)\)，因为最小生成树上的边数固定是 \(n-1\)，故大部分操作都没有求行列式的 \(O(n^3)\) 大。
代码：
#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;
}
考虑式子
对 \(\gcd\) 一类东西，常规思路是莫比乌斯反演掉。此处我们考虑欧拉反演。
把 \(d\) 扔到前面去
我们考虑就枚举这个 \(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}\)。但这只是理论上界，实际应用起来跑得飞快。
代码：
#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;
}
