题解 P6624 【[省选联考 2020 A 卷] 作业题】

· · 题解

更好的阅读体验

题目链接

前置知识：(1) 矩阵树定理。(2) 多项式四则运算（加、减、乘、求逆（牛顿迭代））。

网上的介绍很多。这里就不细讲了。

首先要把这个\gcd给拆掉，才能用我们熟悉的“矩阵树定理”等一系列方法来解题。

众所周知，\sum_{d|x}\varphi(d)=x。本题中，我们就可以利用这个\varphi来拆掉\gcd。具体来说，如果我们用T来表示枚举的一棵树，e_1,e_2,\dots ,e_{n-1}来表示T包含的边的编号，那么首先，答案可以表示为：

\sum_{T}\left(\sum_{i=1}^{n-1}w_{e_i}\right)\times\gcd(w_{e_1},w_{e_2},\dots,w_{e_{n-1}})

我们把\gcd(w_{e_1},w_{e_2},\dots,w_{e_{n-1}})用上述的公式带入，得到：

=\sum_{T}\left(\sum_{i=1}^{n-1}w_{e_i}\right)\times\left(\sum_{d|\gcd(w_{e_1},w_{e_2},\dots,w_{e_{n-1}})}\varphi(d)\right)

然后就可以把这个d，放到前面来，得到：

=\sum_{d=1}^{W}\varphi(d)\times\left(\sum_{\begin{gather*}T\\d|w_{e_1},\dots ,d|w_{e_{n-1}}\end{gather*}}\sum_{i=1}^{n-1}w_{e_i}\right)

其中，W表示最大的权值。后面括号里的部分，相当于【所有【w_{e_i}是d的倍数的边】组成的子图】的所有生成树的边权和。

普通的矩阵树定理，可以用来求所有生成树的边权积之和。其中最为大家所熟知的应用，就是当所有边权都为1时，就相当于是求生成树的数量，也就是生成树计数问题。它的实现，就是求矩阵行列式，因为需要用到高斯消元，所以时间复杂度是O(n^3)的。

但是本题中，要求的是边权和，而不是“边权积之和”。一种朴素的想法是考虑每一条边的贡献，也就是求【包含这条边的生成树数量】。它又等于【原图的生成树数量】减去【原图去掉这条边后的生成树数量】。对每条边都求一次【原图去掉这条边后的生成树数量】，时间复杂度O(mn^3)。因为外层还要枚举d，所以总时间复杂度O(W\cdot mn^3)，无法通过本题。

求所有生成树的边权和，其实有更好的方法。我们把每条边的边权，看做一个多项式：(1+w_{e_i}x)。其中，x不是任何具体的数，它只是一个符号，表示多项式的“一次项”。那么，一个生成树，它的边权之积的“一次项”前的系数，就是这颗生成树的边权和（这可以根据多项式乘法的定义来理解）。

如此以来，求所有生成树的边权和的问题，当我们把边权换成这样一个多项式后，就转化为了求所有生成树的边权积之和的问题。只不过，现在新的边权，不再是一个数，而是一个多项式。所以我们只需要定义出多项式的四则运算，就可以直接用矩阵树定理求解了。另外，我们只关心多项式的“一次项”系数，所以更高的项可以舍去。换句话说，所有的多项式运算，可以在\bmod x^2意义下进行。那么把多项式定义为一个\texttt{pair}就可以了。

多项式的加、减、乘运算都比较简单。不了解的可以看如下代码：

typedef pair<int,int> pii;
#define fi first
#define se second
//为了便于理解，这里就不取模了
pii operator+(const pii& a,const pii& b){
    return mk(a.fi+b.fi,a.se+b.se);
}
pii operator-(const pii& a,const pii& b){
    return mk(a.fi-b.fi,a.se-b.se);
}
pii operator*(const pii& a,const pii& b){
    return mk(a.fi*b.fi,a.se*b.fi+b.se*a.fi);
}

多项式的除法，因为在\bmod x^2意义下进行，所以只需要求出多项式\bmod x^2意义下的逆即可。用经典的倍增法。设原多项式为A(x)，要求的、它的逆为F(x)=A^{-1}(x)\pmod{x^n}。如果已经求出了F_0(x)=A^{-1}(x)\pmod{x^{\lceil\frac{n}{2}\rceil}}，那么：

