Luogu6624 作业题 | Arahc's Home

题面

简化题意

给定一个带权无向图，求其所有生成树的权值和。一棵树的权值的定义为所有边权的和乘上所有边权的最大公约数。

$n\leq 30,w_i\leq 152501$ ，对 $998244353$ 取模。

题解

回顾矩阵树定理可以求什么：一个图所有生成树的边权积的和。

那么我们就要想办法把它变成这样的形式。

考虑欧拉反演 $\varphi\ast 1 = \mathrm{id}$ ：

\left(\sum_{e\in T}^{n-1} w_{e} \right) \times \sum_{ d\mid w_{ e_1 } , w_{ e_2 } , \ldots , w_{ e_{ m-1 }}} \varphi(d)

于是就可以考虑提出公约数（老生长谈了吧），记 $W = \max(w_i)$ ：

\sum_{d}^W \varphi(d) \sum_{ d\mid w_{ e_1 } , w_{ e_2 } , \ldots , w_{ e_{ m-1 }}} \left(\sum_{e\in T}^{n-1} w_{e} \right)

枚举 $d$ ，然后把是 $d$ 的倍数的边加入到图中，然后求这个图的所有生成树边权和的和即可。但这和前面预想的有点不太一样，因为矩阵树定理是用来求所有生成树边权积的和的。

~~如果你全部 exp 然后 ln 回来过了那你挺厉害的，我都不知道这是否可行~~。

那么怎么求生成树边权和的和？我们考虑把和转化为另外的东西的积（当然不是 $\exp,\ln$ ……），我们发现，对于每一个边权 $a$ ，将其写成一个一次函数 $1+ax$ ，那么两个一次函数的积：

(1+ax)(1+bx) = 1 + (a+b)x + abx^2

的一次项系数就是两个一次函数的乘积。

因此在模 $x^2$ 意义下处理对一次函数的四则运算即可，最难的估计是除法：

加： $(k_1x + b_1)+(k_2x + b_2) = (k_1 + k_2)x + (b_1 + b_2)$ 。
减：同加。
乘：上面有了。
除：考虑它的逆元是什么： $\frac{1}{kx+b} = -\frac{k}{b^2}x + \frac{1}{b}$ 。

于是就可以快乐爆枚+矩阵树定理了，注意直接做是 $\mathcal O(n^3W)$ 的，理论上不能通过。但是如果 $d$ 的倍数的边权个数少于 $n-1$ 是没有必要算矩阵的，复杂度为 $\mathcal O(n^3 \frac{\sum \sigma_0 (w_i)} {n-1})$ 。

这个式子看起来比较玄学，根据 OEIS 的序列 A002182， $w$ 全部取 $110880$ （有 $144$ 个约数）可以卡到一个不错的复杂度，可以变成 $n^4 \times 144$ ，在三秒的时限下，可以通过。

代码

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int max_n=33,max_m=200005,mod=998244353;
inline int read(){
    int x=0;bool w=0;char c=getchar();
    while(c<'0' || c>'9') w|=c=='-',c=getchar();
    while(c>='0' && c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return w?-x:x;
}
inline void write(int x){
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10^48);
}
inline int mi(int a,int p=mod-2){
    int res=1;
    while(p){
        if(p&1) res*=a,res%=mod;
        a*=a,a%=mod,p>>=1;
    }
    return res;
}

int pr[max_m],cntp,phi[max_m];
bool isp[max_m];
inline void Getpr(int n){
    memset(isp,1,sizeof(isp)),isp[1]=0;
    phi[1]=1;
    for(register int i=2;i<=n;++i){
        if(isp[i]){
            pr[++cntp]=i,
            phi[i]=i-1;
        }
        for(register int j=1;j<=cntp && i*pr[j]<=n;++j){
            isp[i*pr[j]]=0;
            if(i%pr[j]==0){
                phi[i*pr[j]]=phi[i]*pr[j];
                break;
            }
            phi[i*pr[j]]=phi[i]*phi[pr[j]];
        }
    }
}

struct line{
    int k,b;
    line operator + (const line &b) const{
        line a=*this,res;
        res.k=(a.k+b.k)%mod,
        res.b=(a.b+b.b)%mod;
        return res;
    }
    line operator - (const line &b) const{
        line a=*this,res;
        res.k=(a.k-b.k+mod)%mod,
        res.b=(a.b-b.b+mod)%mod;
        return res;
    }
    line operator * (const line &b) const{
        line a=*this,res;
        res.k=(a.k*b.b+a.b*b.k)%mod,
        res.b=a.b*b.b%mod;
        return res;
    }
    line operator / (const line &b) const{
        line a=*this,res;
        int ivb=mi(b.b);
        res.k=(a.k*b.b%mod-a.b*b.k%mod+mod)%mod*ivb%mod*ivb%mod,
        res.b=a.b*ivb%mod;
        return res;
    }
};
inline line makel(int k,int b){
    line res;
    res.k=k,res.b=b;
    return res;
}

int n,m,W,ans;

struct matrix{
    line a[max_n][max_n];
    inline void clear(){
        for(register int i=0;i<=n;++i)
            for(register int j=0;j<=n;++j)
                a[i][j]=makel(0,0);
    }
    inline void add(int x,int y,line w){
        a[x][x]=a[x][x]+w,
        a[y][y]=a[y][y]+w,
        a[x][y]=a[x][y]-w,
        a[y][x]=a[y][x]-w;
    }
    inline int calc(){
        line res=makel(0,1);
        for(register int i=1;i<n;++i){
            int pos=i;
            while(pos<n && !a[pos][i].b)
                ++pos;
            if(pos!=i){
                swap(a[pos],a[i]);
                res=makel(0,0)-res;
            }
            for(register int j=i+1;j<n;++j){
                line tmp=a[j][i]/a[i][i];
                for(register int k=i;k<n;++k)
                    a[j][k]=a[j][k]-(a[i][k]*tmp);
            }
        }
        for(register int i=1;i<n;++i)
            res=res*a[i][i];
        return (res.k%mod+mod)%mod;
    }
}a;

struct edge{
    int u,v,w;
}e[max_m];
vector<edge> ed[max_m];

signed main(){
    n=read(),m=read();
    for(register int i=1;i<=m;++i){
        e[i].u=read(),e[i].v=read(),e[i].w=read();
        W=max(W,e[i].w);
    }
    Getpr(W);
    for(register int d=1;d<=W;++d)
        for(register int i=1;i<=m;++i)
            if(e[i].w%d==0)
                ed[d].push_back(e[i]);
    for(register int d=1;d<=W;++d){
        int _n=ed[d].size();
        if(_n<n-1) continue;
        a.clear();
        for(register int i=0;i<_n;++i)
            a.add(ed[d][i].u,ed[d][i].v,makel(ed[d][i].w,1));
        ans=(ans+phi[d]*a.calc())%mod;
    }
    write((ans%mod+mod)%mod);
    return 0;
}