P6624 [省选联考 2020 A 卷] 作业题

大力推式子+矩阵树题

简记 $\gcd(w_1,w_2,...,w_{n-1})$ 为 $G(W)$

$$\sum_{T\in E} \left(\sum_i^{n-1} w_i\right)\times G(W)\\ =\sum_g^M g\sum_{Tree\in E} \left(\sum_i^{n-1} w_i\right)\times[G(W)=g]\\ =\sum_g^M g\sum_{Tree\in E} \left(\sum_i^{n-1} w_i\right)\times\left(\prod_i^{n-1}[g\vert w_i]\right)\times[G(\frac{W}{g})=1]\\ =\sum_g^M g\sum_{Tree\in E} \left(\sum_i^{n-1} w_i\right)\times\left(\prod_i^{n-1}[g\vert w_i]\right)\times \sum_{d\vert G(\frac{W}{g})}\mu(d)\\ =\sum_g^M g\sum_d^M\mu(d)\sum_{Tree\in E} \left(\sum_i^{n-1} w_i\right)\times\left(\prod_i^{n-1}[g\vert w_i]\right)\times [d\vert G(\frac{W}{g})]\\ =\sum_g^M \sum_d^Mg\mu(d)\sum_{Tree\in E} \left(\sum_i^{n-1} w_i\right)\times\left(\prod_i^{n-1}[g\vert w_i]\right)\times \left(\prod_i^{n-1}[d\vert \frac{w_i}{g}]\right)\\ =\sum_g^M \sum_d^Mg\mu(d)\sum_{Tree\in E} \left(\sum_i^{n-1} w_i\right)\times\left(\prod_i^{n-1}[gd\vert w_i]\right)\\ =\sum_T^M \sum_{d\vert T}\mu(d)\left(\frac{T}{d}\right)\sum_{Tree\in E} \left(\sum_i^{n-1} w_i\right)\times\left(\prod_i^{n-1}[T\vert w_i]\right)\\ $$

令 $F(x)=\sum_{d\vert x}\mu(x)\dfrac{x}{d}$，显然 $F(x)$ 可以在筛出 $\mu$ 之后 $O(n\ln n)$ 地求得。所以原式变成：

$$\sum_T^MF(T)\times\sum_{Tree\in E} \left(\sum_i^{n-1} w_i\right)\times\left(\prod_i^{n-1}[T\vert w_i]\right)\\ $$

由于边数不多，我们可以对于每个 $T$ 暴力保留那些边权是 $T$ 的倍数的边。如果少于 $n-1$ 条说明不存在任何生成树。

现在转化为给定一个无向图 $G$，如何求出：

$$\sum_{Tree\in E}\left(\sum_i^{n-1}w_i\right)$$

发现我们会做的其实是乘积和：

$$\sum_{Tree\in E}\left(\prod_i^{n-1}w_i\right)$$

事实上存在一种重新更改边权的方式，使得计算乘积之后可以得到和。

考虑如下式子：

$$(1+bx)(1+dx)=1+(b+d)x+bdx^2$$

发现一次项的变化就是 $b+d=(b+d)$，因此我们把每条边的边权赋值为一个多项式 $1+w_ix$，之后正常做乘积合即可。因为我们只需要边权和，所以二次项可以全部舍去，即在 $mod x^2$ 意义下做多项式四则运算。

多项式的加减乘都是简单的，讨论如何做除法，即多项式求逆。

尝试计算下面这个除法算式：

$$\dfrac{a+bx}{c+dx}\ \mod x^2$$

先求出 $\dfrac{1}{c+dx}$，再将其乘上 $a+bx$。

设 $C+Dx=\dfrac{1}{c+dx}$，则根据定义有 $(C+Dx)(c+dx)=1$，对应项拆开后得到 $Cc=1$ 和 $Cd+Dc=0$，解得 $C=\dfrac{1}{c}$，$D=-\frac{d}{c^2}$。再乘上 $a+bx$得到公式：

$$\dfrac{a+bx}{c+dx}=\dfrac{a}{c}+\dfrac{bc-ad}{c^2}x\ \mod x^2$$

如果 $c$ 为 $0$ 怎么办？在高斯消元的过程中，可以先寻找 $i$ 这一列有没有非 $0$ 常数项。如果有，将其交换至上方作为除数；如果没有，则说明所有多项式都无常数项，直接 $\dfrac{bx}{dx}=\dfrac{b}{d}$ 即可。

// Problem: P6624 [省选联考 2020 A 卷] 作业题
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P6624
// Memory Limit: 512 MB
// Time Limit: 3000 ms
// Genshin Impact launch
// ZYB TXDY
// 
// Powered by CP Editor (https://cpeditor.org)

