【BZOJ】3456: 城市规划（多项式求ln）
题解
在我写过分治NTT，多项式求逆之后
我又一次写了多项式求ln
我们定义一个数列的指数型生成函数为
\(\sum_{i = 0}^{n} \frac{A_{i}}{i!} x^{i}\)
然后这个有个很好的性质，是什么呢，就是我们考虑两个排列\(A\)和\(B\)，不改变原来的顺序，把它们合并成一个排列，方案数显然是
\(\binom{|A| + |B|}{|A|}\)
现在每个相同长度的排列\(A\)带有一个价值\(A_i\)，\(B\)同理
\(C_{k} = \sum_{i = 0}^{k} \binom{k}{i}A_{i}B_{k - i}\)
\(\frac{C_{k}}{k!} = \sum_{i = 0}^{k} \frac{A_{i}}{i!} \cdot \frac{B_{k - i}}{(k - i)!}\)
这样的话两个指数型函数的卷积就是我们想要的排列总和除以\(n!\)我们可以卷积之后还原回去
然后我们考虑
\(G(x) = \sum_{i = 0}^{+\infty} \frac{2^{\binom{i}{2}}}{i!}\)表示\(i\)个点带标号的无向图个数
和\(C(x) = \sum_{i = 0}^{+\infty} \frac{c_i}{i!} x^{i}\)表示\(i\)个点带标号的无向联通图个数
然后我们可以列出来
\(G(x) = \frac{C(x)}{1!} + \frac{C^2(x)}{2!} + \frac{C^3(x)}{3!} + .... \frac{C^{n}(x)}{n!} = e^{C(x)}\)
怎么理解呢，以8个点的带标号无向图举个例子
8 可以由3 5拼成，但是相乘的时候5 3会再算一次，所以除上\(2!\)
而由4 4 拼成，虽然4 4只算一次，但是1 2 3 4 5 6 7 8 和 5 6 7 8 1 2 3 4本质上一样，所以也要除上\(2!\)
然后我们可以得到
\(C(x) = ln(G(x))\)
怎么求\(ln(G(x))\)呢
设\(F(x) = ln(G(x))\)
两边分别求导
\(F'(x) = \frac{G'(x)}{G(x)}\)
然后再两边积分起来
\(F(x) = \int \frac{G'(x)}{G(x)}\)
求逆元是\(O(n \log n)\)求导\(O(n)\)求积分\(O(n)\)
代码
#include <bits/stdc++.h>
#define fi first
#define se second
#define pii pair<int,int>
#define mp make_pair
#define pb push_back
#define enter putchar('\n')
#define space putchar(' ')
#define MAXN 130005
#define mo 994711
//#define ivorysi
using namespace std;
typedef unsigned long long int64;
typedef long double db;
typedef unsigned int u32;
template<class T>
void read(T &res) {
res = 0;char c = getchar();T f = 1;
while(c < '0' || c > '9') {
if(c == '-') f = -1;
c = getchar();
}
while(c >= '0' && c <= '9') {
res = res * 10 + c - '0';
c = getchar();
}
res *= f;
}
template<class T>
void out(T x) {
if(x < 0) {putchar('-');x = -x;}
if(x >= 10) out(x / 10);
putchar('0' + x % 10);
}
const int MOD = 1004535809,MAXL = 1 << 18;
int W[(1 << 19) + 5],fac[MAXL + 5],inv[MAXL + 5],invfac[MAXL + 5],N;
int mul(int a,int b) {
return 1LL * a * b % MOD;
}
int inc(int a,int b) {
return a + b >= MOD ? a + b - MOD : a + b;
}
int fpow(int x,int c) {
int res = 1,t = x;
while(c) {
if(c & 1) res = mul(res,t);
t = mul(t,t);
c >>= 1;
}
return res;
}
struct Poly {
vector<int> p;
Poly() {p.clear();}
friend void NTT(Poly &f,int len,int on) {
f.p.resize(len);
for(int i = 1 , j = len >> 1; i < len - 1; ++i) {
if(i < j) swap(f.p[i],f.p[j]);
int k = len >> 1;
while(j >= k) {
j -= k;
k >>= 1;
}
j += k;
}
for(int h = 2 ; h <= len ; h <<= 1) {
int wn = W[(MAXL + on * MAXL / h) % MAXL];
for(int k = 0 ; k < len ; k += h) {
int w = 1;
for(int j = k ; j < k + h / 2 ; ++j) {
int u = f.p[j],t = mul(w,f.p[j + h / 2]);
f.p[j] = inc(u,t);
f.p[j + h / 2] = inc(u,MOD - t);
w = mul(w,wn);
}
}
}
if(on == -1) {
int InvL = fpow(len,MOD - 2);
for(int i = 0 ; i < len ; ++i) f.p[i] = mul(f.p[i],InvL);
}
}
friend Poly operator * (Poly a,Poly b) {
int L = a.p.size() + b.p.size() - 2;
int t = 1;
while(t <= L) t <<= 1;
a.p.resize(t);b.p.resize(t);
NTT(a,t,1);NTT(b,t,1);
Poly c;
for(int i = 0 ; i < t ; ++i) {
c.p.pb(mul(a.p[i],b.p[i]));
}
NTT(c,t,-1);
int s = c.p.size() - 1;
while(s >= 0 && c.p[s] == 0) {c.p.pop_back();--s;}
return c;
}
friend Poly operator - (Poly a,Poly b) {
Poly c;
int L = max(a.p.size(),b.p.size());
a.p.resize(L);b.p.resize(L);
for(int i = 0 ; i < L ; ++i) c.p.pb(inc(a.p[i],MOD - b.p[i]));
return c;
}
friend Poly operator + (Poly a,Poly b) {
Poly c;
int L = max(a.p.size(),b.p.size());
a.p.resize(L);b.p.resize(L);
for(int i = 0 ; i < L ; ++i) c.p.pb(inc(a.p[i],b.p[i]));
return c;
}
friend Poly Inverse(Poly f,int t) {
f.p.resize(t);
if(t == 1) {
Poly g;g.p.pb(fpow(f.p[0],MOD - 2));
return g;
}
Poly g = Inverse(f,t >> 1);
t *= 2;
NTT(f,t,1);NTT(g,t,1);
Poly r;
for(int i = 0 ; i < t; ++i) {
r.p.pb(inc(mul(2,g.p[i]),MOD - mul(mul(g.p[i],g.p[i]),f.p[i])));
}
NTT(r,t,-1);
t >>= 1;
r.p.resize(t);
--t;
while(t >= 0 && r.p[t] == 0) {r.p.pop_back();--t;}
return r;
}
friend Poly Integral(const Poly &f) {
Poly g;
int L = f.p.size();
g.p.resize(L + 1);
for(int i = 1 ; i <= L ; ++i) {
g.p[i] = mul(f.p[i - 1],inv[i]);
}
return g;
}
friend Poly Derivative(const Poly &f) {
Poly g;
int L = f.p.size();
g.p.resize(L - 1);
for(int i = 0 ; i < L - 1; ++i) {
g.p[i] = mul((i + 1),f.p[i + 1]);
}
return g;
}
friend Poly ln(const Poly &f) {
int t = 1;
while(t <= f.p.size() - 1) t <<= 1;
return Integral(Derivative(f) * Inverse(f,t));
}
}g,f;
void Solve() {
read(N);
W[0] = 1;
W[1] = fpow(3,(MOD - 1) / MAXL);
for(int i = 2 ; i < MAXL ; ++i) W[i] = mul(W[i - 1],W[1]);
fac[0] = 1;invfac[0] = 1;
inv[1] = 1;
for(int i = 2 ; i <= MAXL ; ++i) inv[i] = mul(inv[MOD % i],MOD - MOD / i);
for(int i = 1 ; i <= MAXL ; ++i) fac[i] = mul(fac[i - 1],i);
for(int i = 1 ; i <= MAXL ; ++i) invfac[i] = mul(invfac[i - 1],inv[i]);
g.p.resize(N + 1);
g.p[0] = g.p[1] = 1;
for(int i = 2 ; i <= N ; ++i) {
g.p[i] = mul(fpow(2,1LL * i * (i - 1) / 2 % (MOD - 1)),invfac[i]);
}
f = ln(g);
out(mul(fac[N],f.p[N]));enter;
}
int main() {
#ifdef ivorysi
freopen("f1.in","r",stdin);
#endif
Solve();
}
