「BZOJ 3456」城市规划

Problem

Description

n 个点的简单无向连通图个数。

答案对 1004535809 取模。

Constraints

1\le n\le 130000

Solution

Analysis

n 个点的简单无向连通图个数为 f_{n} ,容易有:

\begin{aligned} f_{n} &= 2^{\binom{n}{2} } -\sum_{i = 1}^{n - 1}2^{\binom{n - i}{2}}\binom{n - 1}{i - 1}f_{i}\\ &= 2^{\binom{n}{2} } -(n - 1)!\sum_{i = 1}^{n - 1}\frac{2^{\binom{n - i}{2}}}{(n - i)!}\cdot\frac{f_{i}}{(i - 1)!} \end{aligned}

亦即 n 个点的简单无向图个数减去非连通图个数,后者可通过枚举点 1 所在连通块大小计算。
此递推柿显然可以用分治 FFT 计算,复杂度 O\left(n\log^{2}{n}\right)

尝试寻找更优的解法,发现有:

\begin{aligned} \frac{2^{\binom{n}{2}}}{(n - 1)!}&=\frac{f_{n}}{(n - 1)!}+\sum_{i = 1}^{n - 1}\frac{2^{\binom{n-i}{2}}}{(n - i)!}\cdot\frac{f_{i}}{(i - 1)!}\\ &=\sum_{i=0}^{n}\frac{2^{\binom{n-i}{2}}}{(n - i)!}\cdot\frac{f_{i}}{(i - 1)!} \end{aligned}

发现这是一个卷积的形式。
设多项式 f(x),g(x),h(x)

\begin{aligned} [x^{n}] f(x) &= \frac{2^{\binom{n}{2}}}{n!}\\ [x^{n}] g(x) &= \frac{f_{n}}{(n - 1)!}\\ [x^{n}] h(x) &= \frac{2^{\binom{n}{2}}}{(n - 1)!} \end{aligned}

则有:

f(x) g(x) = h(x)

g(x) = f^{-1}(x)h(x) \pmod{x^{n + 1}}

答案即为 (n - 1)![x^{n}]g(x)

复杂度 O(n \log{n})

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
#include<cstdio>
#include<algorithm>

using i64=long long;
using uchar=unsigned char;

constexpr int mxdg(262144);
constexpr int p(1004535809);
constexpr int proot(3);

template<typename _Tp>
inline void swap(_Tp&x,_Tp&y)
{_Tp z=x;x=y;y=z;}
template<typename _Tp>
inline int fpow(_Tp v,int n){
int pw=1;
for(;n;n>>=1,v=(i64)v*v%p)
if(n&1)pw=(i64)pw*v%p;
return pw;
}
template<typename x_tp,typename y_tp>
inline void inc(x_tp&x,const y_tp&y)
{(p-y<=x)?(x+=y-p):(x+=y);}
template<typename x_tp,typename y_tp>
inline void dec(x_tp&x,const y_tp&y)
{x-=y;(x<0)&&(x+=p);}

template<typename x_tp,typename y_tp>
inline x_tp pls(const x_tp&x,const y_tp&y)
{return(x<p-y)?(x+y):(x-p+y);}
template<typename x_tp,typename y_tp>
inline x_tp sub(const x_tp&x,const y_tp&y)
{return(x<y)?(x-y+p):(x-y);}

namespace IOManager{
constexpr int FILESZ(131072);
char buf[FILESZ];
const char*ibuf=buf,*tbuf=buf;

struct IOManager{
inline char gc()
{return(ibuf==tbuf)&&(tbuf=(ibuf=buf)+fread(buf,1,FILESZ,stdin),ibuf==tbuf)?EOF:*ibuf++;}

template<typename _Tp>
inline operator _Tp(){
_Tp s=0u;char c=gc();
for(;c<48;c=gc());
for(;c>47;c=gc())
s=(_Tp)(s*10u+c-48u);
return s;
}
};
}IOManager::IOManager io;