F(x)=2F_0(x)-F_0^2(x)A(x)\pmod{x^{n}}

对本题来说，我们就不需要递归了。先求出常数项，也就是求逆元。然后把常数项作为F_0(x)带入上式即可。具体的代码：

pii inv(const pii& a){
    int t=pow_mod(a.fi,MOD-2);
    return mk(2*t,0)-mk(t,0)*mk(t,0)*a;
}
pii operator/(const pii& a,const pii& b){
    return a*inv(b);
}

还有一个小细节是，高斯消元求行列式的时候，如果当前行的（要作为被除数）的这个多项式，常数项为0，上述的求逆的方法就不管用了，因为0没有逆元。这种情况下，我们首先考虑，在它下面，找一个常数项不为0的行，与它交换（根据行列式的性质，两行交换后行列式要变号，也就是res:=-res）。如果它下面所有行常数项都为0，那对这一列，我们就不需要管常数项了，直接拿一次项消就行（就把每一项都看成普通的数而不是多项式就行）。

至此，所有四则运算都可以O(1)实现，所以做一次高斯消元，求出行列式的复杂度就是O(n^3)的。外层还要枚举d，所以总时间复杂度O(Wn^3)。这个复杂度仍然无法通过本题，我们还需要再对W这部分做一些优化。

发现对于一个d。如果“是它的倍数”的边，数量小于n-1条，那必不可能有任何生成树。特判这一情况后，我们惊喜地发现，复杂度降为：O(\frac{\sum\sigma_0(w_{e_i})}{n-1}\cdot n^3)=O(n^2\sum\sigma_0(w_{e_i}))。其中\sigma_0(x)表示x的约数数量。这个复杂度很好理解，因为每个d，要作为约数出现n-1次才被计算，所以被计算到的d最多只有\frac{\sum\sigma_0(w_{e_i})}{n-1}个。再看这个复杂度，一个小于等于W的数，约数个数是O(\sqrt{W})级别的，又因为边数是O(n^2)级别的，所以总复杂度就是O(n^4\sqrt{W})。事实上，这个\sqrt{W}只是一个理论上限，在本题的数据范围下，打表发现，w的约数个数最多为144。足以通过本题。

参考代码（在LOJ查看）：

//problem:LOJ3304
#include <bits/stdc++.h>
using namespace std;

#define pb push_back
#define mk make_pair
#define lob lower_bound
#define upb upper_bound
#define fi first
#define se second
#define SZ(x) ((int)(x).size())

typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;

const int MOD=998244353;
inline int mod1(int x){return x<MOD?x:x-MOD;}
inline int mod2(int x){return x<0?x+MOD:x;}
inline void add(int& x,int y){x=mod1(x+y);}
inline void sub(int& x,int y){x=mod2(x-y);}
inline int pow_mod(int x,int i){int y=1;while(i){if(i&1)y=(ll)y*x%MOD;x=(ll)x*x%MOD;i>>=1;}return y;}

pii operator+(const pii& lhs,const pii& rhs){
    return mk(mod1(lhs.fi+rhs.fi),mod1(lhs.se+rhs.se));
}
pii operator-(const pii& lhs,const pii& rhs){
    return mk(mod2(lhs.fi-rhs.fi),mod2(lhs.se-rhs.se));
}
pii operator*(const pii& lhs,const pii& rhs){
    return mk((ll)lhs.fi*rhs.fi%MOD,((ll)lhs.se*rhs.fi+(ll)rhs.se*lhs.fi)%MOD);
}
pii inv(const pii& a){
    //assert(a.fi!=0);
    int t=pow_mod(a.fi,MOD-2);
    return mk(2LL*t%MOD,0)-mk(t,0)*mk(t,0)*a;//牛顿迭代
}
pii operator/(const pii& lhs,const pii& rhs){
    return lhs*inv(rhs);
}
pii operator-(const pii& rhs){
    return mk(mod2(-rhs.fi),mod2(-rhs.se));
}//负号
pii& operator+=(pii& lhs,const pii& rhs){
    lhs=lhs+rhs;
    return lhs;
}
pii& operator-=(pii& lhs,const pii& rhs){
    lhs=lhs-rhs;
    return lhs;
}
pii& operator*=(pii& lhs,const pii& rhs){
    lhs=lhs*rhs;
    return lhs;
}
pii& operator/=(pii& lhs,const pii& rhs){
    lhs=lhs/rhs;
    return lhs;
}

