闲话/社论 22.12.15 两道计数题
闲话
好今天闲话终于有点拿得出手的东西了
向计数巨佬们致以敬意和谢意。
既然代码之后补了，那我这也少写点闲话吧（
一会再补！
下面是社论部分！
两道计数题
首先感谢 jijidawang 和 Alpha1022。
以及还有人记得这题吗？22.10.05写的那道！
给定 \(n\)，\(M\)。你有 \(n\) 台电脑排成一排，你需要依次开启所有电脑。
你可以手动开启一台电脑。在任意时刻，若电脑 \(i-1\) 与电脑 \(i+1\) 都已经开启 \((1<i<n)\)，电脑 \(i\) 将立刻被自动开启。你不能再开启已经开启的电脑。
求你有多少种开启电脑的方案。两个方案不同当且仅当你手动开启的电脑的集合不同，或是手动开启电脑的顺序不同。答案对 \(M\) 取模。
我们首先枚举有 \(k\) 个位置的电脑是自动开启的，并将每一段全手动开启的电脑的方案数用 egf 表示。
容易发现答案可以被表为
\[\sum_{k=0}^n \left[\frac{x^{n-k+1}}{(n - k + 1)!}\right] \left(\sum_{i\ge 1} \frac {2^{i-1}}{i!} x^i \right)^k = \sum_{k=0}^n (n - k + 1)! [x^{n - k + 1}] \left(\frac{e^{2x} - 1}{2} \right)^k
\]
考虑设
\[f = \sqrt{\frac{x(e^{2x} - 1)}{2}}
\]
则 \(f\) 的最低次项为 \(x\)，系数为 \(1\)。对于一个 \(k\)，我们要计算的就是 \([x^{n+1}] f^{2k}\)。考虑拉格朗日反演，设 \(g = f^{\langle -1\rangle}\)，则我们要求的就是
\[[x^{n+1}] f^{2k} = \frac {2k}{n + 1} [x^{n - 2k + 1}]\left(\frac{x}{g}\right)^{n+1}
\]
如果我们能求出 \(g\)，我们就可以迅速得到每个值。
由 \(g\) 的定义得到
\[\sqrt{\frac{g(e^{2g} - 1)}{2}} = x
\]
即
\[\frac{g(e^{2g} - 1)}{2} - x^2 = 0
\]
分治 fft / 多项式牛顿迭代即可。总时间复杂度 \(O(n \log n)/O(\frac{n \log^2 n}{\log \log n})\)。
牛迭
poly newton(int n) {
poly F(2), A, B; F[1] = 1;
for (int k = 2; k < (n << 1); k <<= 1) {
F.redegree(k);
A = (2 * F).exp();
B = ginv(2) * (A - 1) + F * A;
A = ginv(2) * ((A - 1) * F) - "x^2"_p;
while(B[0] == 0) B <<= 1, A <<= 1;
F -= A * B.inv();
} return F.split(n);
}
signed main() {
cin >> n;
poly g = newton(n + 2);
g = (g << 1).inv().pow(n + 1);
for (int k = 0; n - 2 * k + 1 >= 0; ++ k) {
ans = (ans + 1ll * gfac(n - k + 1) * 2 * k % mod * ginv(n + 1) % mod * g[n - 2 * k + 1]) % mod;
} cout << ans << '\n';
}
[PA 2019 Final] Grafy
计数 \(n\) 个点的有标号无向简单图，满足每个点的入度、出度均为 \(2\)，且无重边、自环。
\(n \le 500\)。
不妨设标号为 \(0\sim n - 1\)。
考虑计数边集组成的排列。我们将每条边拆成两条无向边，随后抽象出边为长 \(2n\) 的排列 \(0\sim 2n - 1\)。\(p_u = v\) 表示 \(\left\lfloor u/2 \right\rfloor \to \left\lfloor v/2 \right\rfloor\) 有边。
一张图同一点的两条出边可以换位置，表示 \(v\) 节点的值（\(2v, 2v + 1\)）也可以换位置，这样每张图能导出 \(2^{2n}\) 个排列。
怎么算呢？考虑一手容斥。
思考什么样的排列是不满足条件的。
- 重边。
\(\left\lfloor {p_{2i}}/2 \right\rfloor = \left\lfloor p_{2i + 1}/2 \right\rfloor\) 肯定不满足。 - 自环。
\(\left\lfloor {p_{2i}}/2 \right\rfloor = i\) 或者 \(\left\lfloor p_{2i + 1}/2 \right\rfloor = i\) 都不可以。
因此答案容斥信息可以表示为（这不是很懂）：
\[\prod_{i=0}^{n-1} \left(1 - [重边] - [2i\ 自环] - [2i+1\ 自环] + 2[均自环]\right)
\]
随后可以枚举不满足条件的点的数量，即出现了 \(i\) 个双自环，\(j\) 个单自环，\(k\) 个重边。有答案
\[\small\frac{1}{2^{2n}} \sum_{i + j + k \le n} \binom{n}{i} \binom{n-i}{j} \binom{n-i-j}{k} 2^i(-1)^{j+k}\times 2^i\times 2^{2j} \times (n - i - j)^{\underline k} 2^k \times (2n - 2i - j - 2k)!
\]
四个 \(\times\) 分割的五个部分分别是：容斥系数；双自环的贡献；单自环的贡献；重边的贡献；剩下部分的贡献。
\[\frac {1}{2^{2n}}\sum_{i + j + k \le n} (-1)^{j+k} 2^{2i+2j+k} \frac{n!\times (n - i - j)! \times (n - i - j - k)!}{i!j!k!(n-i-j-k)!^2}
\]
枚举计算即可。复杂度 \(O(n^3)\)。
code
全屏放得开（
int n, ans, pw2[N];
signed main() {
cin >> n >> mod; Mod.init(mod);
init_fac(n << 1);
pw2[0] = 1; rep(i,1,n<<1) pw2[i] = add(pw2[i - 1], pw2[i - 1]);
rep(i,0,n) rep(j,0,n-i) rep(k,0,n-i-j) {
addi(ans, mul(fac[(n - i - k) * 2 - j], ifac[i], ifac[j], ifac[k], ifac[n - i - j - k], pw2[(i + j) * 2 + k], (j + k) & 1 ? mod - 1 : 1, fac[n - i - j], ifac[n - i - j - k]) );
} muli(ans, fac[n]);
muli(ans, qp(pw2[n << 1], mod - 2)); // act as a prime
cout << ans << '\n';
}
能不能再给力点啊？
\(n \le 10^5\)。
换种统计答案的方法。
我们由一张原图建一张新图，新图是一张有重边无自环的边带标号图。假设原图内有 \((i,u), (i,v)\) 两条边，那么在新图中连一条 \((u,v)\) 边，标号为 \(i\)。由于边无法退化成点，因此保证了原图无重边，这样只需要原图任意边的标号不是自己的一个端点，且没有相同标号就行了。再观察我们可以发现新图是一系列环拼合起来的。自然考察单个环。
我们仍然沿用容斥的思路，现在就需要钦定一系列边的标号是自己的一个端点。考虑令 \(G(x, y)\) 为一个环的 egf，令 \([x^n y^k / n!] G(x, y)\) 为 \(n\) 个点的环的计数，满足有 \(k\) 条边未被确定标号。
先不考虑未确定标号的边，留下来的就是一系列链了。可以发现一条链上一定有一个点满足其左侧边需要选左端点，右侧边需要选右端点。因此假如有一条 \(i\) 条边 / \(i + 1\) 个点的链，其对答案的贡献是 \(i + 1\)。然后环可以翻转所以除 \(2\)，各段也可以顺次转动所以乘 \(n / k\)。
所以这部分（\(k > 0\)）的系数就是
\[\frac {n}{2k}{{n}\brack{1}} [x^n]\left(\sum_{i>0} ix^i\right)^k = \frac{n!}{2k} [x^n]\left(\frac{x}{(1 - x)^2}\right)^k
\]
\(k = 0\) 时平凡，就是第一类斯特林数。
于是
\[\begin{aligned}
G(x, y) & \ = \sum_{n \ge 2}\sum_{k \ge 1} \frac{x^n y^k}{n!} \frac{n!}{2k} [x^n]\left(\frac{x}{(1 - x)^2}\right)^k + \sum_{n\ge 2}\frac{x^n}{n!} (n-1)!
\\ & \ = \frac 12 \sum_{n \ge 2}\sum_{k \ge 1} \frac{x^n y^k}{k} [x^n]\left(\frac{x}{(1 - x)^2}\right)^k + \sum_{n\ge 2}\frac{x^n}{n}
\\ & \ = \frac 12\left( \sum_{k \ge 1} \frac{y^k}{k}\left(\frac{x}{(1 - x)^2}\right)^k - xy\right) + \sum_{n\ge 2}\frac{x^n}{n}
\\ & \ = \frac 12\left( \sum_{k \ge 1} \frac{\left(xy / (1 - x)^2\right)^k}{k} - xy\right) + \ln \frac 1 {1 - x} - x
\\ & \ = \frac 12\left( \ln \frac 1{1 - xy / (1 - x)^2} - xy\right) + \ln \frac 1 {1 - x} - x
\end{aligned}\]
这是一个环的 egf。得到原图的 egf \(F(x,y)\) 只需要将它作 \(\exp\)，即
\[\begin{aligned}
F(x, y) & \ = \exp G(x, y)
\\ & \ = \exp\left(\frac 12\left( \ln \frac 1{1 - xy / (1 - x)^2} - xy\right) + \ln \frac 1 {1 - x} - x\right)
\\ & \ = \frac{\exp(-\frac 12 xy - x)}{(1 - x)\sqrt{1 - xy / (1 - x)^2} }
\end{aligned}\]
经典枚举标号随意的边数为 \(i\)，容斥系数为 \((-1)^{n - i}\)。答案即为
\[\begin{aligned}
& \sum_{k=0}^n (-1)^{n - k} \left[\frac{x^ny^k}{n!k!}\right] F(x, y)
\\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {(-1)^{n - k} n!k!}\right] F(x, y)
\\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {n!k!}\right] F(-x, -y)
\\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {n!k!}\right] \frac{\exp(-\frac 12 xy + x)}{(1 + x)\sqrt{1 - xy / (1 + x)^2} }
\\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {n!k!}\right] \frac{e^x}{(1 + x)} \times \exp(-xy / 2) \times \left({1 - xy / (1 + x)^2}\right) ^ {-1/2}
\\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {n!k!}\right] \frac{e^x}{(1 + x)} \times \sum_{i\ge 0}\frac {(-xy / 2)^i}{i!} \times \sum_{j\ge 0} \binom{-1/2}{j} \left(- \frac {xy}{(1 + x)^2}\right) ^ {j}
\\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {n!k!}\right] \frac{e^x}{(1 + x)} \times \sum_{i\ge 0} \sum_{j\ge 0} \frac {(-xy / 2)^i}{i!} \frac{(-1/2)^{\underline j}}{j} \left(- \frac {xy}{(1 + x)^2}\right) ^ {j}
\\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {n!k!}\right] \frac{e^x}{(1 + x)} \times \sum_{i\ge 0} \sum_{j\ge 0} y^{i+j}\frac {(-x / 2)^i}{i!}\frac{(-1/2)^{\underline j}}{j!} \left(- \frac {x}{(1 + x)^2}\right) ^ {j}
\\ =\ & \left[\frac{x^n}{n!}\right] \frac{e^x}{(1 + x)} \times \sum_{i\ge 0} \sum_{j\ge 0} (i + j)! \frac {(-x / 2)^i}{i!}\frac{(-1/2)^{\underline j}}{j!} \left(- \frac {x}{(1 + x)^2}\right) ^ {j} & (k = i + j)
\\ =\ & \left[\frac{x^n}{n!}\right] \frac{e^x}{(1 + x)} \times \sum_{j\ge 0} (-1/2)^{\underline j} \left(- \frac {x}{(1 + x)^2}\right) ^ {j} \sum_{i\ge 0} \frac{(i + j)!}{j!} \frac {(-x / 2)^i}{i!} & (提出\ i\ 相关部分)
\\ =\ & \left[\frac{x^n}{n!}\right] \frac{e^x}{(1 + x)} \times \sum_{j\ge 0} (-1/2)^{\underline j} \left(- \frac {x}{(1 + x)^2}\right) ^ {j} \sum_{i\ge 0} \binom{i+j}{i}(-x / 2)^i
\\ =\ & \left[\frac{x^n}{n!}\right] \frac{e^x}{(1 + x)} \times \sum_{j\ge 0} (-1/2)^{\underline j} \left(- \frac {x}{(1 + x)^2}\right) ^ {j} \sum_{i\ge 0} (-1)^i \binom{-j-1}{i} (-x / 2)^i & (上指标反转)
\\ =\ & \left[\frac{x^n}{n!}\right] \frac{e^x}{(1 + x)} \times \sum_{j\ge 0} (-1/2)^{\underline j} \left(- \frac {x}{(1 + x)^2}\right) ^ {j} (1 + x / 2)^{-j-1}
\\ =\ & \left[\frac{x^n}{n!}\right] \frac{e^x}{(1 + x)(1 + x / 2)} \times \sum_{j\ge 0} (-1/2)^{\underline j} \left(- \frac {x}{(1 + x)^2(1 + x / 2)}\right) ^ {j}
\end{aligned}\]
气势磅礴……来自 EI。这里的推导尽量没跳步。
最后的式子里唯一不像微分有限的是那个 \((-1/2)^{\underline j}\)。然而需要了解的是，
\[P(x) = \sum_{j\ge 0} (-1/2)^{\underline j} (-x)^j = \sum_{j\ge 0} (-1)^j (-1/2)^{\underline j} x^j = \sum_{j\ge 0} (1/2)^{\overline j} x^j = F\left( \left.\begin{aligned} 1, &1, \frac 12 \\ &1 \end{aligned} \right\rvert x \right)
\]
超几何函数形式。其满足的递推形式为
\[x^2 P'(x) = (1 - x / 2) P(x) - 1
\]
因此最后的整个式子都是微分有限的。根据经典结论，递推可做到 \(O(n) / O(\sqrt n \log n)\)。
code
#include <bits/stdc++.h>
using namespace std; using pii = pair<int,int>; using vi = vector<int>; using ll = long long;
using ull = unsigned long long; using db = double; using ld = long double; using lll = __int128_t;
mt19937 rnd(chrono::steady_clock::now().time_since_epoch().count());
template <typename T> T rand(T l, T r) { return uniform_int_distribution<T>(l, r)(rnd); }
template <typename T> void get(T & x) {
x = 0; char ch = getchar(); bool f = false; while (ch < '0' or ch > '9') f = f or ch == '-', ch = getchar();
while ('0' <= ch and ch <= '9') x = (x << 1) + (x << 3) + ch - '0', ch = getchar(); f && (x = -x);
} template <typename T, typename ... Args> void get(T & a, Args & ... b) { get(a); get(b...); }
#define iot ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
template <typename T1, typename T2> T1 min(T1 a, T2 b) { return a < b ? a : b; }
template <typename T1, typename T2> T1 max(T1 a, T2 b) { return a > b ? a : b; }
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define ufile(x)
#define rep(i,s,t) for (register int i = (s), i##_ = (t) + 1; i < i##_; ++ i)
#define pre(i,s,t) for (register int i = (s), i##_ = (t) - 1; i > i##_; -- i)
#define Aster(i, s) for (int i = head[s], v; i; i = e[i].next)
#define all(s) s.begin(), s.end()
#define eb emplace_back
#define pb pop_back
#define em emplace
const int N = 1e7 + 5;
const int inf = 0x3f3f3f3f;
const ll infll = 0x3f3f3f3f3f3f3f3fll;
int mod, inv2;
// const int mod = 10007;
// const int mod = 469762049, g = 3;
// const int mod = 998244353; // const int g = 3;
// const int mod = 1004535809, g = 3;
// const int mod = 1e9 + 7;
// const int mod = 1e9 + 9;
// const int mod = 1e9 + 3579, bse = 131;
/* add / sub */ template <typename T1, typename T2> T1 add(T1 a, T2 b) { return (a += b) >= mod ? a - mod : a; } template <typename T1, typename ...Args> T1 add(T1 a, Args ... b) { return add(a, add(b...)); } template <typename T1, typename T2> T1 sub(T1 a, T2 b) { return (a -= b) < 0 ? a + mod : a; } template <typename T1, typename ...Args> T1 sub(T1 a, Args ... b) { return add(a, add(b...)); } template <typename T1, typename T2> void addi(T1 & a, T2 b) { (a += b) >= mod ? (a -= mod) : true; } template <typename T1, typename ...Args> void addi(T1 & a, Args ...b) { addi(a, add(b...)); } template <typename T1, typename T2> void subi(T1 & a, T2 b) { (a -= b) < 0 ? (a += mod) : true; } template <typename T1, typename ...Args> void subi(T1 & a, Args ...b) { subi(a, add(b...)); }
/* Fastmod / mul */ struct FastMod { int m; ll b; void init(int _m) { m = _m; if (m == 0) m = 1; b = ((lll)1<<64) / m; } FastMod(int _m) { init(_m); } int operator() (ll a) {ll q = ((lll)a * b) >> 64; a -= q * m; if (a >= m) a -= m; return a; } } Mod(mod); int mul(int a, int b) { return Mod(1ll * a * b); } template <typename T1, typename T2> int mul(T1 a, T2 b) { return Mod((long long)(1ll * a * b)); } template <typename T, typename ...Args> int mul(T a, Args ...b) { return mul(a, mul(b...)); } template <typename T1, typename T2> void muli(T1 & a, T2 b) { a = Mod(1ll * a * b); } template <typename T1, typename ...Args> void muli(T1 & a, Args ...b) { muli(a, mul(b ...)); } // /* trivial multiple function(idk whether its useful*/ int mul(int a, int b) { return 1ll * a * b % mod; } template <typename T1, typename T2> int mul(T1 a, T2 b) { return (long long)(1ll * a * b) % mod; } template <typename T, typename ...Args> int mul(T a, Args ...b) { return mul(a, mul(b...)); }
int n, m, fac[N], inv[N], f[N], g[N], t1, ret;
int qp(int a, int b) {
int ret = 1;
while (b) {
if (b & 1) ret = mul(ret, a);
a = mul(a, a);
b >>= 1;
} return ret;
}
vi operator*(vi a, vi b) {
static long long tem[100001];
int size = a.size() + b.size() - 1;
vi c(size);
memset(tem, 0, size << 3);
for (int i = 0; i < a.size(); i++) for (int j = 0; j < b.size(); j++)
tem[i + j] = add(tem[i + j], mul(a[i], b[j]));
for (int i = 0; i < size; i++) c[i] = tem[i] % mod;
return c;
}
vi operator-(vi a, vi b) {
if (a.size() < b.size()) a.resize(b.size());
for (int i = 0; i < b.size(); i++) a[i] = sub(a[i], b[i]);
return a;
}
vi deri(vi a) {
for (int i = 0; i < a.size() - 1; i++) a[i] = mul(a[i + 1], i + 1);
a.pop_back();
return a;
}
vi shift(vi a, int n) {
a.resize(a.size() + n);
for (int i = a.size() - 1; i >= n; i--) a[i] = a[i - n];
rep(i,0,n-1) a[i] = 0;
return a;
}
void div(int a[], int v) {
rep(i,1,n) a[i] = sub(a[i], mul(v, a[i - 1]));
}
int C(int n, int m) {
if (n < m or m < 0) return 0;
return mul(fac[n], inv[m], inv[n - m]);
}
int main() {
get(n, mod); Mod.init(mod); inv2 = mod + 1 >> 1;
fac[0] = fac[1] = inv[0] = inv[1] = 1;
rep(i,2,n) fac[i] = mul(fac[i - 1], i), inv[i] = mul(mod - mod / i, inv[mod % i]);
rep(i,2,n) inv[i] = mul(inv[i - 1], inv[i]);
rep(i,0,n) g[i] = inv[i];
div(g, 1), div(g, inv2);
vi H = {1, 2 + inv2, 2, inv2}, a = H, b = vi() - H, c = (H - shift(deri(H), 1));
a = shift(a, 2), b[1] = add(b[1], inv2);
b = b * c, c = c * H;
f[0] = 1;
rep(i,1,n) {
t1 = 0;
for (int j = 0; j < a.size(); j++) if (i + 1 - j >= 0)
t1 = add(t1, mul(i + 1 - j, f[i + 1 - j], a[j]));
for (int j = 0; j < b.size(); j++) if (i - j >= 0)
t1 = add(t1, mul(f[i - j], b[j]));
if (c.size() > i) t1 = add(t1, c[i]);
f[i] = mul(mod - qp(b[0], mod - 2), t1);
}
rep(i,0,n) ret = add(ret, mul(f[i], g[n - i]));
printf("%d\n", mul(ret, fac[n]));
// printf("%.6lf\n", 1. * clock() / CLOCKS_PER_SEC);
}
以下是博客签名，与正文无关。
请按如下方式引用此页：
本文作者 joke3579，原文链接：https://www.cnblogs.com/joke3579/p/editorchat221215.html。
遵循 CC BY-NC-SA 4.0 协议。
请读者尽量不要在评论区发布与博客内文完全无关的评论，视情况可能删除。
