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

传送门

好题啊，真的好题啊！

莫比乌斯反演+矩阵树定理（甚至还算加了一点点的生成函数的原理）

---

【分析】

根据题意，求对所有的树 (T) 的值：

(displaystyle ans=sum_{T}val(T)=sum_{T}(sum_{i=1}^{n-1}w_{e_i})gcd_{i=1}^{n-1} w_{e_i})

(e_i) 表示树 (T) 的边

由欧拉反演得：(displaystyle n=sum_{dmid n}oldsymbol varphi(d))

因此 (displaystyle gcd_{i=1}^{n-1} w_{e_i}=sum_{dmid gcd_{i=1}^{n-1}w_{e_i}}oldsymbol varphi(d)=sum_{forall i o dmid w_{e_i}}oldsymbol varphi(d))

因此 (displaystyle ans=sum_T(sum_{i=1}^{n-1}w_{e_i})sum_{forall i o dmid w_{e_i}}oldsymbol varphi(d))

调换枚举顺序得 (displaystyle ans=sum_{dgeq 1}oldsymbol varphi(d)sum_T[forall i o dmid w_{e_i}](sum_{i=1}^{n-1}w_{e_i}))

于是对于每个 (w_{e_i}) ，我们将它塞进它各因数的 vector 中。之后，每个数 (d) 的 vector 中，存着所有边权为它倍数的边

现考虑给定一个 vector ，里面存着若干边，用这些边生成树的权值和为多少，即 (displaystyle f(d)=sum_{Tin v_d}(sum_{i=1}^{n-1}w_{e_i}))

(ecause displaystyle sum_{i=1}^{n-1}w_{e_i}=sum_{i=1}^{n-1}w_{e_i}cdot 1)

若令 (P_i(x)=1+w_{e_i}x)

则 (displaystyle sum_{i=1}^{n-1}w_{e_i}cdot 1=sum_{c_1+c_2+c_3+cdots+c_{n-1}=1}prod_{i=1}^{n-1}P_i[c_i])

后面这玩意儿本质上就是卷积，而且还是 (1) 次项的系数

( herefore displaystyle sum_{i=1}^{n-1}w_{e_i}=(prod_{i=1}^{n-1}P_i(x))[1])

( herefore displaystyle f(d)=sum_{Tin v_d}((prod_{i=1}^{n-1}P_i(x))[1])=(sum_{Tin v_d}prod_{i=1}^{n-1}P_i(x))[1])

于是这就构成了矩阵树的形式。我们将 vector 中的边，构造成对应的一次函数，放进矩阵中，求一下行列式在模 (x^2) 意义下的值。再把一次项系数返回，就是 (f(d)) 的结果

消元方法见下

这样一来，(displaystyle ans=sum_{dgeq 1}oldsymbol varphi(d)f(d))

如果每个 (d) 都求解，复杂度将达到 (O(mcdot n^3)) ，其中 (m) 为边的最大值。复杂度直接爆炸。

考虑将 vector 中的边跑一下并查集。如果最后没有并成一个连通分量，那一定不会生成树，贡献为 (0) ，直接不跑高斯消元。

这样的剪枝相当于额外增加一个 (O(mncdot alpha(n))) ，但感觉是能扣除一大部分高斯消元。试试能不能过（最后当然是过了）。

---

简单讲一下多项式怎么在模 (x^2) 意义下高斯消元。将 (P(x)=a+bx) 记为 ((a,b)) 则：

((a,b)pm(c,d)=(apm c,bpm d))

((a+bx)(c+dx)=ac+(ad+bc)x+bdx^2equiv ac+(ad+bc)xpmod {x^2})

即 ((a,b)cdot (c,d)=(ac,ad+bc))

若 ((a+bx)(c+dx)equiv 1pmod {x^2})

则 (ac=1,ad+bc=0)

故 (c=a^{-1},d=-(bc)cdot a^{-1}=-bc^2)

因此 ((a,b)^{-1}=(a^{-1},-a^{-2}b))

故 ((a,b)/(c,d)=(a,b)cdot (c,d)^{-1})

这样一来就可以重载运算符了：

typedef pair<int,int> pii;
#define fi first
#define se second
inline pii operator ~ (const pii &p) { ll tmp=inv(p.fi); return pii(tmp,dis(0,tmp*tmp%MOD*p.se%MOD)); }
//逆元
inline pii operator - (const pii &p) { return pii( dis(0,p.fi) , dis(0,p.se) ); }
inline pii operator + (const pii &a,const pii &b) { return pii( add(a.fi,b.fi) , add(a.se,b.se) ); }
inline pii operator += (pii &a,const pii &b) { return a=a+b; }
inline pii operator - (const pii &a,const pii &b) { return pii( dis(a.fi,b.fi) , dis(a.se,b.se) ); }
inline pii operator -= (pii &a,const pii &b) { return a=a-b; }
inline pii operator * (const pii &a,const pii &b) { return pii( (ll)a.fi*b.fi%MOD , ((ll)a.fi*b.se+(ll)a.se*b.fi)%MOD ); }
inline pii operator *= (pii &a,const pii &b) { return a=a*b; }
inline pii operator / (const pii &a,const pii &b) { return a*(~b); }
inline pii operator /= (pii &a,const pii &b) { return a=a/b; }

