P6624 - 作业题 题解

如何评价一年多以后再来补省选题

考虑莫反，枚举 gcd（以下 \(T\) 为生成树，\(G\) 为原图，\(G[x]\) 为只保留边权为 \(x\) 的倍数得到的子图，\(T\) 为 \(P\) 的生成树不严格地记为 \(T\in P\)）：

\[\begin{aligned} ans&=\sum_{T\in G}\sum\limits_{i=1}^{n-1}w_{T_i}\mathop{\gcd}_{j=1}^{n-1}w_{T_j}\\&=\sum_kk\sum_{T\in G[k]}\sum_{i=1}^{n-1}w_{T_i}\left[\mathop{\gcd}_{j=1}^{n-1}\dfrac{w_{T_j}}k=1\right]\\&=\sum_kk\sum_{T\in G[k]}\sum_{i=1}^{n-1}w_{T_i}\sum_{o\mid\frac{w_{T_1}}k,\frac{w_{T_2}}k,\cdots,\frac{w_{T_{n-1}}}k}\mu(o)\\&=\sum_o\mu(o)\sum_kk\sum_{T\in G[ko]}\sum_{i=1}^{n-1}w_{T_i}\\&=\sum_o\mu(o)\sum_kk\operatorname{calc}_{ko} \end{aligned} \]

由于值域有 \(v=152501\) 的良心限制，考虑计算出 \(1\sim v\) 内每个 \(\operatorname{calc}_i\)。那这就是个生成树权值和的和，把一次多项式塞进去跑矩阵树即可。总复杂度 \(\mathrm O\!\left(vn^3\right)\)，4e9 了已经，考虑优化。

显然只要计算那些满足有至少 \(n-1\) 条边边权为 \(i\) 的倍数的 \(\operatorname{calc}_i\)，其它的一个生成树都没有一定是 \(0\)。那这些有多少呢？到 SF 主页里查一下 \(\max\mathrm d\) 大概是 \(200\) 以内，\(m=\mathrm O\!\left(n^2\right)\)，所以把所有边权的因数放到一个可重集里，里面所有重数 \(\geq n-1\) 的数最多有 \(\mathrm O\!\left(\dfrac{200n^2}n\right)=\mathrm O(200n)\) 种。总复杂度 \(\mathrm O\!\left(200n^4\right)\)，1e8 差不多行，实际上不是很跑的满。

最后线对计算 \(ans\) 即可。

code

#include<bits/stdc++.h>
using namespace std;
#define pb push_back
const int mod=998244353;
int qpow(int x,int y){
	int res=1;
	while(y){
		if(y&1)res=1ll*res*x%mod;
		x=1ll*x*x%mod;
		y>>=1;
	}
	return res;
}
int inv(int x){return qpow(x,mod-2);}
const int N=40,M=N*N,V=152510;
vector<int> prm;
bool vis[V];
int mu[V],dld[V+1];
void sieve(){
	mu[1]=1;
	for(int i=2;i<V;i++){
		if(!vis[i]){
			prm.pb(i);
			mu[i]=-1;
			dld[i]=1;
		}
		for(int j=0;j<prm.size();j++){
			int x=prm[j];
			if(i*x>V)break;
			vis[i*x]=true;
			if(i%x){
				mu[i*x]=mu[i]*mu[x];
				dld[i*x]=i;
			}
			else{
				dld[i*x]=dld[i];
				if(dld[i*x]==1)mu[i*x]=0;
				else mu[i*x]=mu[dld[i*x]]*mu[i*x/dld[i*x]];
				break;
			}
		}
	}
}
int n,m;
int fr[M],to[M],w[M];
int cnt[V];
int det[V];
struct poly{
	int a,b;
	poly(int x=0,int y=0){a=x,b=y;}
	friend poly operator+(poly x,poly y){return poly((x.a+y.a)%mod,(x.b+y.b)%mod);}
	friend poly operator-(poly x,poly y){return poly((x.a-y.a)%mod,(x.b-y.b)%mod);}
	friend poly operator*(poly x,poly y){
		return poly((1ll*x.a*y.b+1ll*x.b*y.a)%mod,1ll*x.b*y.b%mod);
	}
	friend poly operator/(poly x,poly y){
		int a=x.a,b=x.b,c=y.a,d=y.b;
		int e=(1ll*a*d-1ll*c*b)%mod*inv(1ll*d*d%mod)%mod,f=1ll*b*inv(d)%mod;
		return poly(e,f);
	}
	void operator+=(poly x){*this=*this+x;}
	void operator-=(poly x){*this=*this-x;}
	void operator*=(poly x){*this=*this*x;}
	void operator/=(poly x){*this=*this/x;}
}a[N][N];
void swp(int x,int y){
	for(int i=1;i<=n;i++)swap(a[x][i],a[y][i]);
}
void tadd(int x,poly v,int y){
	for(int i=1;i<=n;i++)a[y][i]+=v*a[x][i];
}
int gauss(){
	poly ans(0,1);
	for(int i=1;;i++){
		int row=0,col=0;
		for(int j=1;j<=n;j++){
			for(int k=i;k<=n;k++)if(a[k][j].b){row=k,col=j;break;}
			if(row)break;
		}
		if(!row)break;
		swp(i,row);
		poly iv=poly(0,1)/a[i][col];
		for(int j=i+1;j<=n;j++)tadd(i,poly(0,-1)*a[j][col]*iv,j);
	}
	for(int i=1;i<=n;i++)ans=ans*a[i][i];
	return ans.a;
}
int main(){
	sieve();
	cin>>n>>m;
	for(int i=1;i<=m;i++)scanf("%d%d%d",fr+i,to+i,w+i),cnt[w[i]]++;
	for(int i=1;i<=152501;i++){
		int cnt0=0;
		for(int j=i;j<=152501;j+=i)cnt0+=cnt[j];
		if(cnt0<n-1)continue;
		memset(a,0,sizeof(a));
		for(int j=1;j<=m;j++)if(w[j]%i==0){
			int x=fr[j],y=to[j],v=w[j];
			a[x][x]+=poly(v,1),a[y][y]+=poly(v,1);
			a[x][y]-=poly(v,1),a[y][x]-=poly(v,1);
		}
		n--,det[i]=gauss(),n++;
//		cout<<i<<":"<<det[i]<<"!!\n";
	}
	int ans=0;
	for(int o=1;o<V;o++)for(int k=1;k*o<V;k++)(ans+=1ll*mu[o]*k*det[k*o]%mod)%=mod;
	cout<<(ans+mod)%mod;
	return 0;
}