[BZOJ 5093] 图的价值

5093: [Lydsy1711月赛]图的价值

Time Limit: 30 Sec  Memory Limit: 256 MB
Submit: 544  Solved: 273
[Submit][Status][Discuss]

Description

“简单无向图”是指无重边、无自环的无向图（不一定连通）。

一个带标号的图的价值定义为每个点度数的k次方的和。

给定n和k，请计算所有n个点的带标号的简单无向图的价值之和。

因为答案很大，请对998244353取模输出。

 

Input

第一行包含两个正整数n,k(1<=n<=10^9,1<=k<=200000)。

 

Output

 输出一行一个整数，即答案对998244353取模的结果。

 

Sample Input

6 5

Sample Output

67584000

开局一个式子, 装备全靠混凝土

首先我们发现贡献可以非常容易地根据每个点的情况分开来算, 我们有:
$$
F(n)=nsum_{i=1}^{n-1}{n-1 choose i}i^k2^{{nchoose 2}-(n-1)}
$$
含义就是枚举某个点的度数, 算出和它相关的连边方案, 其他边随便连. 点和点之间等价, 直接乘 $n$.

接着我们发现和式内部有些和 $i$ 无关的东西, 我们把它扥出来:
$$
F(n)=n2^{{nchoose 2}-(n-1)}sum_{i=1}^{n-1}{n-1choose i}i^k
$$
前面部分显然非常好求, 我们把重点放在后面那个和式上, 设:
$$
G(n)=sum_{i=1}^n{nchoose i}i^k
$$
这个求和式的项数是和 $n$ 有关的, 然而 $n$ 巨大无比根本跑不出来...只有一个 $k$ 的数据范围还有救...

我们考虑如何把求和上界向 $k$ 的方向转化.

翻阅混凝土数学, 我们找到了第二类Stirling反演公式 $(6.10)$:
$$
x^k=sum_{i=0}^kegin{Bmatrix}k\iend{Bmatrix}x^{underline i}
$$
扔进去试试看:
$$
G(n)=sum_{i=1}^n{nchoose i}sum_{r=0}^kegin{Bmatrix}k\rend{Bmatrix}i^{underline r}
$$
按照求和式套路把所有东西放在最里层, 对外层求和指标乱搞再把和里层指标无关的东西抻回来:
$$
egin{aligned}
G(n)&=sum_{i=1}^nsum_{r=0}^k{nchoose i}egin{Bmatrix}k\rend{Bmatrix}i^{underline r} \
&=sum_{r=0}^ksum_{i=1}^n{nchoose i}egin{Bmatrix}k\rend{Bmatrix}i^{underline r} \
&=sum_{r=0}^kegin{Bmatrix}k\rend{Bmatrix}sum_{i=1}^n{nchoose i}i^{underline r}
end{aligned}
$$
我们看那个下降阶乘幂在一个二项式系数旁边就非常不爽, 我们把它换成二项式系数来方便比对公式:
$$
G(n)=sum_{r=0}^kegin{Bmatrix}k\rend{Bmatrix}sum_{i=1}^n{nchoose i}{ichoose r}r!
$$
然后我们果然在混凝土中匹配到了一个公式 $(5.21)$!
$$
{rchoose m}{mchoose k}={rchoose k}{r-kchoose m-k}
$$
代进去化一化:
$$
egin{aligned}
G(n)&=sum_{r=0}^kegin{Bmatrix}k\rend{Bmatrix}sum_{i=1}^n{nchoose r}{n-rchoose i-r}r!\
&=sum_{r=0}^kegin{Bmatrix}k\rend{Bmatrix}r!{nchoose r}sum_{i=1}^n{n-rchoose i-r}\
&=sum_{r=0}^kegin{Bmatrix}k\rend{Bmatrix}r!{nchoose r}2^{n-r}
end{aligned}
$$
完美!

现在问题是怎么求 $egin{Bmatrix}k\rend{Bmatrix}$. 它可能并没有什么封闭形式.

回想 $egin{Bmatrix}k\rend{Bmatrix}$ 到底表示什么:

符号 $egin{Bmatrix}n\kend{Bmatrix}$ 表示将一个有 $n$ 个物品的集合划分成 $k$ 个非空子集的方案数

如果没有集合非空的限制而且集合带有 1~k 的标号, 那么我们可以发现答案就是 $k^n$(给 $n$ 个元素分配一个集合标号)

枚举有多少集合是空的, 容斥一下并消序得到:
$$
egin{Bmatrix}n\kend{Bmatrix}=frac 1 {k!}sum_{i=0}^k{kchoose i}(-1)^i(k-i)^n
$$
我们发现后面的和式就是数列 $left langle (-1)^k ight angle$ 和 $left langle k^n ight angle$ 的二项卷积, 法法塔一发即可 $O(nlog n)$ 求出所有 $egin{Bmatrix} n \ k end{Bmatrix}$.