---

【代码】

最后奉上本蒟蒻的代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<int,int> pii;
#define fi first
#define se second
const int MOD=998244353,MAXN=32,Lim=152501;
inline ll fpow(ll a,ll x) { ll ans=1; for(;x;x>>=1,a=a*a%MOD) if(x&1) ans=ans*a%MOD; return ans; }
//快速幂
inline ll inv(ll a) { return fpow(a,MOD-2); }
inline int add(int a,int b) { return (a+b>=MOD)?(a+b-MOD):(a+b); }
inline int dis(int a,int b) { return (a-b<0)?(a-b+MOD):(a-b); }
//求余意义下的逆元、加减法

inline pii operator ~ (const pii &p) { ll tmp=inv(p.fi); return pii(tmp,dis(0,tmp*tmp%MOD*p.se%MOD)); }
//逆元
inline pii operator - (const pii &p) { return pii( dis(0,p.fi) , dis(0,p.se) ); }
inline pii operator + (const pii &a,const pii &b) { return pii( add(a.fi,b.fi) , add(a.se,b.se) ); }
inline pii operator += (pii &a,const pii &b) { return a=a+b; }
inline pii operator - (const pii &a,const pii &b) { return pii( dis(a.fi,b.fi) , dis(a.se,b.se) ); }
inline pii operator -= (pii &a,const pii &b) { return a=a-b; }
inline pii operator * (const pii &a,const pii &b) { return pii( (ll)a.fi*b.fi%MOD , ((ll)a.fi*b.se+(ll)a.se*b.fi)%MOD ); }
inline pii operator *= (pii &a,const pii &b) { return a=a*b; }
inline pii operator / (const pii &a,const pii &b) { return a*(~b); }
inline pii operator /= (pii &a,const pii &b) { return a=a/b; }
//重载运算符

struct Matrix{
    pii Num[MAXN][MAXN];
    inline pii det(int N){//求行列式
        pii res=pii(1,0);
        for(register int i=1;i<=N;++i){
            int pos=i;
            while(pos<=N&&Num[pos][i]==pii(0,0)) ++pos;
            if(pos>N) return pii(0,0);
            if(pos!=i) swap(Num[pos],Num[i]),res=-res;
            res*=Num[i][i];
            for(register int j=N;j>=i;j--) Num[i][j]/=Num[i][i];

            for(register int k=i+1;k<=N;++k)
                for(register int j=N;j>=i;--j)
                    Num[k][j]-=Num[i][j]*Num[k][i];
        }
        return res;
    }
}M;
struct Edge{
    int u,v,w;
    Edge(int u_=0,int v_=0,int w_=0):u(u_),v(v_),w(w_) {}
};
int n,Pa[MAXN],CntSet;
int find(int u) { return (u==Pa[u])?u:(Pa[u]=find(Pa[u])); }
inline void merge(int u,int v){
    if(find(u)==find(v)) return ;
    Pa[find(u)]=find(v);
    --CntSet;
}
//并查集相关操作
inline int calc(const vector<Edge> &v){
    //给定 vector ，求答案
    if(v.size()<n-1) return 0;
    //一定不成树
    iota(Pa+1,Pa+n+1,1); CntSet=n;
    memset(M.Num,0,sizeof(M.Num));
    //并查集和矩阵初始化
    for(const auto &e : v)
        merge(e.u,e.v),
        M.Num[e.u][e.v]-=pii(1,e.w),
        M.Num[e.v][e.u]-=pii(1,e.w),
        M.Num[e.u][e.u]+=pii(1,e.w),
        M.Num[e.v][e.v]+=pii(1,e.w);
    if(CntSet!=1) return 0;//无法生成树
    return M.det(n-1).se;
}
vector<Edge> e[Lim+1];
ll fc[Lim+1],phi[Lim+1],prime[Lim+1],cntprime;
inline void sieve(){
    //欧拉筛积性函数
    phi[1]=1;
    for(int i=2;i<=Lim;i++){
        if(!fc[i]) phi[i]=(fc[i]=prime[++cntprime]=i)-1;
        for(int j=1;j<=cntprime;j++)
            if(i*prime[j]>Lim) break;
            else if(fc[i]>prime[j]) fc[i*prime[j]]=prime[j],phi[i*prime[j]]=phi[i]*(prime[j]-1);
            else{
                fc[i*prime[j]]=prime[j];
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
    }
}
inline void pre(){
    int m;
    cin>>n>>m;
    for(int i=1,u,v,w;i<=m;i++){
        cin>>u>>v>>w;
        for(int d=1;d*d<=w;d++) if(w%d==0) {
            e[d].push_back( Edge(u,v,w) );
            if(w/d!=d) e[w/d].push_back( Edge(u,v,w) );
        }
        //给因数打上标记
    }
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    pre();
    sieve();
    ll ans=0;
    for(int i=1;i<=Lim;i++) ans+=phi[i]*calc(e[i])%MOD;
    ans%=MOD;
    cout<<ans;
    cout.flush();
    return 0;
}