[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)=sumlimits_{k=0}^infty 2^{{kchoose 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
",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 }

BZOJ 3456

日常图包