正题

题目链接:https://www.luogu.com.cn/problem/P6624


题目大意

\(n\)个点的一张图,每条边有权值,一棵生成树的权值是所有边权和乘上边权的\(gcd\),即

\[val(T)=\left(\sum\limits_{i=1}^{n-1} w_{e_i}\right) \times \gcd(w_{e_1},w_{e_2},\dots,w_{e_{n-1}})
\]

求所有生成树的权值和


解题思路

首先要知道一个东西\(\varphi*I=id\),于是我们就有

\[gcd(w_1,w_2,\dots,w_n)=\sum_{d|w_1,d|w_2,\dots,d|w_n}\varphi(d)
\]

然后化进这个式子里就是

\[\left(\sum\limits_{i=1}^{n-1} w_{i}\right) \times \sum_{d|w_1,d|w_2,\dots,d|w_n}\varphi(d)
\]

然后把\(\varphi(d)\)丢出去就是

\[\sum_{d=1}^n\varphi(d)\times\left(\sum\limits_{i=1,d|w_i}^{n-1} w_{i}\right)
\]

之后就是怎么快速算后面那个东西的事情了。

一个比较朴素的做法是枚举每一条边产生的贡献,然后固定这条边然后矩阵树求一次答案乘上这条边就好了。

但是这个比较慢,我们可以利用生成函数的思想,每一条边权变为\(F_i(x)=w_ix+1\)的一个多项式,然后所有多项式乘起来之后的一次项系数就是答案了。因为这样只会选择一条边的\(w_i\)

多项式除法的话用的比较少,这里化一下

\[\frac{ax+b}{cx+d}=zx+w\Rightarrow ax+b=zcx^2+(zd+cw)x+wd
\]

然后二次项我们直接丢掉就有

\[b=wd\Rightarrow w=\frac{b}{d}
\]

带进后面那个去就有

\[a=zd+c\frac{b}{d}\Rightarrow z=\frac{ad-cb}{d^2}
\]

所以

\[\frac{ax+b}{cx+d}=\frac{ad-cb}{d^2}x+\frac{b}{d}
\]

这样时间复杂度是\(O(\ max\{w_i\}(m+n^3)\ )\)的,好像过不了,可以把边数不够\(n-1\)的直接不管,这样可以省去很多无用情况


code

#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const ll N=2e5+10,P=998244353;
struct fuc{
ll x,y;
fuc(ll xx=0,ll yy=0)
{x=xx;y=yy;return;}
};
ll power(ll x,ll b){
ll ans=1;x%=P;
while(b){
if(b&1)ans=ans*x%P;
x=x*x%P;b>>=1;
}
return ans;
};
ll inv(ll x){return power(x,P-2);}
fuc operator+(fuc a,fuc b)
{return fuc((a.x+b.x)%P,(a.y+b.y)%P);}
fuc operator-(fuc a,fuc b)
{return fuc((a.x-b.x+P)%P,(a.y-b.y+P)%P);}
fuc operator*(fuc a,fuc b)
{return fuc((a.x*b.y+a.y*b.x)%P,a.y*b.y%P);}
fuc operator/(fuc a,fuc b)
{return fuc((a.x*b.y-b.x*a.y+P)%P*inv(b.y*b.y)%P,a.y*inv(b.y)%P);}
ll n,m,ans,x[N],y[N],w[N],phi[N],pri[N],cnt;
bool v[N];fuc a[31][31];
fuc det(){
fuc ans(0,1);ll f=1;
for(ll i=1;i<n;i++){
ll w=i;
for(ll j=i;j<n;j++)
if(a[i][j].y){
if(i!=j)f=-f;
w=j;break;
}
swap(a[i],a[w]);
if(!a[i][i].y)return fuc(0,0);
ans=ans*a[i][i];
fuc inv=fuc(0,1)/a[i][i];
for(ll j=n-1;j>=i;j--)
a[i][j]=a[i][j]*inv;
for(ll j=i+1;j<n;j++){
fuc rate=a[j][i];
for(ll k=i;k<n;k++)
a[j][k]=a[j][k]-rate*a[i][k];
}
}
return (ans.x*f+P)%P;
}
void init(ll n){
phi[1]=1;
for(ll i=2;i<=n;i++){
if(!v[i])pri[++cnt]=i,phi[i]=i-1;
for(ll j=1;j<=cnt&&i*pri[j]<=n;j++){
v[i*pri[j]]=1;
if(i%pri[j]==0){
phi[i*pri[j]]=phi[i]*pri[j];
break;
}
phi[i*pri[j]]=phi[i]*phi[pri[j]];
}
}
return;
}
signed main()
{
scanf("%lld%lld",&n,&m);
ll mx=0;
for(ll i=1;i<=m;i++){
scanf("%lld%lld%lld",&x[i],&y[i],&w[i]);
x[i]--;y[i]--;mx=max(mx,w[i]);
}
init(mx);
for(ll i=1;i<=mx;i++){
ll cnt=0;
for(ll j=1;j<=m;j++)
if(w[j]%i==0)cnt++;
if(cnt<n-1)continue;
memset(a,0,sizeof(a));
for(ll j=1;j<=m;j++)
if(w[j]%i==0){
a[x[j]][x[j]]=a[x[j]][x[j]]+fuc(w[j],1);
a[y[j]][y[j]]=a[y[j]][y[j]]+fuc(w[j],1);
a[x[j]][y[j]]=a[x[j]][y[j]]-fuc(w[j],1);
a[y[j]][x[j]]=a[y[j]][x[j]]-fuc(w[j],1);
}
(ans+=phi[i]*det().x%P)%=P;
}
printf("%lld\n",ans);
return 0;
}

最新文章

  1. display:inline-block的间隙问题和解决办法
  2. SQL server存储过程语法及实例(转)
  3. Kibana源码剖析 —— savedSearch从读取到跳转
  4. 基础篇之 Create Type
  5. IOS开发关于测试的好的网址资源
  6. Recover Rotated Sorted Array
  7. [DevExpress]设置RepositoryItemComboBox只可下拉选择不可编辑
  8. dynamic 和var
  9. MsSql省市联动表
  10. CDZSC_2015寒假新人(1)——基础 a
  11. 批处理命令篇--配置免安装mysql 5.6.22, 以及1067错误的一个解决方法
  12. poj 3013 Big Christmas Tree
  13. hive导出查询文件到本地文件的2种办法
  14. maven的配置-2019-4-13
  15. VS2008安装“Visual Studio Web 创作组件”安装失败的解决方法
  16. XSS-HTML&amp;javaSkcript&amp;CSS&amp;jQuery&amp;ajax-CSS
  17. Caffe源码阅读(1) 全连接层
  18. 【UNIX网络编程(三)】TCP客户/server程序演示样例
  19. 小橙书阅读指南(十三)——连通性算法(union-find)
  20. SVN的基本操作

热门文章

  1. 多线程Synchronized的两种锁
  2. jwt《token》
  3. java.lang.NullPointerException at org.apache.jsp.index_jsp._jspInit(index_jsp.java:40)
  4. &#127942;【JVM技术专区】「难点-核心-遗漏」TLAB内存分配+锁的碰撞(技术串烧)!
  5. T-SQL - query02_查看数据库信息|查看服务器名称|查看实例名
  6. 1day漏洞反推技巧实战(2)
  7. vue 输入框禁止输入空格 ,只能输入数字,禁止输入数字
  8. Tars | 第3篇 Tars中期汇报测试文档(Java语言实现Subset路由规则)
  9. Git 系列教程(6)- 查看 commit 提交历史
  10. Asp.Net 5上传文件 (Core API方式)