【省选2020A卷】作业题【矩阵树】【扩域】【莫比乌斯反演】

传送门

为什么世界上会有这么傻的题啊……我佛了

很显然就是矩阵树强行缝合莫反

f ( n ) f(n) f(n)表示所有边权 gcd ⁡ \gcd gcd n n n的生成树权值和, g ( d ) g(d) g(d)表示所有边权都是 d d d的倍数的生成树权值和

g ( d ) = ∑ d ∣ n f ( n ) g(d)=\sum_{d\mid n}f(n) g(d)=dnf(n)

f ( d ) = ∑ d ∣ n g ( n ) μ ( n d ) f(d)=\sum_{d\mid n}g(n)\mu(\frac nd) f(d)=dng(n)μ(dn)

做完了

开权值大小个vector,把每条边压到边权的约数的vector中

如果一个vector中边数小于 n − 1 n-1 n1就跳过,否则对每条边构造 ( w x + 1 ) (wx+1) (wx+1)算矩阵树得到 g g g,然后无脑算出 f f f统计答案就好了

显然矩阵树次数不会超过 O ( m d ( V ) n ) = O ( n d ( V ) ) O(\frac {md(V)}n)=O(nd(V)) O(nmd(V))=O(nd(V)),总复杂度 O ( n 4 d ( v ) ) O(n^4d(v)) O(n4d(v))

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <vector>
using namespace std;
struct edge{int u,v,w;};
typedef long long ll;
const int MOD=998244353,N=152501;
inline int qpow(int a,int p)
{
	int ans=1;
	while (p)
	{
		if (p&1) ans=(ll)ans*a%MOD;
		a=(ll)a*a%MOD;p>>=1;
	}
	return ans;
}
#define inv(x) qpow(x,MOD-2)
int np[N+5],pl[N+5],mu[N+5],cnt;
inline void init()
{
	np[1]=mu[1]=1;
	for (int i=2;i<=N;i++)
	{
		if (!np[i]) mu[pl[++cnt]=i]=-1;
		for (int j=1,x;(x=i*pl[j])<=N;j++)
		{
			np[x]=1;
			if (i%pl[j]==0) break;
			mu[x]=(MOD-mu[i])%MOD;
		}
	}
}
inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;}
inline int dec(const int& x,const int& y){return x<y? x-y+MOD:x-y;}
struct poly{int x,y;poly(const int& x=0,const int& y=0):x(x),y(y){}};
inline poly operator +(const poly& a,const poly& b){return poly(add(a.x,b.x),add(a.y,b.y));}
inline poly operator -(const poly& a,const poly& b){return poly(dec(a.x,b.x),dec(a.y,b.y));}
inline poly operator *(const poly& a,const poly& b){return poly(((ll)a.x*b.y+(ll)a.y*b.x)%MOD,(ll)a.y*b.y%MOD);}
inline poly operator /(const poly& a,const poly& b){ll t=inv(b.y);return poly(dec((ll)a.x*b.y%MOD,(ll)a.y*b.x%MOD)*t%MOD*t%MOD,a.y*t%MOD);}
int n,mx;
poly g[35][35];
inline int det()
{
	poly ans(0,1);
	for (int i=1;i<n;i++)
	{
		ans=ans*g[i][i];
		for (int j=i+1;j<n;j++)
		{
			poly t=g[j][i]/g[i][i];
			for (int k=i;k<n;k++)
				g[j][k]=g[j][k]-g[i][k]*t;
		}		
	}
	return ans.x;
}
vector<edge> lis[N+5];
int F[N+5],G[N+5];
int main()
{
	init();
	int m;
	scanf("%d%d",&n,&m);
	for (int i=1;i<=m;i++)
	{
		edge e;
		scanf("%d%d%d",&e.u,&e.v,&e.w);
		mx=max(mx,e.w);
		for (int i=1;i*i<=e.w;i++)
			if (e.w%i==0)
			{
				lis[i].push_back(e);
				if (i*i<e.w) lis[e.w/i].push_back(e);
			}
	}
	for (int k=1;k<=mx;k++)
	{
		if ((int)lis[k].size()<n-1) continue;
		for (int i=1;i<=n;i++)
			for (int j=1;j<=n;j++)
				g[i][j]=poly(0,0);
		for (vector<edge>::iterator it=lis[k].begin();it!=lis[k].end();++it)
		{
			poly t(it->w,1);
			int u=it->u,v=it->v;
			g[u][u]=g[u][u]+t,g[v][v]=g[v][v]+t;
			g[u][v]=g[u][v]-t,g[v][u]=g[v][u]-t;	
		}	
		G[k]=det();
	}
	for (int d=1;d<=mx;d++)
		for (int i=1;i*d<=mx;i++)
			F[d]=(F[d]+(ll)G[i*d]*mu[i])%MOD;
	int ans=0;
	for (int i=1;i<=mx;i++) ans=(ans+(ll)i*F[i])%MOD;
	printf("%d\n",ans);
	return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值