【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();

}

巴特西

【BZOJ】3456: 城市规划（多项式求ln）

题解

代码

最新文章

热门文章