【洛谷4841】城市规划（多项式）
题目：
分析：
一句话题意：求\(n\)个点的带标号简单（无重边无自环）无向连通图数目。（以下提到的所有图都是带标号简单无向图）
正难则反。设\(f[n]\)表示\(n\)个点的连通图数目，\(g[n]\)表示\(n\)个点的图数目。那么从\(g[n]\)中减去不连通图的数目就是\(f[n]\)。
对于两个点之间的边只有选和不选两种选择，所以显然\(g[n]=2^{C_n^2}=2^\frac{n(n-1)}{2}\)。现在问题变成了不连通图的数目。
任意钦定一个点，比如说\(1\)号点。枚举\(1\)号点所在连通块的大小\(i\)，则相当于先从\(n-1\)个点中选出\(i-1\)个点和\(i\)组成连通块，再将剩下的\(n-i\)个点组成一张图，即：
\[f[n]=g[n]-\sum_{i=1}^{n-1} C_{n-1}^{i-1}f[i]g[n-i]
\]
由于钦定了一个具体的点（\(1\)号点）所在的连通块来讨论，所以不会重复。
\(g[n]=2^{C_n^2}\)，很好算。下面来解决后面的求和怎么算。首先把组合数按定义展开：
\[\sum_{i=1}^{n-1}\frac{g[n-i](n-1)!}{(i-1)!(n-i)!}f[i]
\]
可以看出，可以把\((n-1)!\)提到前面，然后有些项和\(i\)有关，有些项和\(n-i\)有关。把它们分别搞到一块，就可以……卷积？
尝试搞到一块……
\[(n-1)!\sum_{i=1}^{n-1}\frac{f[i]}{(i-1)!}\cdot \frac{g[n-i]}{(n-i)!}
\]
右边看起来非常像卷积，于是设两个多项式：
\[F(x)=\sum_{i=0} \frac{f[i]x^i}{(i-1)!}
\]
\[G(x)=\sum_{i=0} \frac{g[n-i]x^i}{i!}
\]
至此看起来似乎可以分治FFT？博主不会，告辞
想想原来那个式子；
\[f[n]=g[n]-(n-1)!\sum_{i=1}^{n-1}\frac{f[i]}{(i-1)!}\cdot \frac{g[n-i]}{(n-i)!}
\]
把右边第二项移到左边，再两边除以\((n-1)!\)：
\[\frac{f[n]}{(n-1)!}+\sum_{i=1}^{n-1}\frac{f[i]}{(i-1)!}\cdot \frac{g[n-i]}{(n-i)!}=\frac{g[n]}{(n-1)!}
\]
发现\(\frac{f[n]}{(n-1)!}\)和求和里面关于\(f[i]\)那一项长得很像，并且当\(i=n\)时\(\frac{g[n-i]}{(n-i)!}=g[0]=1\)，所以可以直接把它扔到求和里面，即：
\[\sum_{i=1}^n\frac{f[i]}{(i-1)!}\cdot \frac{g[n-i]}{(n-i)!}=\frac{g[n]}{(n-1)!}
\]
然后设多项式\(H(x)=\sum_i \frac{2^{C_i^2}x^i}{(i-1)!}\)，则\(FG=H\)。\(G\)和\(H\)都是已知的，则\(F=G^{-1}H\)，求出\(F\)后第\(n\)项系数乘上\((n-1)!\)就是答案。
对于最终的式子还有另一种理解：\(1\)所在的连通块大小从\(1\)到\(n\)的所有情况之和就是\(n\)个点的图的方案数，即：
\[g[n]=\sum_{i=1}^n C_{n-1}^{i-1}f[i]g[n-i]
\]
同样把组合数展开，就能得出上面的式子。
代码：
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cctype>
using namespace std;
namespace zyt
{
template<typename T>
inline bool read(T &x)
{
char c;
bool f = false;
x = 0;
do
c = getchar();
while (c != EOF && c != '-' && !isdigit(c));
if (c == EOF)
return false;
if (c == '-')
f = true, c = getchar();
do
x = x * 10 + c - '0', c = getchar();
while (isdigit(c));
if (f)
x = -x;
return true;
}
template<typename T>
inline void write(T x)
{
static char buf[20];
char *pos = buf;
if (x < 0)
putchar('-'), x = -x;
do
*pos++ = x % 10 + '0';
while (x /= 10);
while (pos > buf)
putchar(*--pos);
}
typedef long long ll;
const int N = 1.3e5 + 10, LEN = N << 2, p = 1004535809, g = 3;
inline int power(int a, int b)
{
int ans = 1;
while (b)
{
if (b & 1)
ans = (ll)ans * a % p;
a = (ll)a * a % p;
b >>= 1;
}
return ans;
}
inline int inv(const int a)
{
return power(a, p - 2);
}
namespace Polynomial
{
int omega[LEN], winv[LEN], rev[LEN];
void init(const int n, const int lg2)
{
int w = power(g, (p - 1) / n), wi = inv(w);
omega[0] = winv[0] = 1;
for (int i = 1; i < n; i++)
{
omega[i] = (ll)omega[i - 1] * w % p;
winv[i] = (ll)winv[i - 1] * wi % p;
}
for (int i = 0; i < n; i++)
rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) << (lg2 - 1)));
}
void ntt(int *a, const int *w, const int n)
{
for (int i = 0; i < n; i++)
if (i < rev[i])
swap(a[i], a[rev[i]]);
for (int l = 1; l < n; l <<= 1)
for (int i = 0; i < n; i += (l << 1))
for (int k = 0; k < l; k++)
{
int tmp = (a[i + k] - (ll)w[n / (l << 1) * k] * a[i + l + k] % p + p) % p;
a[i + k] = (a[i + k] + (ll)w[n / (l << 1) * k] * a[i + l + k] % p) % p;
a[i + l + k] = tmp;
}
}
void _inv(const int *a, int *b, const int n)
{
if (n == 1)
b[0] = ::zyt::inv(a[0]);
else
{
static int tmp[LEN];
_inv(a, b, (n + 1) >> 1);
int m = 1, lg2 = 0;
while (m < (n << 1))
m <<= 1, ++lg2;
init(m, lg2);
memcpy(tmp, a, sizeof(int[n]));
memset(tmp + n, 0, sizeof(int[m - n]));
memset(b + ((n + 1) >> 1), 0, sizeof(int[m - ((n + 1) >> 1)]));
ntt(tmp, omega, m), ntt(b, omega, m);
for (int i = 0; i < m; i++)
b[i] = (2LL * b[i] % p - (ll)tmp[i] * b[i] % p * b[i] % p + p) % p;
ntt(b, winv, m);
int invm = ::zyt::inv(m);
for (int i = 0; i < n; i++)
b[i] = (ll)b[i] * invm % p;
}
}
void inv(const int *a, int *b, const int n)
{
static int tmp[N];
memcpy(tmp, a, sizeof(int[n]));
_inv(tmp, b, n);
}
void mul(const int *a, const int *b, int *c, const int n)
{
static int x[LEN], y[LEN];
memcpy(x, a, sizeof(int[n]));
memcpy(y, b, sizeof(int[n]));
int m = 1, lg2 = 0;
while (m < (n << 1))
m <<= 1, ++lg2;
memset(x + n, 0, sizeof(int[m - n]));
memset(y + n, 0, sizeof(int[m - n]));
init(m, lg2);
ntt(x, omega, m), ntt(y, omega, m);
for (int i = 0; i < m; i++)
x[i] = (ll)x[i] * y[i] % p;
ntt(x, winv, m);
int invm = ::zyt::inv(m);
for (int i = 0; i < n; i++)
c[i] = (ll)x[i] * invm % p;
}
}
int fac[N], finv[N], n;
void init()
{
fac[0] = 1;
for (int i = 1; i <= n; i++)
fac[i] = (ll)fac[i - 1] * i % p;
finv[n] = inv(fac[n]);
for (int i = n; i > 0; i--)
finv[i - 1] = (ll)finv[i] * i % p;
}
int F[LEN], G[LEN], H[LEN];
int work()
{
using Polynomial::inv;
using Polynomial::mul;
read(n);
init();
for (int i = 0; i <= n; i++)
{
G[i] = (ll)power(2, (ll)i * (i - 1) / 2 % (p - 1)) * finv[i] % p;
if (i)
H[i] = (ll)power(2, (ll)i * (i - 1) / 2 % (p - 1)) * finv[i - 1] % p;
}
inv(G, G, n + 1);
mul(G, H, F, n + 1);
write((ll)F[n] * fac[n - 1] % p);
return 0;
}
}
int main()
{
return zyt::work();
}
