【LOJ6569】仙人掌计数

【题目链接】

点击打开链接

【思路要点】

考虑对有根仙人掌的圆方树计数，定义子树大小为子树内圆点的个数。

令子树大小为 i i i 的圆点和方点各有 r i , s i r_i,s_i ri,si 个，则其指数型生成函数分别为
R ( x ) = ∑ i = 1 r i i ! x i , S ( x ) = ∑ i = 1 s i i ! x i R(x)=\sum_{i=1}\frac{r_i}{i!}x^i,S(x)=\sum_{i=1}\frac{s_i}{i!}x^i R(x)=i=1∑i!rixi,S(x)=i=1∑i!sixi

稍加推导可得
R ( x ) = x e S ( x ) , S ( x ) = R ( x ) + R ( x ) 2 2 + R ( x ) 3 2 + ⋯ = R ( x ) + R ( x ) 2 2 − 2 R ( x ) R(x)=xe^{S(x)},S(x)=R(x)+\frac{R(x)^2}{2}+\frac{R(x)^3}{2}+\dots=R(x)+\frac{R(x)^2}{2-2R(x)} R(x)=xeS(x),S(x)=R(x)+2R(x)2+2R(x)3+⋯=R(x)+2−2R(x)R(x)2

因此有
R ( x ) = x e R ( x ) + R ( x ) 2 2 − 2 R ( x ) R(x)=xe^{R(x)+\frac{R(x)^2}{2-2R(x)}} R(x)=xeR(x)+2−2R(x)R(x)2
x e R ( x ) + R ( x ) 2 2 − 2 R ( x ) − R ( x ) = 0 xe^{R(x)+\frac{R(x)^2}{2-2R(x)}}-R(x)=0 xeR(x)+2−2R(x)R(x)2−R(x)=0

令 G ( R ( x ) ) = x e R ( x ) + R ( x ) 2 2 − 2 R ( x ) − R ( x ) G(R(x))=xe^{R(x)+\frac{R(x)^2}{2-2R(x)}}-R(x) G(R(x))=xeR(x)+2−2R(x)R(x)2−R(x) ，则 G ( R ( x ) ) = 0 G(R(x))=0 G(R(x))=0

考虑牛顿迭代解上述方程，令
R ( x ) ≡ R 0 ( x ) ( m o d x n ) R(x)\equiv R_0(x)\quad(mod\ x^n) R(x)≡R0(x)(mod xn)

则有 G ( R ( x ) ) ≡ G ( R 0 ( x ) ) + ( R ( x ) − R 0 ( x ) ) G ′ ( R 0 ( x ) ) ≡ 0 ( m o d x 2 n ) G(R(x))\equiv G(R_0(x))+(R(x)-R_0(x))G'(R_0(x))\equiv0\quad(mod\ x^{2n}) G(R(x))≡G(R0(x))+(R(x)−R0(x))G′(R0(x))≡0(mod x2n)

