[矩阵树][生成函数][欧拉反演]luogu P6624 [省选联考 2020 A 卷] 作业题

题面

https://www.luogu.com.cn/problem/P6624

分析

后面那坨 gcd 直接欧拉反演掉，然后枚举 d ，原式形式变成

$\sum_{d=1}^{w_{max}} \sum_{T,d|w_{e1},d|w_{e2},...,d|w_{en-1}}\sum_{i=1}^{n}e_i$

后面的东西显然是矩阵树定理，但是我们只学过求边权积，没学过求边权和

因为矩阵树是乘积转移，考虑引入生成函数的思想

初始矩阵每条边为 $f(ei)=1+w_{ei}x$

那么考虑相乘完的结果， n 次项表示选择了 n 个边的边权和，所以最终边权积只用取一次项

又因为没有负数项，相乘时不会降次，所以运算过程中只用保留前两项即可，重载一下运算符便于计算

时间复杂度 $O(w_{max}n^3)$ 会炸，那只对于至少能选出 n-1 条边的 d 做即可，时间复杂度约为 $O(w_{max}n^2)$ 跑不满

代码

//矩阵树求所有生成树边权和 
//莫比乌斯反演 
#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
typedef long long ll;
const ll P=998244353;
const int N=31;
const int W=152502;
int n,m,mx,u[N*N],v[N*N],w[N*N],phi[W],pri[W],num[W];
bool unc[W];
ll ans;

ll Pow(ll x,ll y) {ll ans=1;for (;y;y>>=1,x=x*x%P) if (y&1) ans=ans*x%P;return ans;}
struct GF {
    ll a0,a1;
    friend GF operator + (GF a,GF b) {return (GF){(a.a0+b.a0)%P,(a.a1+b.a1)%P};}
    friend GF operator - (GF a,GF b) {return (GF){(a.a0-b.a0+P)%P,(a.a1-b.a1)%P};}
    friend GF operator * (GF a,GF b) {return (GF){a.a0*b.a0%P,(a.a0*b.a1%P+a.a1*b.a0%P)%P};}
    friend GF operator / (GF a,GF b) {return (GF){a.a0*Pow(b.a0,P-2)%P,(a.a1*b.a0%P-a.a0*b.a1+P)%P*Pow(b.a0*b.a0%P,P-2)%P};}
}g[N][N];

void Prime() {
    phi[1]=1;
    for (int i=2;i<=mx;i++) {
        if (!unc[i]) pri[++pri[0]]=i,phi[i]=i-1;
        for (int j=1;j<=pri[0]&&1ll*i*pri[j]<=mx;j++) {
            unc[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]];
        }
    }
}

void Add(int u,int v,int w) {
    GF e=(GF){1,w};
    g[u][u]=g[u][u]+e;g[v][v]=g[v][v]+e;g[u][v]=g[u][v]-e;g[v][u]=g[v][u]-e;
}

ll Gauss() {
    int c=1;
    GF ans=(GF){1,0};
    for (int i=1;i<n;i++) {
        if (!g[i][i].a0)
            for (int j=i+1;j<n;j++)
                if (g[j][i].a0) {swap(g[i],g[j]),c*=-1;break;}
        for (int j=i+1;j<n;j++) {
            GF k=g[j][i]/g[i][i];
            for (int l=i;l<n;l++) g[j][l]=g[j][l]-g[i][l]*k;
        }
        ans=ans*g[i][i];
    }
    return ans.a1*c;
}

int main() {
    scanf("%d%d",&n,&m);
    for (int i=1;i<=m;i++) {
        scanf("%d%d%d",&u[i],&v[i],&w[i]);mx=max(mx,w[i]);
        for (int j=1;j*j<=w[i];j++)
            if (w[i]%j==0) num[j]++,num[w[i]/j]+=(j*j!=w[i]);
    }
    Prime();
    for (int i=1;i<=mx;i++)
        if (num[i]>=n-1) {
            memset(g,0,sizeof g);
            for (int j=1;j<=m;j++)
                if (w[j]%i==0) Add(u[j],v[j],w[j]);
            (ans+=Gauss()*phi[i]%P)%=P;
        }
    printf("%lld\n",(ans+P)%P);
}