const int MAXW=152501;
int p[MAXW+5],cnt,phi[MAXW+5];
bool v[MAXW+5];
void sieve(int lim){
    phi[1]=1;
    for(int i=2;i<=lim;++i){
        if(!v[i]){
            phi[i]=i-1;
            p[++cnt]=i;
        }
        for(int j=1;j<=cnt && p[j]*i<=lim;++j){
            v[i*p[j]]=1;
            if(i%p[j]==0){
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
            phi[i*p[j]]=phi[i]*phi[p[j]];
        }
    }
}//线性筛,求phi

const int MAXN=30;
const int MAXM=MAXN*(MAXN-1)/2;
int n,m;
struct Edge{
    int u,v,w;
}e[MAXM+5];

struct Matrix{
    int n;
    pii a[MAXN+5][MAXN+5];
    Matrix(){
        n=0;
        memset(a,0,sizeof(a));
    }
};

pii get_det2(Matrix A,int i){
    //常数项全部为0
    pii res=mk(1,0);
    int line=0;
    for(int j=i;j<=A.n;++j){
        if(A.a[j][i].se!=0){
            line=j;
            break;
        }
    }
    if(!line)return mk(0,0);
    if(line!=i){
        for(int j=1;j<=A.n;++j)swap(A.a[i][j],A.a[line][j]);
        res=-res;
    }
    int _inv=pow_mod(A.a[i][i].se,MOD-2);
    for(int j=i+1;j<=A.n;++j){
        int t=(ll)A.a[j][i].se*_inv%MOD;
        for(int k=1;k<=A.n;++k)
            sub(A.a[j][k].se,(ll)A.a[i][k].se*t%MOD);
    }
    res*=A.a[i][i];
    return res;
}
pii get_det(Matrix A){
    pii res=mk(1,0);
    for(int i=1;i<=A.n;++i){
        int line=0;
        for(int j=i;j<=A.n;++j){
            if(A.a[j][i].fi!=0){
                line=j;
                break;
            }
        }
        if(line==0){
            res*=get_det2(A,i);
            continue;
        }
        if(line!=i){
            for(int j=1;j<=A.n;++j)swap(A.a[i][j],A.a[line][j]);
            res=-res;
        }
        pii _inv=inv(A.a[i][i]);
        for(int j=i+1;j<=A.n;++j){
            pii t=A.a[j][i]*_inv;
            for(int k=1;k<=A.n;++k)
                A.a[j][k]-=A.a[i][k]*t;
        }
        res*=A.a[i][i];
    }
    return res;
}

int main() {
    cin>>n>>m;
    int max_w=0;
    for(int i=1;i<=m;++i){
        cin>>e[i].u>>e[i].v>>e[i].w;
        max_w=max(max_w,e[i].w);
    }
    sieve(max_w);
    int ans=0;
    for(int i=1;i<=max_w;++i){
        int cnt_e=0;
        for(int j=1;j<=m;++j)if(e[j].w%i==0)cnt_e++;
        if(cnt_e<n-1)continue;

        Matrix mat;
        for(int j=1;j<=m;++j)if(e[j].w%i==0){
            mat.a[e[j].u][e[j].u]+=mk(1,e[j].w);
            mat.a[e[j].v][e[j].v]+=mk(1,e[j].w);
            mat.a[e[j].u][e[j].v]-=mk(1,e[j].w);
            mat.a[e[j].v][e[j].u]-=mk(1,e[j].w);
        }
        mat.n=n-1;
        pii det=get_det(mat);
        ans=((ll)ans+(ll)phi[i]*det.se)%MOD;
    }
    cout<<ans<<endl;
    return 0;
}