[BZOJ 5093] 圖的價值
阿新 • • 發佈：2019-01-09
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 5Sample 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
