（以下假设$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条边边权相同，那么最多有有$2\sqrt{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}$
[loj3304]作业题_数学-多项式[loj3304]作业题_数学-莫比乌斯反演_02
  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 }
View Code