求完Stirling数回代即可. 复杂度 $O(klog k)$.

代码实现

就用了个多项式乘法还套全家桶真是吃枣药丸

#include <bits/stdc++.h>

const int G=3;
const int DFT=1;
const int IDFT=-1;
const int MAXN=550000;
const int MOD=998244353;
const int INV2=(MOD+1)>>1;
const int PHI=MOD-1;

typedef std::vector<int> Poly;

Poly Sqrt(Poly);
void Read(Poly&);
Poly Inverse(Poly);
Poly Ln(const Poly&);
Poly Exp(const Poly&);
void Print(const Poly&);
void NTT(Poly&,int,int);
Poly Pow(const Poly&,int);
Poly Integral(const Poly&);
Poly Derivative(const Poly&);
Poly operator*(Poly,Poly);
Poly operator/(Poly,Poly);
Poly operator%(Poly,Poly);
Poly operator+(const Poly&,const Poly&);
Poly operator-(const Poly&,const Poly&);

int Cn[MAXN];
int rev[MAXN];
int inv[MAXN];
int fact[MAXN];

int NTTPre(int);
int Sqrt(int,int);
int Pow(int,int,int);
int Log(int,int,int);
int ExGCD(int,int,int&,int&);

int main(){
    int n,k;
    scanf("%d%d",&n,&k);
    int ans=1ll*n*Pow(2,1ll*(n-1)*(n-2)/2%(MOD-1),MOD)%MOD;
    --n;
    Poly a(k+1),b(k+1);
    fact[0]=a[0]=Cn[0]=1;
    for(int i=1;i<=k;i++){
        fact[i]=1ll*fact[i-1]*i%MOD;
        inv[i]=Pow(fact[i],MOD-2,MOD);
        a[i]=1ll*(i&1?MOD-1:1)*inv[i]%MOD;
        b[i]=1ll*Pow(i,k,MOD)*inv[i]%MOD;
        Cn[i]=1ll*Cn[i-1]*(n-i+1)%MOD*Pow(i,MOD-2,MOD)%MOD;
    }
    Poly s=a*b;
    int g=0;
    for(int i=0;i<=k;i++)
        (g+=1ll*s[i]*fact[i]%MOD*Cn[i]%MOD*Pow(2,n-i,MOD)%MOD)%=MOD;
    ans=1ll*g*ans%MOD;
    printf("%d
",ans);
    return 0;
}

void Read(Poly& a){
    for(auto& i:a)
        scanf("%d",&i);
}

void Print(const Poly& a){
    for(auto i:a)
        printf("%d ",i);
    puts("");
}

Poly Pow(const Poly& a,int k){
    Poly log=Ln(a);
    for(auto& i:log)
        i=1ll*i*k%MOD;
    return Exp(log);
}

Poly Sqrt(Poly a){
    int len=a.size();
    if(len==1)
        return Poly(1,Sqrt(a[0],MOD));
    else{
        Poly b=a;
        b.resize((len+1)>>1);
        b=Sqrt(b);
        b.resize(len);
        Poly inv=Inverse(b);
        int bln=NTTPre(inv.size()+a.size());
        NTT(a,bln,DFT);
        NTT(inv,bln,DFT);
        for(int i=0;i<bln;i++)
            a[i]=1ll*a[i]*INV2%MOD*inv[i]%MOD;
        NTT(a,bln,IDFT);
        for(int i=0;i<len;i++)
            b[i]=(1ll*b[i]*INV2%MOD+a[i])%MOD;
        return b;
    }
}

Poly Exp(const Poly& a){
    size_t len=1;
    Poly ans(1,1),one(1,1);
    while(len<(a.size()<<1)){
        len<<=1;
        Poly b=a;
        b.resize(len);
        ans=ans*(one-Ln(ans)+b);
        ans.resize(len);
    }
    ans.resize(a.size());
    return ans;
}

Poly Ln(const Poly& a){
    Poly ans=Integral(Derivative(a)*Inverse(a));
    ans.resize(a.size());
    return ans;
}

Poly Integral(const Poly& a){
    int len=a.size();
    Poly ans(len+1);
    for(int i=1;i<len;i++)
        ans[i]=1ll*a[i-1]*Pow(i,MOD-2,MOD)%MOD;
    return ans;
}

Poly Derivative(const Poly& a){
    int len=a.size();
    Poly ans(len-1);
    for(int i=1;i<len;i++)
        ans[i-1]=1ll*a[i]*i%MOD;
    return ans;
}

Poly operator/(Poly a,Poly b){
    int n=a.size()-1,m=b.size()-1;
    Poly ans(1);
    if(n>=m){
        std::reverse(a.begin(),a.end());
        std::reverse(b.begin(),b.end());
        b.resize(n-m+1);
        ans=Inverse(b)*a;
        ans.resize(n-m+1);
        std::reverse(ans.begin(),ans.end());
    }
    return ans;
}

Poly operator%(Poly a,Poly b){
    int n=a.size()-1,m=b.size()-1;
    Poly ans;
    if(n<m)
        ans=a;
    else
        ans=a-(a/b)*b;
    ans.resize(m);
    return ans;
}

Poly operator*(Poly a,Poly b){
    int len=a.size()+b.size()-1;
    int bln=NTTPre(len);
    NTT(a,bln,DFT);
    NTT(b,bln,DFT);
    for(int i=0;i<bln;i++)
        a[i]=1ll*a[i]*b[i]%MOD;
    NTT(a,bln,IDFT);
    a.resize(len);
    return a;
}

Poly operator+(const Poly& a,const Poly& b){
    Poly ans(std::max(a.size(),b.size()));
    std::copy(a.begin(),a.end(),ans.begin());
    for(size_t i=0;i<b.size();i++)
        ans[i]=(ans[i]+b[i])%MOD;
    return ans;
}

Poly operator-(const Poly& a,const Poly& b){
    Poly ans(std::max(a.size(),b.size()));
    std::copy(a.begin(),a.end(),ans.begin());
    for(size_t i=0;i<b.size();i++)
        ans[i]=(ans[i]+MOD-b[i])%MOD;
    return ans;
}

Poly Inverse(Poly a){
    int len=a.size();
    if(len==1)
        return Poly(1,Pow(a[0],MOD-2,MOD));
    else{
        Poly b(a);
        b.resize((len+1)>>1);
        b=Inverse(b);
        int bln=NTTPre(b.size()*2+a.size());
        NTT(a,bln,DFT);
        NTT(b,bln,DFT);
        for(int i=0;i<bln;i++)
            b[i]=(2ll*b[i]%MOD-1ll*b[i]*b[i]%MOD*a[i]%MOD+MOD)%MOD;
        NTT(b,bln,IDFT);
        b.resize(len);
        return b;
    }
}

void NTT(Poly& a,int len,int opt){
    a.resize(len);
    for(int i=0;i<len;i++)
        if(rev[i]>i)
            std::swap(a[i],a[rev[i]]);
    for(int i=1;i<len;i<<=1){
        int step=i<<1;
        int wn=Pow(G,(PHI+opt*PHI/step)%PHI,MOD);
        for(int j=0;j<len;j+=step){
            int w=1;
            for(int k=0;k<i;k++,w=1ll*w*wn%MOD){
                int x=a[j+k];
                int y=1ll*a[j+k+i]*w%MOD;
                a[j+k]=(x+y)%MOD;
                a[j+k+i]=(x-y+MOD)%MOD;
            }
        }
    }
    if(opt==IDFT){
        int inv=Pow(len,MOD-2,MOD);
        for(int i=0;i<len;i++)
            a[i]=1ll*a[i]*inv%MOD;
    }
}

int NTTPre(int n){
    int bln=1,bct=0;
    while(bln<n){
        bln<<=1;
        ++bct;
    }
    for(int i=0;i<bln;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1));
    return bln;
}

inline int Pow(int a,int n,int p){
    int ans=1;
    while(n>0){
        if(n&1)
            ans=1ll*a*ans%p;
        a=1ll*a*a%p;
        n>>=1;
    }
    return ans;
}

int ExGCD(int a,int b,int& x,int& y){
    if(b==0){
        x=1,y=0;
        return a;
    }
    else{
        int gcd=ExGCD(b,a%b,y,x);
        y-=x*(a/b);
        return gcd;
    }
}

inline int Sqrt(int a,int p){
    int s=Pow(G,Log(G,a,p)>>1,p);
    return std::min(s,MOD-s);
}

inline int Log(int a,int x,int p){
    int s=sqrt(p+0.5);
    int inv=Pow(Pow(a,s,p),p-2,p);
    std::unordered_map<int,int> m;
    m[1]=0;
    int pow=1;
    for(int i=1;i<s;i++){
        pow=1ll*a*pow%p;
        if(!m.count(pow))
            m[pow]=i;
    }
    for(int i=0;i<s;i++){
        if(m.count(x))
            return i*s+m[x];
        x=1ll*x*inv%MOD;
    }
    return -1;
}