#include <bits/stdc++.h>
#define INF 0x7fffffff
#define MAXN 200005
#define eps 1e-9
#define foru(a,b,c)	for(int a=b;a<=c;a++)
#define RT return 0;
#define LL long long
#define LXF int
#define RIN rin()
#define HH printf("\n")
#define all(x) (x).begin(),(x).end()
#define fi first
#define se second
using namespace std;
inline LXF rin(){
	LXF x=0,w=1;
	char ch=0;
	while(ch<'0'||ch>'9'){ 
	if(ch=='-') w=-1;
	ch=getchar();
	}
	while(ch>='0'&&ch<='9'){
	x=x*10+(ch-'0');
	ch=getchar();
	}
	return x*w;
}
const LL mod=998244353;
LL qpow(LL a,LL b){
	a%=mod;
	LL ret=1;
	while(b){
		if(b&1ll)	ret=ret*a%mod;
		b>>=1ll,a=a*a%mod;
	}
	return ret;
}
int n,m,M;
vector<int> prim;
int mu[MAXN];
bool vis[MAXN];
LL F[MAXN];
bool ok[MAXN];
void pre(){
	mu[1]=1;
	foru(i,2,M){
		if(!vis[i]){
			prim.push_back(i);
			mu[i]=-1;
		}
		for(int j=0;j<prim.size() && i*prim[j]<=M;j++){
			vis[i*prim[j]]=1;
			if(i%prim[j]==0){
				mu[i*prim[j]]=0;
				break;
			}
			mu[i*prim[j]]=-mu[i];
		}
	}
	for(int d=1;d<=M;d++){
		for(int T=d;T<=M;T+=d){
			ok[d]|=ok[T];
			F[T]=(F[T]+mu[d]*(T/d)%mod)%mod;
		}
	}
}
class Edge{
	public:
		int u,v,w;
}eg[905];
class Function{
	public:
		Function(LL _a=0,LL _b=0):a(_a),b(_b){}
		LL a,b;
		Function operator + (const Function x)const{
			return Function((a+x.a)%mod,(b+x.b)%mod);
		}
		Function operator - (const Function x)const{
			return Function((a-x.a+mod)%mod,(b-x.b+mod)%mod);
		}
		Function operator * (const Function x)const{
			return Function(a*x.a%mod,(a*x.b%mod+b*x.a%mod)%mod);
		}
		Function operator / (const Function x)const{
			if(a!=0 && x.a==0)	exit(-1);
			if(a==0 && x.a==0)	return Function(0,b*qpow(x.b,mod-2)%mod);
			LL C=qpow(x.a,mod-2);
			// cout<<a<<' '<<a*C%mod<<endl;
			return Function(a*C%mod,(b*C%mod-a*x.b%mod*C%mod*C%mod+mod)%mod);
		}
}a[35][35];
LL det(){
	Function ret=Function(1,0);
	// cout<<a[2][1].a<<endl;
	foru(i,1,n-1){
		int p=0;
		foru(j,i,n-1){
			if(a[j][i].a!=0){
				p=j;
				break;
			}
		}
		if(p>i){
			swap(a[i],a[p]);
			ret=Function(0,0)-ret;
		}
		foru(j,i+1,n-1){
			auto t=a[j][i]/a[i][i];
			// cout<<a[j][i].a<<' '<<a[j][i].b<<endl;
			foru(k,i,n-1){
				// cout<<k<<' '<<(t*a[i][k]).a<<endl;
				a[j][k]=a[j][k]-(t*a[i][k]);
			}
		}
		ret=ret*a[i][i];
	}
	return ret.b;
}
signed main(){
	n=RIN,m=RIN;
	foru(i,1,m){
		int u=RIN,v=RIN,w=RIN;
		eg[i]={u,v,w};
		M=max(M,w);
		ok[w]=1;
	}
	pre();
	LL ans=0;
	foru(T,1,M){
		// if(!ok[T])	continue;
		foru(i,1,n){
			vis[i]=0;
			foru(j,1,n)	a[i][j]=Function();
		}
		int ecnt=0;
		foru(i,1,m){
			auto [u,v,w]=eg[i];
			if(w%T!=0)	continue;
			vis[u]=vis[v]=1;
			ecnt++;
			a[u][u]=a[u][u]+Function(1,w);
			a[v][v]=a[v][v]+Function(1,w);
			a[u][v]=a[u][v]-Function(1,w);
			a[v][u]=a[v][u]-Function(1,w);
		}
		bool ct=1;
		foru(i,1,n)	ct&=vis[i];
		if(!ct || ecnt<n-1)	continue;
		
		LL res=det();
		// cout<<ecnt<<' '<<T<<' '<<F[T]<<' '<<res<<endl;
		ans=(ans+F[T]*res%mod)%mod;
	}
	cout<<ans;
	return 0;
}