int wn[mxdg],
iwn[mxdg];

inline void init(){
const int wnv=fpow(proot,(p-1)/mxdg);
const int iwnv=fpow(wnv,p-2);
wn[0]=iwn[0]=1;
for(int i=1,v=1,iv=1;i!=mxdg;++i){
v=(i64)v*wnv%p;wn[i]=v;
iv=(i64)iv*iwnv%p;iwn[i]=iv;
}
}

template<typename _Tp>
void DFT(_Tp*A,const int&deg){
for(int i=0,j=0;i!=deg;++i){
if(i<j)swap(A[i],A[j]);
for(int k=deg>>1;(j^=k)<k;k>>=1);
}
_Tp*const td=A+deg;
for(int l=1,tp=mxdg>>1;l<deg;l+=l,tp>>=1)
for(_Tp*i=A,z;i<td;i+=l+l)
for(int j=0,*w=wn;j!=l;++j,w+=tp)
z=(i64)i[j+l]**w%p,
i[j+l]=sub(i[j],z),
inc(i[j],z);
}
template<typename _Tp>
void IDFT(_Tp*A,const int&deg){
for(int i=0,j=0;i!=deg;++i){
if(i<j)swap(A[i],A[j]);
for(int k=deg>>1;(j^=k)<k;k>>=1);
}
_Tp*const td=A+deg;
for(int l=1,tp=mxdg>>1;l<deg;l+=l,tp>>=1)
for(_Tp*i=A,z;i<td;i+=l+l)
for(int j=0,*w=iwn;j!=l;++j,w+=tp)
z=(i64)i[j+l]**w%p,
i[j+l]=sub(i[j],z),
inc(i[j],z);

const int invdeg=fpow(deg,p-2);
for(int i=0;i!=deg;++i)
A[i]=(i64)A[i]*invdeg%p;
}

int tmp[1][mxdg+1];

inline void cp(const int*const&srcs,const int*const&srct,int*const&dsts,int*const&dstt){
std::copy(srcs,srct,dsts);
if(srct-srcs<dstt-dsts)
std::fill(dsts+(srct-srcs),dstt,0);
}

void polyinv(const int*const&F,int*const&invF,const int&n){
std::fill(invF,invF+n,0);
static int*const&inv_tmp=tmp[0];
invF[0]=fpow(F[0],p-2);
int cb=0;
for(;(1<<cb)<=n;++cb);
if(n==(1<<--cb))
--cb;
for(int deg=4;~cb;--cb,deg+=deg){
const int cn=(n+(1<<cb)-1)>>cb;
cp(F,F+cn,inv_tmp,inv_tmp+deg);

DFT(inv_tmp,deg);DFT(invF,deg);
for(int i=0;i!=deg;++i)
invF[i]=(i64)invF[i]*sub(2ll,(i64)inv_tmp[i]*invF[i]%p)%p;
IDFT(invF,deg);

std::fill(invF+cn,invF+deg,0);
}
}

int F[mxdg+1],
G[mxdg+1],
H[mxdg+1];
int inv[mxdg+1];

int main(){
const int n=io;

inv[0]=inv[1]=F[0]=F[1]=H[1]=1;
for(int i=2,pw=2;i<=n;++i,inc(pw,pw))
inv[i]=(i64)(p-p/i)*inv[p%i]%p,
H[i]=(i64)F[i-1]*pw%p,
F[i]=(i64)H[i]*inv[i]%p;

int deg=1;
for(;deg<=n+n+3;deg+=deg);
init_wn(deg);
polyinv(F,G,n+1);

DFT(G,deg);DFT(H,deg);
for(int i=0;i!=deg;++i)
G[i]=(i64)G[i]*H[i]%p;
IDFT(G,deg);

int ans=G[n];
for(int i=2;i!=n;++i)
ans=(i64)ans*i%p;
printf("%d",ans);

return 0;
}