[BZOJ 3456] 城市规划
3456. 城市规划
Description
刚刚解决完电力网络的问题, 阿狸又被领导的任务给难住了.
刚才说过, 阿狸的国家有n个城市, 现在国家需要在某些城市对之间建立一些贸易路线, 使得整个国家的任意两个城市都直接或间接的连通. 为了省钱, 每两个城市之间最多只能有一条直接的贸易路径. 对于两个建立路线的方案, 如果存在一个城市对, 在两个方案中是否建立路线不一样, 那么这两个方案就是不同的, 否则就是相同的. 现在你需要求出一共有多少不同的方案.
好了, 这就是困扰阿狸的问题. 换句话说, 你需要求出n个点的简单(无重边无自环)无向连通图数目.
由于这个数字可能非常大, 你只需要输出方案数mod 1004535809(479 * 2 ^ 21 + 1)即可.
Input
仅一行一个整数n(<=130000)
Output
仅一行一个整数, 为方案数 mod 1004535809.
Sample Input
3
Sample Output
4
Hint
对于 100%的数据, n <= 130000
题解
首先我们可以发现, 含有 $n$ 个点的有标号简单无向图的个数是 $2^{{n \choose 2}}$.
然后我们有标号简单无向图和有标号简单连通无向图的关系: 每个有标号简单无向图都可以划分为若干联通分量, 而这些联通分量就是有标号简单联通无向图.
如何体现这个划分的过程? 首先我们使用EGF来保持标号运算, 若 $n$ 个点的有标号简单联通无向图的个数的生成函数为 $F(x)$, 则 $\frac {F(x)^k} {k!}$ 即为含有 $k$ 个联通分量的简单无向图个数的生成函数.
为啥最后有一个分母 $k!$ 呢? 我们把 $k$ 个 $F(x)$ 卷起来, 每个划分方案的每个排列都会被统计一次. 因为直接卷积时每个分量是有序合并的. 而我们在用联通无向图构造无向图的时候显然不能考虑联通分量的顺序, 所以我们除一个 $k!$.
设 $n$ 个点的简单无向图的个数的EGF为 $G(x)$, 我们就有了:
$$ G(x) = \sum_{k=0}^\infty \frac {F(x)^k}{k!} $$
又根据函数 $\exp(x)=e^x$ 的泰勒展开式:
$$ e^x=\sum_{k=0}^\infty \frac{x^k}{k!} $$
我们将 $F(x)$ 代入 $x$ 中, 得到:
$$ G(x) = \sum_{k=0}^\infty \frac {F(x)^k}{k!} = \exp(F(x)) $$
于是我们就有了:
$$ F(x)=\ln(G(x)) $$
所以我们构造出 $G(x)=\sum\limits_{k=0}^\infty 2^{{k\choose 2}} \frac {x^k}{k!} $ 之后求它的 $\ln$ 就可以得到 $F(x)$ 了.
代码实现
直接上全家桶了...然而并不全都用得到...所以码长比较怪异
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=1004535809; 8 const int INV2=(MOD+1)>>1; 9 const int PHI=MOD-1; 10 11 typedef long long int64; 12 typedef std::vector<int> Poly; 13 14 Poly Sqrt(Poly); 15 void Read(Poly&); 16 Poly Inverse(Poly); 17 Poly Ln(const Poly&); 18 Poly Exp(const Poly&); 19 void Print(const Poly&); 20 void NTT(Poly&,int,int); 21 Poly Pow(const Poly&,int); 22 Poly Integral(const Poly&); 23 Poly Derivative(const Poly&); 24 Poly operator*(Poly,Poly); 25 Poly operator/(Poly,Poly); 26 Poly operator%(Poly,Poly); 27 Poly operator+(const Poly&,const Poly&); 28 Poly operator-(const Poly&,const Poly&); 29 30 int n; 31 int rev[MAXN]; 32 int fact[MAXN]; 33 34 int NTTPre(int); 35 int Pow(int,int,int); 36 int Pow(int,int64,int); 37 38 int main(){ 39 scanf("%d",&n); 40 fact[0]=1; 41 for(int i=1;i<=n;i++) 42 fact[i]=1ll*i*fact[i-1]%MOD; 43 Poly a(n+1); 44 a[0]=1; 45 for(int i=1;i<=n;i++) 46 a[i]=1ll*Pow(2,1ll*i*(i-1)/2,MOD)*Pow(fact[i],MOD-2,MOD)%MOD; 47 a=Ln(a); 48 printf("%d\n",int(1ll*a[n]*fact[n]%MOD)); 49 return 0; 50 } 51 52 void Read(Poly& a){ 53 for(auto& i:a) 54 scanf("%d",&i); 55 } 56 57 void Print(const Poly& a){ 58 for(auto i:a) 59 printf("%d ",i); 60 puts(""); 61 } 62 63 Poly Pow(const Poly& a,int k){ 64 Poly log=Ln(a); 65 for(auto& i:log) 66 i=1ll*i*k%MOD; 67 return Exp(log); 68 } 69 70 Poly Sqrt(Poly a){ 71 int len=a.size(); 72 if(len==1) 73 return Poly(1,int(sqrt(a[0]))); 74 else{ 75 Poly b=a; 76 b.resize((len+1)>>1); 77 b=Sqrt(b); 78 b.resize(len); 79 Poly inv=Inverse(b); 80 int bln=NTTPre(inv.size()+a.size()); 81 NTT(a,bln,DFT); 82 NTT(inv,bln,DFT); 83 for(int i=0;i<bln;i++) 84 a[i]=1ll*a[i]*INV2%MOD*inv[i]%MOD; 85 NTT(a,bln,IDFT); 86 for(int i=0;i<len;i++) 87 b[i]=(1ll*b[i]*INV2%MOD+a[i])%MOD; 88 return b; 89 } 90 } 91 92 Poly Exp(const Poly& a){ 93 size_t len=1; 94 Poly ans(1,1),one(1,1); 95 while(len<(a.size()<<1)){ 96 len<<=1; 97 Poly b=a; 98 b.resize(len); 99 ans=ans*(one-Ln(ans)+b); 100 ans.resize(len); 101 } 102 ans.resize(a.size()); 103 return ans; 104 } 105 106 Poly Ln(const Poly& a){ 107 Poly ans=Integral(Derivative(a)*Inverse(a)); 108 ans.resize(a.size()); 109 return ans; 110 } 111 112 Poly Integral(const Poly& a){ 113 int len=a.size(); 114 Poly ans(len+1); 115 for(int i=1;i<len;i++) 116 ans[i]=1ll*a[i-1]*Pow(i,MOD-2,MOD)%MOD; 117 return ans; 118 } 119 120 Poly Derivative(const Poly& a){ 121 int len=a.size(); 122 Poly ans(len-1); 123 for(int i=1;i<len;i++) 124 ans[i-1]=1ll*a[i]*i%MOD; 125 return ans; 126 } 127 128 Poly operator/(Poly a,Poly b){ 129 int n=a.size()-1,m=b.size()-1; 130 Poly ans(1); 131 if(n>=m){ 132 std::reverse(a.begin(),a.end()); 133 std::reverse(b.begin(),b.end()); 134 b.resize(n-m+1); 135 ans=Inverse(b)*a; 136 ans.resize(n-m+1); 137 std::reverse(ans.begin(),ans.end()); 138 } 139 return ans; 140 } 141 142 Poly operator%(Poly a,Poly b){ 143 int n=a.size()-1,m=b.size()-1; 144 Poly ans; 145 if(n<m) 146 ans=a; 147 else 148 ans=a-(a/b)*b; 149 ans.resize(m); 150 return ans; 151 } 152 153 Poly operator*(Poly a,Poly b){ 154 int len=a.size()+b.size()-1; 155 int bln=NTTPre(len); 156 NTT(a,bln,DFT); 157 NTT(b,bln,DFT); 158 for(int i=0;i<bln;i++) 159 a[i]=1ll*a[i]*b[i]%MOD; 160 NTT(a,bln,IDFT); 161 a.resize(len); 162 return a; 163 } 164 165 Poly operator+(const Poly& a,const Poly& b){ 166 Poly ans(std::max(a.size(),b.size())); 167 std::copy(a.begin(),a.end(),ans.begin()); 168 for(size_t i=0;i<b.size();i++) 169 ans[i]=(ans[i]+b[i])%MOD; 170 return ans; 171 } 172 173 Poly operator-(const Poly& a,const Poly& b){ 174 Poly ans(std::max(a.size(),b.size())); 175 std::copy(a.begin(),a.end(),ans.begin()); 176 for(size_t i=0;i<b.size();i++) 177 ans[i]=(ans[i]+MOD-b[i])%MOD; 178 return ans; 179 } 180 181 Poly Inverse(Poly a){ 182 int len=a.size(); 183 if(len==1) 184 return Poly(1,Pow(a[0],MOD-2,MOD)); 185 else{ 186 Poly b(a); 187 b.resize((len+1)>>1); 188 b=Inverse(b); 189 int bln=NTTPre(b.size()*2+a.size()); 190 NTT(a,bln,DFT); 191 NTT(b,bln,DFT); 192 for(int i=0;i<bln;i++) 193 b[i]=(2ll*b[i]%MOD-1ll*b[i]*b[i]%MOD*a[i]%MOD+MOD)%MOD; 194 NTT(b,bln,IDFT); 195 b.resize(len); 196 return b; 197 } 198 } 199 200 void NTT(Poly& a,int len,int opt){ 201 a.resize(len); 202 for(int i=0;i<len;i++) 203 if(rev[i]>i) 204 std::swap(a[i],a[rev[i]]); 205 for(int i=1;i<len;i<<=1){ 206 int step=i<<1; 207 int wn=Pow(G,(PHI+opt*PHI/step)%PHI,MOD); 208 for(int j=0;j<len;j+=step){ 209 int w=1; 210 for(int k=0;k<i;k++,w=1ll*w*wn%MOD){ 211 int x=a[j+k]; 212 int y=1ll*a[j+k+i]*w%MOD; 213 a[j+k]=(x+y)%MOD; 214 a[j+k+i]=(x-y+MOD)%MOD; 215 } 216 } 217 } 218 if(opt==IDFT){ 219 int inv=Pow(len,MOD-2,MOD); 220 for(int i=0;i<len;i++) 221 a[i]=1ll*a[i]*inv%MOD; 222 } 223 } 224 225 int NTTPre(int n){ 226 int bln=1,bct=0; 227 while(bln<n){ 228 bln<<=1; 229 ++bct; 230 } 231 for(int i=0;i<bln;i++) 232 rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1)); 233 return bln; 234 } 235 236 inline int Pow(int a,int n,int p){ 237 int ans=1; 238 while(n>0){ 239 if(n&1) 240 ans=1ll*a*ans%p; 241 a=1ll*a*a%p; 242 n>>=1; 243 } 244 return ans; 245 } 246 247 inline int Pow(int a,int64 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 }
本博客已弃用, 新个人主页: https://rvalue.moe, 新博客: https://blog.rvalue.moe