\(\frak{Description}\)
\(\frak{Solution}\)
看到 \(\gcd\) 考虑用 \(\varphi*1=\rm id\) 来化简
\[\begin{align}
val(T)&=\left(\sum_{i=1}^{n-1}w_{e_i}\right)\cdot \gcd(w_{e_1},w_{e_2},\dots,w_{e_{n-1}})\\
&=\left(\sum_{i=1}^{n-1}w_{e_i}\right)\cdot \sum_{d\mid w_{e_1},\dots,d\mid w_{e_{n-1}}}\varphi(d)\\
&=\sum_{d}\varphi(d)\cdot \sum_{T}\left(\sum_{i=1}^{n-1}w_{e_i}\right)\cdot [d\mid w_{e_1},\dots,d\mid w_{e_{n-1}}]
\end{align}
\]
枚举因数，后面的一坨可以用矩阵树定理。但是矩阵树定理是求 "权值积之和" 而不是 "权值和之和"，我们还需要做一个转化：将权值 \(w\) 变成 \(wx+1\)（模 \(x^2\) 意义下），这样最后求出 "权值积之和"，而 "权值和之和" 就是一次项系数！
还有一个问题，计算时涉及到高斯消元，会有元素的加减乘除。如何推导 \(\frac{ax+b}{cx+d}\)（在模 \(x^2\) 意义下）？我们设 \(\frac{1}{cx+d}=fx+g\)，可以推出 \(g=1/d\)，从而得到 \(f=\frac{-c}{d^2}\)，代回去就可以得知 \(\frac{ax+b}{cx+d}=\frac{b}{d}+\frac{ad-bc}{d^2}\cdot x\).
最后还可以做一些剪枝：只计算边数大于等于 \(n-1\) 的矩阵。复杂度就是 \(\mathcal O\left(n^3\cdot \frac{\sum \sigma_0(d)}{n-1}\right)=\mathcal O(n^2\cdot \sum\sigma_0(d))\).
另外这题还有一个拓展：给定一张图，图上每条边是红色或蓝色，求恰有 \(k\) 条红边的生成树个数。\(n\le 50\).
类似地，对于限制条件也可以利用生成函数，把红边边权设为 \(x+0\)，蓝边为 \(0x+1\). 那么最后得到的多项式中 \(x^k\) 的系数就是答案。
\(\frak{Code}\)
# include <cctype>
# include <cstdio>
# define print(x,y) write(x), putchar(y)
template <class T>
inline T read(const T sample) {
T x=0; char s; bool f=0;
while(!isdigit(s=getchar())) f|=(s=='-');
for(; isdigit(s); s=getchar()) x=(x<<1)+(x<<3)+(s^48);
return f? -x: x;
}
template <class T>
inline void write(T x) {
static int writ[50], w_tp=0;
if(x<0) putchar('-'), x=-x;
do writ[++w_tp]=x-x/10*10, x/=10; while(x);
while(putchar(writ[w_tp--]^48), w_tp);
}
# include <cstring>
# include <iostream>
using namespace std;
typedef long long ll;
const int maxv = 152511;
const int mod = 998244353;
inline int inc(int x,int y) { return x+y>=mod?x+y-mod:x+y; }
inline int adj(int x,int y) { return x+y>=mod?x+y-mod:(x+y<0?x+y+mod:x+y); }
inline ll inv(ll x,int y=mod-2) {
ll r=1;
for(; y; y>>=1, x=x*x%mod)
if(y&1) r=r*x%mod;
return r;
}
struct node {
ll x,y;
node() {}
node(int X,int Y):x(X),y(Y) {}
node operator + (const node& t) const { return node((x+t.x)%mod,(y+t.y)%mod); }
node operator - (const node& t) const { return node((x-t.x)%mod,(y-t.y)%mod); }
node operator * (const node& t) const { return node(x*t.x%mod,(x*t.y+y*t.x)%mod); }
node operator / (const node& t) const {
ll Inv = inv(t.x);
return node(x*Inv%mod, (y*t.x-x*t.y)%mod*Inv%mod*Inv%mod);
}
inline void Print() { printf("(%lld, %lld) ",x,y); }
} a[33][33];
struct edge { int u,v,w; } E[900];
bool is[maxv];
int n,m,maxw,cnt[maxv],pc,p[maxv],phi[maxv];
void sieve(int n) {
phi[1]=1;
for(int i=2;i<=n;++i) {
if(!is[i]) p[++pc]=i, phi[i]=i-1;
for(int j=1; j<=pc && i*p[j]<=n; ++j) {
is[i*p[j]] = true;
if(i%p[j]==0) {
phi[i*p[j]] = phi[i]*p[j];
break;
} else phi[i*p[j]] = phi[i]*(p[j]-1);
}
}
}
void div(int x) {
for(int i=1; i*i<=x; ++i)
if(x%i==0) {
++ cnt[i];
if(i*i!=x) ++ cnt[x/i];
}
}
void consGraph(int val) {
memset(a,0,sizeof a);
node t; int x,y;
for(int i=1;i<=m;++i) {
if(E[i].w%val) continue;
t = node(1,E[i].w);
x=E[i].u, y=E[i].v;
a[x][x] = a[x][x]+t;
a[y][y] = a[y][y]+t;
a[x][y] = a[x][y]-t;
a[y][x] = a[y][x]-t;
}
}
node gauss() {
int j; bool f=0; node ret=node(1,0), Inv, mul;
for(int i=1;i<=n-1;++i) {
for(j=i;j<=n-1;++j) if(a[j][i].x) break;
if(j>n-1) return node(0,0);
if(i^j) swap(a[i],a[j]), f^=1;
ret = ret*a[i][i];
Inv = node(1,0)/a[i][i];
for(j=i+1;j<=n-1;++j) {
mul = a[j][i]*Inv;
for(int k=i;k<=n-1;++k)
a[j][k] = a[j][k]-mul*a[i][k];
}
}
return f? node(0,0)-ret: ret;
}
int main() {
n=read(9), m=read(9);
for(int i=1;i<=m;++i) {
E[i].u=read(9), E[i].v=read(9);
div(E[i].w=read(9));
maxw = max(maxw,E[i].w);
}
sieve(maxw); int ans=0;
for(int i=1;i<=maxw;++i) if(cnt[i]>=n-1)
consGraph(i), ans = adj(ans,gauss().y%mod*phi[i]%mod);
print(ans,'\n');
return 0;
}
