【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 ivorysifreopen("f1.in","r",stdin); #endifSolve(); }轉載于:https://www.cnblogs.com/ivorysi/p/9756961.html
總結
以上是生活随笔為你收集整理的【BZOJ】3456: 城市规划(多项式求ln)的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 用MYSQL的存储过程创建百万级测试数据
- 下一篇: BZOJ2648: SJY摆棋子