（以下假设$T=(V,{e_{1},e_{2},...,e_{n-1} })$是一棵树）
根据莫比乌斯反演，有$gcd(w_{1},w_{2},...,w_{e_{n-1}})=sum_{d|w_{e_{i}}}varphi(d)$
容易想到枚举$d$，之后相当于求$sum_{d}varphi(d)sum_{d|w_{e_{i}}}sum_{i=1}^{n-1} w_{e_{i}}$
那么即找出所有权值为$d$的倍数的边，然后求所有生成树的边权和的和，但矩阵树定理只能求边权积之和（将邻接矩阵权值改为边权，将度数矩阵的度数改为出边的边权和，即【bzoj3534】）
考虑令边i的边权为$1+w_{i}x$，那么容易发现所有边权乘积的和的1次项系数即为答案
同时多项式只需要维护0次项和1次项（即模$x^{2}$意义下运算），因此复杂度为$o(dn^{3})$，仍然无法通过
考虑一个优化：如果$d$的倍数的边小于$n-1$条，显然无解；那么最坏的构造方式为：选出n组边，每组n条边边权相同，那么最多有有$2sqrt{d}n$种合法的d，复杂度降为$o(sqrt{d}n^{4})$，且其中会产生许多重复边权，并不会跑满
还有一件事情，这样的高斯消元会有很多细节：
1.多项式除法，由于在模$x^{2}$意义下，所以直接手动解方程推出$(ax+b)^{-1}equiv frac{1}{b}-frac{a}{b^{2}}x(mod x^{2})$
2.若存在常数项非0，将该点与当前行交换；若不存在，那么直接用一次项系数消除（因为常数项都为0），即$inv(ax)=a^{-1}$
View Code
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define N 35 4 #define D 152501 5 #define mod 998244353 6 #define mp make_pair 7 #define pii pair<int,int> 8 #define fi first 9 #define se second 10 struct ji{ 11 int nex,to,len; 12 }edge[N*N]; 13 int E,n,m,x,y,z,ans,head[N],p[D+10],phi[D+10]; 14 pii d[N][N]; 15 int ksm(int n,int m){ 16 if (!m)return 1; 17 int s=ksm(n,m>>1); 18 s=1LL*s*s%mod; 19 if (m&1)s=1LL*s*n%mod; 20 return s; 21 } 22 pii add(pii x,pii y){ 23 return mp((x.fi+y.fi)%mod,(x.se+y.se)%mod); 24 } 25 pii mul(pii x,pii y){ 26 return mp(1LL*x.fi*y.fi%mod,(1LL*x.fi*y.se+1LL*x.se*y.fi)%mod); 27 } 28 pii inv(pii x){ 29 if (!x.fi)return mp(ksm(x.se,mod-2),0); 30 int s=ksm(x.fi,mod-2); 31 return mp(s,mod-1LL*x.se*s%mod*s%mod); 32 } 33 void add(int x,int y,int z){ 34 edge[E].nex=head[x]; 35 edge[E].to=y; 36 edge[E].len=z; 37 head[x]=E++; 38 } 39 int find1(int k){ 40 for(int i=k;i<n;i++) 41 if (d[i][k].fi)return i; 42 return 0; 43 } 44 int find2(int k){ 45 for(int i=k;i<n;i++) 46 if (d[i][k].se)return i; 47 return 0; 48 } 49 int main(){ 50 phi[1]=1; 51 for(int i=2;i<=D;i++){ 52 if (!phi[i]){ 53 p[++p[0]]=i; 54 phi[i]=i-1; 55 } 56 for(int j=1;(j<=p[0])&&(i*p[j]<=D);j++){ 57 if (i%p[j])phi[i*p[j]]=phi[i]*phi[p[j]]; 58 else{ 59 phi[i*p[j]]=phi[i]*p[j]; 60 break; 61 } 62 } 63 } 64 scanf("%d%d",&n,&m); 65 memset(head,-1,sizeof(head)); 66 for(int i=1;i<=m;i++){ 67 scanf("%d%d%d",&x,&y,&z); 68 add(x,y,z); 69 add(y,x,z); 70 } 71 for(int i=1;i<=D;i++){ 72 x=0; 73 for(int j=0;j<E;j+=2) 74 if (edge[j].len%i==0)x++; 75 if (x<n-1)continue; 76 memset(d,0,sizeof(d)); 77 for(int j=1;j<n;j++) 78 for(int k=head[j];k!=-1;k=edge[k].nex) 79 if (edge[k].len%i==0){ 80 d[j][j]=add(d[j][j],mp(1,edge[k].len)); 81 if (edge[k].to<n)d[j][edge[k].to]=mp(mod-1,mod-edge[k].len); 82 } 83 pii s=mp(1,0); 84 for(int j=1;j<n;j++){ 85 x=find1(j); 86 if (!x)y=find2(j); 87 if (!x){ 88 s.se=0; 89 break; 90 } 91 if (x!=j){ 92 s=mp(mod-s.fi,mod-s.se); 93 for(int k=j;k<n;k++)swap(d[j][k],d[x][k]); 94 } 95 s=mul(s,d[j][j]); 96 pii t=inv(d[j][j]); 97 for(int k=j;k<n;k++)d[j][k]=mul(d[j][k],t); 98 for(int k=j+1;k<n;k++){ 99 t=mp(mod-d[k][j].fi,mod-d[k][j].se); 100 for(int l=j;l<n;l++)d[k][l]=add(d[k][l],mul(t,d[j][l])); 101 } 102 } 103 ans=(ans+1LL*s.se*phi[i])%mod; 104 } 105 printf("%d",ans); 106 }