「省选联考2020A卷」作业题

给定一张带边权的简单无向连通图，定义一棵生成树 $T$ 的权值为 $$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}})$$ 求所有生成树的权值之和。 $n\leqslant 30 , m\leqslant\frac{n(n-1)}{2},1\leqslant w_i\leqslant 152501$ 。

首先我们把碍事的 $\gcd$ 消掉，这里通过欧拉反演，即 $\varphi*I=id$ （证明的话把 $1\dots n$ 的数按与 $n$ 的 $\gcd$ 分类）。当然也可以直接减掉倍数好像还更方便。那么简单推下式子： $$\begin{aligned} &\sum\limits_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}})（现在开始设\gcd 为 g）\\ =&\sum\limits_T\left(\sum\limits_{i=1}^{n-1}w_{e_i}\right)\times\sum\limits_{d|g}\varphi(d)\\ =&\sum\limits_{d=1}^{w_{max}}\varphi(d)\sum\limits_{ T\wedge\forall i,d|w_{e_i}}\left(\sum\limits_{i=1}^{n-1}w_{e_i}\right) \end{aligned}$$ 我们枚举 $d$ ，每次找出所有 $w_i$ 为 $d$ 的倍数的边，我们现在需要求出仅由这些边构成的生成树的边权和的和。对于这个，我们构造一个模 $x^2$ 意义下的多项式矩阵，对每条边重新定义边权为 $1+w_ix$ 。基尔霍夫矩阵的主余子式等于 $\sum\limits_T\prod\limits_{i=1}^{n-1}w_i$ （这个我不会证明），那么我们把重新定义的边权带进去发现一次项系数就是我们要求的边权和之和。

然后我们考虑 $\bmod x^2$ 意义下的多项式四则运算，加减法没有区别，乘法：$(a+bx)(c+dx)=ac+(ad+bc)x$，除法：$\frac{a+bx}{c+dx}=\frac{a}{c}+\frac{bc-ad}{c^2}x$ （证明的话乘一个 $c+dx$ 就好了）。当 $c=0$ 怎么办？对于选取每一行的主元，我们找到不为0的次数最低的那一项，如果次数为0那么就直接做，如果为1说明这一列在当前行后面的常数项都为0，此时除法应该特判为 $\frac{bx}{dx}=\frac{b}{d}$（我不是很确定这里的必要性，代码是否特判都过了）。这样的话对 $1-w_{max}$ 都做一遍矩阵树复杂度是 $O(wn^3)$ 的，我们可以做一个剪枝，当添加入的边大于等于 $n-1$ 时才去做矩阵树，这样的复杂度是 $O(n^3\max\sigma_0w)$ 即 $144\times n^3$ ，可以通过本题。

代码：

#include<bits/stdc++.h>
using namespace std;
int const N=152501,p=998244353;
int n,u[2505],v[2505],w[2505],phi[N+5],pr[50005],cnt;
bool flag[N+5];
void init()
{
    phi[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!flag[i])pr[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&i*pr[j]<=N;j++)
        {
            flag[i*pr[j]]=1;
            if(!(i%pr[j])){phi[i*pr[j]]=phi[i]*pr[j];break;}
            phi[i*pr[j]]=phi[i]*(pr[j]-1);
        }
    }
}
inline int mod(const int &x){return x>=p?x-p:x;}
int pw(int x,int y)
{
    int res=1;
    while(y)
    {
        if(y&1)res=1ll*res*x%p;
        x=1ll*x*x%p;
        y>>=1;
    }
    return res;
}
struct poly
{
    int x0,x1;
    poly(){x0=x1=0;}
    poly(int a,int b){x0=a,x1=b;}
    friend poly operator +(poly a,poly b)
        {return poly(mod(a.x0+b.x0),mod(a.x1+b.x1));}
    friend poly operator -(poly a,poly b)
        {return poly(mod(a.x0-b.x0+p),mod(a.x1-b.x1+p));}
    friend poly operator *(poly a,poly b)
        {return poly(1ll*a.x0*b.x0%p,mod(1ll*a.x0*b.x1%p+1ll*a.x1*b.x0%p));}
    friend poly operator /(poly a,poly b)
    {
        if(a.x0==0&&b.x0==0)
        {
            int inv=pw(b.x1,p-2);
            return poly(1ll*a.x1*inv%p,0);
        }
        int inv=pw(b.x0,p-2);
        return poly(1ll*a.x0*inv%p,((1ll*a.x1*b.x0-1ll*a.x0*b.x1)%p*inv%p*inv%p+p)%p);
    }
}mat[31][31];
int deg(poly a){return a.x0?0:(a.x1?1:2);}
int det()
{
    int f=1;
    for(int i=1;i<=n;i++)
    {
        int c=i;
        for(int j=i+1;j<=n;j++)
            if(deg(mat[j][i])<deg(mat[c][i]))c=j;
        if(c!=i)swap(mat[i],mat[c]),f*=-1;
        if(deg(mat[i][i])==2)return 0;
        for(int j=i+1;j<=n;j++)
        {
            poly t=mat[j][i]/mat[i][i];
            for(int k=i;k<=n;k++)
                mat[j][k]=mat[j][k]-mat[i][k]*t;
        }
    }
    poly res;
    if(f==1)res=poly(1,0);
    else res=poly(p-1,0);
    for(int i=1;i<=n;i++)res=res*mat[i][i];
    return res.x1;
}
int main()
{
    init();
    int m,ans=0;
    scanf("%d%d",&n,&m);
    n--;
    for(int i=1;i<=m;i++)scanf("%d%d%d",&u[i],&v[i],&w[i]);
    for(int i=1;i<=N;i++)
    {
        memset(mat,0,sizeof(mat));
        int cnt=0;
        for(int j=1;j<=m;j++)
            if(!(w[j]%i))
            {
                mat[u[j]][v[j]]=mat[u[j]][v[j]]-poly(1,w[j]);
                mat[v[j]][u[j]]=mat[v[j]][u[j]]-poly(1,w[j]);
                mat[u[j]][u[j]]=mat[u[j]][u[j]]+poly(1,w[j]);
                mat[v[j]][v[j]]=mat[v[j]][v[j]]+poly(1,w[j]);
                cnt++;
            }
        if(cnt<n)continue;
        ans=mod(ans+1ll*det()*phi[i]%p);
    }
    printf("%d",ans);
    return 0;
}