因此 R ( x ) ≡ R 0 ( x ) − G ( R 0 ( x ) ) G ′ ( R 0 ( x ) ) ( m o d x 2 n ) R(x)\equiv R_0(x)-\frac{G(R_0(x))}{G'(R_0(x))}\quad(mod\ x^{2n}) R(x)≡R0(x)−G′(R0(x))G(R0(x))(mod x2n)

其中 G ′ ( R ( x ) ) = d G ( R ( x ) ) d R ( x ) = x e R ( x ) + R ( x ) 2 2 − 2 R ( x ) ( 1 + 4 R ( x ) − 2 R ( x ) 2 ( 2 − 2 R ( x ) ) 2 ) − 1 G'(R(x))=\frac{dG(R(x))}{dR(x)}=xe^{R(x)+\frac{R(x)^2}{2-2R(x)}}(1+\frac{4R(x)-2R(x)^2}{(2-2R(x))^2})-1 G′(R(x))=dR(x)dG(R(x))=xeR(x)+2−2R(x)R(x)2(1+(2−2R(x))24R(x)−2R(x)2)−1

按此迭代，时间复杂度 O ( N L o g N ) O(NLogN) O(NLogN) 。

注意到我们计算的是有根仙人掌的个数，因此在最后应当将答案除去 N N N 。

【代码】

#include<bits/stdc++.h>
using namespace std;
const int MAXN = 262144;
const int P = 998244353;
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
template <typename T> void chkmax(T &x, T y) {x = max(x, y); }
template <typename T> void chkmin(T &x, T y) {x = min(x, y); }
template <typename T> void read(T &x) {x = 0; int f = 1;char c = getchar();for (; !isdigit(c); c = getchar()) if (c == '-') f = -f;for (; isdigit(c); c = getchar()) x = x * 10 + c - '0';x *= f;
}
template <typename T> void write(T x) {if (x < 0) x = -x, putchar('-');if (x > 9) write(x / 10);putchar(x % 10 + '0');
}
template <typename T> void writeln(T x) {write(x);puts("");
}
namespace Poly {const int MAXN = 262144;const int P = 998244353;const int LOG = 25;const int G = 3;int power(int x, int y) {if (y == 0) return 1;int tmp = power(x, y / 2);if (y % 2 == 0) return 1ll * tmp * tmp % P;else return 1ll * tmp * tmp % P * x % P;}int invn[MAXN], tmpa[MAXN], tmpb[MAXN];int N, Log, home[MAXN]; bool initialized;int forward[MAXN], bckward[MAXN], inv[LOG];void init() {initialized = true;forward[0] = bckward[0] = inv[0] = invn[1] = 1;for (int len = 2, lg = 1; len <= MAXN; len <<= 1, lg++)inv[lg] = power(len, P - 2);for (int i = 2; i < MAXN; i++)invn[i] = (P - 1ll * (P / i) * invn[P % i] % P) % P;int delta = power(G, (P - 1) / MAXN);for (int i = 1; i < MAXN; i++)forward[i] = bckward[MAXN - i] = 1ll * forward[i - 1] * delta % P;}void NTTinit() {for (int i = 0; i < N; i++) {int ans = 0, tmp = i;for (int j = 1; j <= Log; j++) {ans <<= 1;ans += tmp & 1;tmp >>= 1;}home[i] = ans;}}void NTT(int *a, int mode) {assert(initialized);for (int i = 0; i < N; i++)if (home[i] < i) swap(a[i], a[home[i]]);int *g;if (mode == 1) g = forward;else g = bckward;for (int len = 2, lg = 1; len <= N; len <<= 1, lg++) {for (int i = 0; i < N; i += len) {for (int j = i, k = i + len / 2; k < i + len; j++, k++) {int tmp = a[j];int tnp = 1ll * a[k] * g[MAXN / len * (j - i)] % P;a[j] = (tmp + tnp > P) ? (tmp + tnp - P) : (tmp + tnp);a[k] = (tmp - tnp < 0) ? (tmp - tnp + P) : (tmp - tnp);}}}if (mode == -1) {for (int i = 0; i < N; i++)a[i] = 1ll * a[i] * inv[Log] % P;}}void times(vector <int> &a, vector <int> &b, vector <int> &c) {assert(a.size() >= 1), assert(b.size() >= 1);int goal = a.size() + b.size() - 1;N = 1, Log = 0;while (N < goal) {N <<= 1;Log++;}for (unsigned i = 0; i < a.size(); i++)tmpa[i] = a[i];for (int i = a.size(); i < N; i++)tmpa[i] = 0;for (unsigned i = 0; i < b.size(); i++)tmpb[i] = b[i];for (int i = b.size(); i < N; i++)tmpb[i] = 0;NTTinit();NTT(tmpa, 1);NTT(tmpb, 1);for (int i = 0; i < N; i++)tmpa[i] = 1ll * tmpa[i] * tmpb[i] % P;NTT(tmpa, -1);c.resize(goal);for (int i = 0; i < goal; i++)c[i] = tmpa[i];}void timesabb(vector <int> &a, vector <int> &b, vector <int> &c) {assert(a.size() >= 1), assert(b.size() >= 1);int goal = a.size() + b.size() * 2 - 2;N = 1, Log = 0;while (N < goal) {N <<= 1;Log++;}for (unsigned i = 0; i < a.size(); i++)tmpa[i] = a[i];for (int i = a.size(); i < N; i++)tmpa[i] = 0;for (unsigned i = 0; i < b.size(); i++)tmpb[i] = b[i];for (int i = b.size(); i < N; i++)tmpb[i] = 0;NTTinit();NTT(tmpa, 1);NTT(tmpb, 1);for (int i = 0; i < N; i++)tmpa[i] = 1ll * tmpa[i] * tmpb[i] % P * tmpb[i] % P;NTT(tmpa, -1);c.resize(goal);for (int i = 0; i < goal; i++)c[i] = tmpa[i];}void getinv(vector <int> &a, vector <int> &b) {assert(a.size() >= 1), assert(a[0] != 0);b.clear(), b.push_back(power(a[0], P - 2));while (b.size() < a.size()) {vector <int> c, ta = a;ta.resize(b.size() * 2);timesabb(ta, b, c);b.resize(b.size() * 2);for (unsigned i = 0; i < b.size(); i++)b[i] = (2ll * b[i] - c[i] + P) % P;}b.resize(a.size());}void getder(vector <int> &a, vector <int> &b) {assert(a.size() >= 1);if (a.size() == 1) {b.clear();b.resize(1);} else {b.resize(a.size() - 1);for (unsigned i = 0; i < b.size(); i++)b[i] = (i + 1ll) * a[i + 1] % P;}}void getint(vector <int> &a, vector <int> &b) {b.resize(a.size() + 1), b[0] = 0;for (unsigned i = 0; i < a.size(); i++)b[i + 1] = 1ll * invn[i + 1] * a[i] % P;}void getlog(vector <int> &a, vector <int> &b) {assert(a.size() >= 1), assert(a[0] == 1);vector <int> da, inva, db;getder(a, da), getinv(a, inva);times(da, inva, db), getint(db, b);b.resize(a.size());}void getexp(vector <int> &a, vector <int> &b) {assert(a.size() >= 1), assert(a[0] == 0);b.clear(), b.push_back(1);while (b.size() < a.size()) {vector <int> lnb, res;b.resize(b.size() * 2), getlog(b, lnb);for (unsigned i = 0; i < lnb.size(); i++)if (i == 0) lnb[i] = (P + 1 + a[i] - lnb[i]) % P;else if (i < a.size()) lnb[i] = (P + a[i] - lnb[i]) % P;else lnb[i] = (P - lnb[i]) % P;times(lnb, b, res);res.resize(b.size());swap(res, b);}b.resize(a.size());}void getshl(vector <int> &a, vector <int> &b, ull bits) {if (a.size() < bits) bits = a.size();b.clear(), b.resize(bits);for (unsigned i = 0; b.size() < a.size(); i++)b.push_back(a[i]);}void getshr(vector <int> &a, vector <int> &b, ull bits) {if (a.size() < bits) bits = a.size(); b.clear();for (unsigned i = bits; i < a.size(); i++)b.push_back(a[i]);b.resize(a.size());}void getpowk(vector <int> &a, vector <int> &b, int k) {assert(k >= 1);unsigned pos = a.size();for (unsigned i = 0; i < a.size(); i++)if (a[i]) {pos = i;break;}if (pos == a.size()) {b = a;return;}int val = power(a[pos], k), inv = power(a[pos], P - 2);vector <int> lntmp, tmp;getshr(a, tmp, pos);for (unsigned i = 0; i < tmp.size(); i++)tmp[i] = 1ll * tmp[i] * inv % P;getlog(tmp, lntmp);for (unsigned i = 0; i < lntmp.size(); i++)lntmp[i] = 1ll * lntmp[i] * k % P;getexp(lntmp, tmp);for (unsigned i = 0; i < tmp.size(); i++)tmp[i] = 1ll * tmp[i] * val % P;getshl(tmp, b, 1ull * pos * k);}int getinvfunc(vector <int> &f, unsigned n) {assert(f[0] == 0 && f[1] != 0);assert(n >= 1 && n <= f.size());int inv = power(n, P - 2);vector <int> tmp;getshr(f, tmp, 1);vector <int> invf;getinv(tmp, invf);vector <int> res;getpowk(invf, res, n);return 1ll * inv * res[n - 1] % P;}
}
int fac[MAXN], inv[MAXN];
int power(int x, int y) {if (y == 0) return 1;int tmp = power(x, y / 2);if (y % 2 == 0) return 1ll * tmp * tmp % P;else return 1ll * tmp * tmp % P * x % P;
}
void init(int n) {Poly :: init();fac[0] = 1;for (int i = 1; i <= n; i++)fac[i] = 1ll * fac[i - 1] * i % P;inv[n] = power(fac[n], P - 2);for (int i = n - 1; i >= 0; i--)inv[i] = inv[i + 1] * (i + 1ll) % P;
}
void update(int &x, int y) {x += y;if (x >= P) x -= P;
}
int main() {init(131072);vector <int> r;r.push_back(0);r.push_back(1);while (r.size() < 131072) {vector <int> d, f, invf, l, e;Poly :: times(r, r, d);r.resize(r.size() * 2);d.resize(r.size());f.resize(r.size());for (unsigned i = 0; i < r.size(); i++) {d[i] = (2ll * r[i] - d[i] + P) % P;f[i] = 1ll * (P - 2) * r[i] % P;if (i == 0) update(f[i], 2);}Poly :: getinv(f, invf);Poly :: times(d, invf, l);l.resize(r.size());Poly :: getexp(l, e);vector <int> g, dg, invdg, res;Poly :: getshl(e, g, 1);for (unsigned i = 0; i < r.size(); i++)d[i] = 2ll * d[i] % P;Poly :: times(invf, invf, f);f.resize(r.size());Poly :: times(d, f, e);e.resize(r.size());update(e[0], 1);Poly :: times(g, e, dg);dg.resize(r.size());update(dg[0], P - 1);Poly :: getinv(dg, invdg);for (unsigned i = 0; i < r.size(); i++)update(g[i], P - r[i]);Poly :: times(g, invdg, res);res.resize(r.size());for (unsigned i = 0; i < r.size(); i++)update(r[i], P - res[i]);}int T; read(T);while (T--) {int x; read(x);writeln(1ll * fac[x - 1] * r[x] % P);}return 0;
}

【LOJ6569】仙人掌计数相关推荐

无标号有根仙人掌计数
下面的做法(是个常数有点大的O(nlog⁡2n)\mathcal O(n\log^2 n)O(nlog2n))我只是口胡了一下,没写过,所以可能有错首先,假设ana_nan表示有nnn个点的无标号 ...
2019hdu多校第一场H Desert [多项式计数] [仙人掌计数]
首先注意到n个点的仙人掌其实和n个圆点的圆方树是等价的,那么转化为求圆方树的统计会简单一点. 圆方树要满足的性质有以下几条: 1.奇数层为圆点,偶数层为方点,即只存在圆点和方点连的边. 2.圆点下的子 ...
[学习笔记]多项式与有标号简单图计数
学了一天的有标号无向图计数真的自闭了- 本篇文章是基于2019WC汪乐平大佬的讲课课件<生成函数,多项式算法与图的计数>编写的. 注意:文中所有生成函数都规定为指数型生成函数(EGF),请 ...
[TravelNotes] WC 2017 THUWC 2017 游记
Pre() 春节在家赶紧临时抱佛脚一发学习杜教筛和SAM 刷刷模板题水水水还想起自己数据结构的姿势水平不高啊恶补了一下替罪羊树线段树分裂二进制分组然后重新学习了LCT的姿势再看看自 ...
[日常] NOIWC2019 冬眠记
NOIWC 2019 冬眠记辣鸡rvalue天天写意识流流水账 Day 0 早上没有跑操(极度舒服.png) 和春哥在博客颓图的时候突然被来送笔电的老爹查水表(捂脸) 母上大人骗我说这功能机不能放存 ...
WXHRound#14被虐记
T2:无标号有根仙人掌计数,不会倒是搞懂了O(n^2log n)无标号无根树计数先考虑无标号有根树的计数记dp[k]为当我dp到i时用1~i大小的树可以凑出k的方案数则每次就拿dp[i]去更新 ...
●洛谷P3687 [ZJOI2017]仙人掌
题链: https://www.luogu.org/problemnew/show/P3687 题解: 计数DP,树形DP. (首先对于这个图来说,如果初始就不是仙人掌,那么就直接输出0) 然后由于本 ...
[BZOJ4784][UOJ290][ZJOI017]仙人掌-树形DP
仙人掌 Description 如果一个无自环无重边无向连通图的任意一条边最多属于一个简单环,我们就称之为仙人掌.所谓简单环即不经过重复的结点的环. 现在九条可怜手上有一张无自环无重边的无向连通图,但 ...
业界毒瘤仙人掌一条龙服务
本slide是为了NJU集训队准备....未完待续... 正经定义 : 无向图中的每条边至多位于一个简单环上,且任意两点可达. 由此可知仙人掌的构造方式很"优美",即生成一棵树,把 ...

【LOJ6569】仙人掌计数

【LOJ6569】仙人掌计数相关推荐

最新文章

热门文章