「联合省选 2020 A」作业题 做题心得

「联合省选 2020 A」作业题 做题心得

第一次遇到这样套了好几层的数学题。代码我写了6K多,比ds还多(

但其实真的全都是套路,就好像把怎么做明摆着写出来了一样。顺着解下去,一层一层的套,就搞出来了。

这题的代码难度也很高,对于调代码也有很大锻炼

符号与定义

我瞎写的

对于图 (G)

  • (tree(G)) 表示 (G) 的生成树集合
  • (E(G)) 表示 (G) 的边集,(V(G)) 表示 (G) 的点集

对于树 (T)

  • (g|T)(T) 的每条边都是 (g) 的倍数

  • (T/g)(T) 的每条边除以 (g) 之后得到的树;它当且仅当 (g|T) 时有意义

  • $sum(T) $ :(T) 的边权和

  • (gcd(T))(T) 的边权的 (gcd)

题解

这题分两步:反演,矩阵树定理转化

反演

首先看到这个gcd第一反应显然是反演

[F(g)=sumlimits_{Tin tree(G)} [gcd(T)=g] sum(T) \ G(g)=sumlimits_{Tin tree(G)} [g|T] sum(T) ]

那么显然有 (G(g)=sumlimits_{g|h} F(h))

反演得 (F(g)=sumlimits_{g|h} G(h)mu(dfrac{h}{g}))

答案为

[sumlimits_{g=1}^{n} gF(g)\ =sumlimits_{g=1}^{n} gsumlimits_{g|h} G(h) mu(dfrac{h}{g})\ =sumlimits_{h=1}^{n} G(h) sumlimits_{g|h} g*mu(dfrac{h}{g}) ]

注意到后面那个玩意就是 (mu*Id),即 (varphi)

那答案就是 (sumlimits_{i=1}^{n} G(i)varphi(i)) ,能把 (G) 求出来即可

(G) 听起来很好求,就是把图中 (i) 的倍数的边拎出来,求所有生成树 (T) 的边权和的和

矩阵树定理转化

肯定考虑矩阵树定理,但是这个玩意解决的是所有生成树的边权的积的和

那现在相当于,要想办法把“乘法”转换成“加法”,并且保持普通的加法还是加法。

我一开始想到的是,都 (exp) 一遍。但是这样相当于要求所有生成树的积的积然后 (ln) 回来,可依然不是矩阵树定理的形式

同样考虑一个乘起来能变成和的形式:((1+ax)(1+bx)equiv 1+(a+b)x pmod{x^2})

换句话说,我们把每条边(权为 (w))看成一个形式幂级数:(1+wx)

然后把它们在模 (x^2) 意义下乘起来,就能得到所有生成树的边权和了,并且形式幂级数加还是照样加,没有问题 (不像 (exp) 那样都转换掉了)

总结

关于思维过程:看到gcd就反演是常规操作吧,然后接下来就是一个经典题,全jb是套路

那这个题有什么收获呢?其实主要还是调代码的能力的提升...详情见后面

#include <bits/stdc++.h>
using namespace std;
namespace Flandre_Scarlet
{
    #define N 100
    #define V 200005
    #define int long long
    #define mod 998244353
    #define F(i,l,r) for(int i=l;i<=r;++i)
    #define D(i,r,l) for(int i=r;i>=l;--i)
    #define Fs(i,l,r,c) for(int i=l;i<=r;c)
    #define Ds(i,r,l,c) for(int i=r;i>=l;c)
    #define MEM(x,a) memset(x,a,sizeof(x))
    #define FK(x) MEM(x,0)
    #define Tra(i,u) for(int i=G.st(u),v=G.to(i);~i;i=G.nx(i),v=G.to(i))
    #define p_b push_back
    #define sz(a) ((int)a.size())
    #define all(a) a.begin(),a.end()
    #define iter(a,p) (a.begin()+p)
    #define PUT(a,n) F(i,1,n) printf("%d ",a[i]); puts("");
    int I() {char c=getchar(); int x=0; int f=1; while(c<'0' or c>'9') f=(c=='-')?-1:1,c=getchar(); while(c>='0' and c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar(); return ((f==1)?x:-x);}
    template <typename T> void Rd(T& arg){arg=I();}
    template <typename T,typename...Types> void Rd(T& arg,Types&...args){arg=I(); Rd(args...);}
    void RA(int *p,int n) {F(i,1,n) *p=I(),++p;}
    int n,m;
    int g[N][N];
    vector<int> edges;
    void Input()
    {
        Rd(n,m);
        F(i,1,m)
        {
            int u,v,w; Rd(u,v,w);
            g[u][v]=g[v][u]=w;
            edges.p_b(w);
        }
    }
    int qpow(int a,int b,int m=mod) {int r=1; while(b) {if (b&1) r=r*a%m; a=a*a%m,b>>=1;} return r;}
    int primes[V]; bool notp[V]; int phi[V];
    void Init()
    {
        int n=2e5;
        int& cnt=primes[0];
        notp[1]=phi[1]=1;
        F(i,2,n)
        {
            if (!notp[i])
            {
                primes[++cnt]=i;
                phi[i]=i-1;
            }
            for(int j=1;j<=cnt and i*primes[j]<=n;++j)
            {
                int u=primes[j];
                notp[i*u]=1;
                if (i%u==0) 
                {
                    phi[i*u]=phi[i]*u;
                    break;
                }
                else
                {
                    phi[i*u]=phi[i]*(u-1);
                }
            }
        }
    }
    struct node
    {
        int a,b;
        node(int w=0){a=1,b=w;}
        node(int x,int y) {a=x,b=y;}
    }; // poly: a+bx, %(x^2)
    node operator-(node a) {return (node){mod-a.a,mod-a.b};}
    node operator+(node a,node b) {return (node){(a.a+b.a)%mod,(a.b+b.b)%mod};}
    node operator-(node a,node b) {return a+(-b);}
    node operator*(node a,node b) {return (node){(a.a*b.a)%mod,(a.a*b.b+a.b*b.a)%mod};}
    node inv(node a) {int iv=qpow(a.a,mod-2,mod); return (node){iv,(mod-a.b)*iv%mod*iv%mod};}
    node operator/(node a,node b) {return a*inv(b);}
    namespace Linear_Algebra
    {
        
        class Matrix
        {
        public:
            node a[N][N]; int n;
            node* operator[](int id) {return a[id];}
            void lswap(int i,int j) {F(x,1,n) swap(a[i][x],a[j][x]);}
            void ltime(int i,node k) {F(x,1,n) a[i][x]=(a[i][x]*k);}
            void laddt(int i,node k,int j) {F(x,1,n) a[j][x]=(a[j][x]+a[i][x]*k);}
            void put()
            {
                F(i,1,n)
                {
                    F(j,1,n) printf("(%lld,%lld) ",a[i][j].a,a[i][j].b);
                    puts("");
                }
            }
            int det()
            {
                node ans=(node){1,0};
                F(i,1,n) 
                {
                    bool all0=1;
                    F(j,i,n) if (a[j][i].a!=0) {all0=0; break;}
                    if (all0) // 按 b 消元
                    {
                        if (a[i][i].b==0)
                        {
                            F(j,i+1,n) if (a[j][i].b!=0) 
                            {
                                lswap(i,j); ans=(-ans);
                                break;
                            }
                        }
                        F(j,i+1,n)
                        {
                            int k=(mod-a[j][i].b)*qpow(a[i][i].b,mod-2,mod)%mod;
                            laddt(i,(node){k,0},j);
                        }
                    }
                    else // 按 a 消元
                    {
                        if (a[i][i].a==0)
                        {
                            F(j,i+1,n) if (a[j][i].a!=0)
                            {
                                lswap(i,j); ans=(-ans);
                                break;
                            }
                        }
                        F(j,i+1,n)
                        {
                            node k=(-a[j][i])/(a[i][i]);
                            laddt(i,k,j);
                        }
                    }
                }
                F(i,1,n) ans=ans*a[i][i];
                return ans.b;
            }
        }Lg;
        node d[N];
        int EdgeSum(node g2[N][N]) // edge sum of spanning tree
        {
            F(i,1,n)
            {
                d[i]=node(0,0);
                F(j,1,n) if (j!=i)
                {
                    d[i]=d[i]+g2[i][j];
                }
            }
            F(i,1,n-1) Lg[i][i]=d[i];
            F(i,1,n-1) F(j,1,n-1) if (i!=j) Lg[i][j]=-g2[i][j];
            Lg.n=n-1;
            return Lg.det();
        }
    }
    bool cxk[V];
    node g2[N][N];
    struct small_union_find
    {
        int fa[N]; 
        void clear() {F(i,1,n) fa[i]=i;}
        int find(int u) {return u==fa[u]?u:fa[u]=find(fa[u]);}
        void merge(int u,int v) {fa[find(u)]=find(v);}
        bool con() {F(i,1,n) if (find(i)!=find(1)) return false; return true;}
    }un;
    void Sakuya()
    {
        Init();
        F(d,1,V) 
        {
            for(auto e:edges) if (e%d==0) 
            {
                cxk[d]=1;
            }
        }
        int ans=0;
        F(d,1,V) if (cxk[d])
        {
            un.clear();
            F(i,1,n) F(j,1,n) 
            {
                if (g[i][j] and g[i][j]%d==0) 
                {
                    un.merge(i,j);
                    g2[i][j]=node(1,g[i][j]);
                }
                else
                {
                    g2[i][j]=node(0,0);
                }
            }
            
            if (!un.con()) 
            {
                continue;
            }
            int es=Linear_Algebra::EdgeSum(g2);
            ans=(ans+phi[d]*es%mod)%mod;
        }
        printf("%lld
",ans);
    }
    void IsMyWife()
    {
        Input();
        Sakuya();
    }
}
#undef int //long long
int main()
{
    Flandre_Scarlet::IsMyWife();
    getchar();
    return 0;
}

后记:关于考场调试

  1. 不会吧不会吧还有人筛 (varphi) 会写错?

  2. 把那个 (pmod{x^2}) 的形式幂级数打成 struct 先给写了,然后写好矩阵

    此时可以把样例的图输入进去试试,不去处理因数啥的,看看它是否输出是 44 (所有生成树的边权的和的和)

  3. 加上处理倍数的代码,看看样例能否跑对

    可以选择对于每个 (i)(G(i)) 输出来,看看是否和期望的一样

    样例跑对了,手造几个小数据,看看能否跑对 (我偷懒了,直接去下了数据)

还有就是一个细节要注意:写高斯消元及相关玩意的时候,一定要注意关于 (0) 的讨论,尽管非常麻烦,但不能不讨论 —— 不然会 WA 自闭

原文地址:https://www.cnblogs.com/LightningUZ/p/14635318.html