[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)=n\sum_{i=1}^{n-1}{n-1 \choose i}i^k2^{{n\choose 2}-(n-1)}
$$
含义就是枚举某个点的度数, 算出和它相关的连边方案, 其他边随便连. 点和点之间等价, 直接乘 $n$.

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

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

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

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

回想 $\begin{Bmatrix}k\\r\end{Bmatrix}$ 到底表示什么:

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

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

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

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

代码实现

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

  1 #include <bits/stdc++.h>
  2 
  3 const int G=3;
  4 const int DFT=1;
  5 const int IDFT=-1;
  6 const int MAXN=550000;
  7 const int MOD=998244353;
  8 const int INV2=(MOD+1)>>1;
  9 const int PHI=MOD-1;
 10 
 11 typedef std::vector<int> Poly;
 12 
 13 Poly Sqrt(Poly);
 14 void Read(Poly&);
 15 Poly Inverse(Poly);
 16 Poly Ln(const Poly&);
 17 Poly Exp(const Poly&);
 18 void Print(const Poly&);
 19 void NTT(Poly&,int,int);
 20 Poly Pow(const Poly&,int);
 21 Poly Integral(const Poly&);
 22 Poly Derivative(const Poly&);
 23 Poly operator*(Poly,Poly);
 24 Poly operator/(Poly,Poly);
 25 Poly operator%(Poly,Poly);
 26 Poly operator+(const Poly&,const Poly&);
 27 Poly operator-(const Poly&,const Poly&);
 28 
 29 int Cn[MAXN];
 30 int rev[MAXN];
 31 int inv[MAXN];
 32 int fact[MAXN];
 33 
 34 int NTTPre(int);
 35 int Sqrt(int,int);
 36 int Pow(int,int,int);
 37 int Log(int,int,int);
 38 int ExGCD(int,int,int&,int&);
 39 
 40 int main(){
 41     int n,k;
 42     scanf("%d%d",&n,&k);
 43     int ans=1ll*n*Pow(2,1ll*(n-1)*(n-2)/2%(MOD-1),MOD)%MOD;
 44     --n;
 45     Poly a(k+1),b(k+1);
 46     fact[0]=a[0]=Cn[0]=1;
 47     for(int i=1;i<=k;i++){
 48         fact[i]=1ll*fact[i-1]*i%MOD;
 49         inv[i]=Pow(fact[i],MOD-2,MOD);
 50         a[i]=1ll*(i&1?MOD-1:1)*inv[i]%MOD;
 51         b[i]=1ll*Pow(i,k,MOD)*inv[i]%MOD;
 52         Cn[i]=1ll*Cn[i-1]*(n-i+1)%MOD*Pow(i,MOD-2,MOD)%MOD;
 53     }
 54     Poly s=a*b;
 55     int g=0;
 56     for(int i=0;i<=k;i++)
 57         (g+=1ll*s[i]*fact[i]%MOD*Cn[i]%MOD*Pow(2,n-i,MOD)%MOD)%=MOD;
 58     ans=1ll*g*ans%MOD;
 59     printf("%d\n",ans);
 60     return 0;
 61 }
 62 
 63 void Read(Poly& a){
 64     for(auto& i:a)
 65         scanf("%d",&i);
 66 }
 67 
 68 void Print(const Poly& a){
 69     for(auto i:a)
 70         printf("%d ",i);
 71     puts("");
 72 }
 73 
 74 Poly Pow(const Poly& a,int k){
 75     Poly log=Ln(a);
 76     for(auto& i:log)
 77         i=1ll*i*k%MOD;
 78     return Exp(log);
 79 }
 80 
 81 Poly Sqrt(Poly a){
 82     int len=a.size();
 83     if(len==1)
 84         return Poly(1,Sqrt(a[0],MOD));
 85     else{
 86         Poly b=a;
 87         b.resize((len+1)>>1);
 88         b=Sqrt(b);
 89         b.resize(len);
 90         Poly inv=Inverse(b);
 91         int bln=NTTPre(inv.size()+a.size());
 92         NTT(a,bln,DFT);
 93         NTT(inv,bln,DFT);
 94         for(int i=0;i<bln;i++)
 95             a[i]=1ll*a[i]*INV2%MOD*inv[i]%MOD;
 96         NTT(a,bln,IDFT);
 97         for(int i=0;i<len;i++)
 98             b[i]=(1ll*b[i]*INV2%MOD+a[i])%MOD;
 99         return b;
100     }
101 }
102 
103 Poly Exp(const Poly& a){
104     size_t len=1;
105     Poly ans(1,1),one(1,1);
106     while(len<(a.size()<<1)){
107         len<<=1;
108         Poly b=a;
109         b.resize(len);
110         ans=ans*(one-Ln(ans)+b);
111         ans.resize(len);
112     }
113     ans.resize(a.size());
114     return ans;
115 }
116 
117 Poly Ln(const Poly& a){
118     Poly ans=Integral(Derivative(a)*Inverse(a));
119     ans.resize(a.size());
120     return ans;
121 }
122 
123 Poly Integral(const Poly& a){
124     int len=a.size();
125     Poly ans(len+1);
126     for(int i=1;i<len;i++)
127         ans[i]=1ll*a[i-1]*Pow(i,MOD-2,MOD)%MOD;
128     return ans;
129 }
130 
131 Poly Derivative(const Poly& a){
132     int len=a.size();
133     Poly ans(len-1);
134     for(int i=1;i<len;i++)
135         ans[i-1]=1ll*a[i]*i%MOD;
136     return ans;
137 }
138 
139 Poly operator/(Poly a,Poly b){
140     int n=a.size()-1,m=b.size()-1;
141     Poly ans(1);
142     if(n>=m){
143         std::reverse(a.begin(),a.end());
144         std::reverse(b.begin(),b.end());
145         b.resize(n-m+1);
146         ans=Inverse(b)*a;
147         ans.resize(n-m+1);
148         std::reverse(ans.begin(),ans.end());
149     }
150     return ans;
151 }
152 
153 Poly operator%(Poly a,Poly b){
154     int n=a.size()-1,m=b.size()-1;
155     Poly ans;
156     if(n<m)
157         ans=a;
158     else
159         ans=a-(a/b)*b;
160     ans.resize(m);
161     return ans;
162 }
163 
164 Poly operator*(Poly a,Poly b){
165     int len=a.size()+b.size()-1;
166     int bln=NTTPre(len);
167     NTT(a,bln,DFT);
168     NTT(b,bln,DFT);
169     for(int i=0;i<bln;i++)
170         a[i]=1ll*a[i]*b[i]%MOD;
171     NTT(a,bln,IDFT);
172     a.resize(len);
173     return a;
174 }
175 
176 Poly operator+(const Poly& a,const Poly& b){
177     Poly ans(std::max(a.size(),b.size()));
178     std::copy(a.begin(),a.end(),ans.begin());
179     for(size_t i=0;i<b.size();i++)
180         ans[i]=(ans[i]+b[i])%MOD;
181     return ans;
182 }
183 
184 Poly operator-(const Poly& a,const Poly& b){
185     Poly ans(std::max(a.size(),b.size()));
186     std::copy(a.begin(),a.end(),ans.begin());
187     for(size_t i=0;i<b.size();i++)
188         ans[i]=(ans[i]+MOD-b[i])%MOD;
189     return ans;
190 }
191 
192 Poly Inverse(Poly a){
193     int len=a.size();
194     if(len==1)
195         return Poly(1,Pow(a[0],MOD-2,MOD));
196     else{
197         Poly b(a);
198         b.resize((len+1)>>1);
199         b=Inverse(b);
200         int bln=NTTPre(b.size()*2+a.size());
201         NTT(a,bln,DFT);
202         NTT(b,bln,DFT);
203         for(int i=0;i<bln;i++)
204             b[i]=(2ll*b[i]%MOD-1ll*b[i]*b[i]%MOD*a[i]%MOD+MOD)%MOD;
205         NTT(b,bln,IDFT);
206         b.resize(len);
207         return b;
208     }
209 }
210 
211 void NTT(Poly& a,int len,int opt){
212     a.resize(len);
213     for(int i=0;i<len;i++)
214         if(rev[i]>i)
215             std::swap(a[i],a[rev[i]]);
216     for(int i=1;i<len;i<<=1){
217         int step=i<<1;
218         int wn=Pow(G,(PHI+opt*PHI/step)%PHI,MOD);
219         for(int j=0;j<len;j+=step){
220             int w=1;
221             for(int k=0;k<i;k++,w=1ll*w*wn%MOD){
222                 int x=a[j+k];
223                 int y=1ll*a[j+k+i]*w%MOD;
224                 a[j+k]=(x+y)%MOD;
225                 a[j+k+i]=(x-y+MOD)%MOD;
226             }
227         }
228     }
229     if(opt==IDFT){
230         int inv=Pow(len,MOD-2,MOD);
231         for(int i=0;i<len;i++)
232             a[i]=1ll*a[i]*inv%MOD;
233     }
234 }
235 
236 int NTTPre(int n){
237     int bln=1,bct=0;
238     while(bln<n){
239         bln<<=1;
240         ++bct;
241     }
242     for(int i=0;i<bln;i++)
243         rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1));
244     return bln;
245 }
246 
247 inline int Pow(int a,int n,int p){
248     int ans=1;
249     while(n>0){
250         if(n&1)
251             ans=1ll*a*ans%p;
252         a=1ll*a*a%p;
253         n>>=1;
254     }
255     return ans;
256 }
257 
258 int ExGCD(int a,int b,int& x,int& y){
259     if(b==0){
260         x=1,y=0;
261         return a;
262     }
263     else{
264         int gcd=ExGCD(b,a%b,y,x);
265         y-=x*(a/b);
266         return gcd;
267     }
268 }
269 
270 inline int Sqrt(int a,int p){
271     int s=Pow(G,Log(G,a,p)>>1,p);
272     return std::min(s,MOD-s);
273 }
274 
275 inline int Log(int a,int x,int p){
276     int s=sqrt(p+0.5);
277     int inv=Pow(Pow(a,s,p),p-2,p);
278     std::unordered_map<int,int> m;
279     m[1]=0;
280     int pow=1;
281     for(int i=1;i<s;i++){
282         pow=1ll*a*pow%p;
283         if(!m.count(pow))
284             m[pow]=i;
285     }
286     for(int i=0;i<s;i++){
287         if(m.count(x))
288             return i*s+m[x];
289         x=1ll*x*inv%MOD;
290     }
291     return -1;
292 }
BZOJ 5093

日常图包

posted @ 2019-01-09 09:20  rvalue  阅读(415)  评论(0)  编辑  收藏  举报