随机的艺术

随机的艺术

引子

随机,万恶之源。你做的题的数据多半是随机造的,也有能被随机乱搞过去的题,有的题标程就是随机,有的题就是出题人随机想到的 idea… 不得不说,OI 和随机还真扯不开关系。现在不妨让我们一起探索其不为大多数人所知的一面。

定义

这里先给出一些下文会用到的定义:

均匀随机:一个随机变量 $x$ 是均匀随机的当且仅当对于任何可能的取值 $v\in V$ 都有 $P(x=v)=\frac1{|V|}$。即每种情况出现的概率均等。

单射:一个映射 $f(x):X\rightarrow Y$ 是单射当且仅当不存在 $x,y\in X$ 使得 $x\ne y$ 且 $f(x)=f(y)$。

满射:一个映射 $f(x):X\rightarrow Y$ 是满射当且仅当对于任意 $y\in Y$ 都有 $x\in X$ 使得 $f(x)=y$。

双射:一个映射 $f(x):X\rightarrow Y$ 是双射当且仅当它既是单射又是满射。实际上简单地说就是 $X$ 中每一个元素和 $Y$ 中每一个元素一一对应。

主体

随机数

在开始这场旅途之前,你或许需要一把引导你进入随机世界的钥匙:随机数。有了随机的整数,我们几乎可以生成任何随机目标。随机数的生成在 这篇洛谷日报 里面已经有过一些阐释,这里仅做一个简要的总结。

朴素地,你可以用 srand(seed) 来初始化伪随机数引擎,并使用 rand() 来生成一个值域为 [0,RAND_MAX] 的随机整数。但可惜的是,标准只要求了 RAND_MAX 至少为 32767,于是乎生成的整数有时并不能满足我们的需求。而且标准对 rand() 的实现没有具体要求,因此生产的循环数质量可能并不是很高(例如,有一定循环节)。因此 C++ 来救世了!

C++11(C++11!C++11!C++11!)提供了 random 库,可以在 cppreference 上查看详细。简要来说,它提供了许多可供选择的随机数生成引擎,其中最常用的应该是 std::mt19937 梅森旋转算法,它生成的随机数达到了 $(2^{19937}-1)$ 的周期,并且分布均匀,性能较好。下面是一个使用 std::mt19937 的简单例子:

1
2
3
4
5
6
7
8
9
10
11
12
13
#include <iostream>
#include <random>

std::mt19937 engine(19260817);
int main() {
// 一个均匀随机的整数分布,使用方式为 dist(engine)
std::uniform_int_distribution<int> dist(0, 4);
int cnt[5] = { 0, 0, 0, 0, 0 };
for (int i = 0; i < 100000; ++i)
++cnt[dist(engine)];
for (int i = 0; i < 5; ++i)
std::cout << i << ": " << cnt[i] << std::endl;
}

下面是在我电脑上运行一次的输出:

1
2
3
4
5
0: 20024
1: 19821
2: 20078
3: 20133
4: 19944

可以看到分布是比较均匀的。但值得注意的是,mt19937 并不安全,也就是说攻击者可以根据 mt19937 的几个输出来推测出其内部状态,进而预测它的下一个输出。还可以参见我出的 这道题,大意为根据 mt19937 生成的随机数来反推初始化 mt19937 使用的种子。不过安全性似乎和 OIer 无关,我们就直接跳过吧。

准备工作

C++11 的随机库好是好,但调用并不是很方便,这里先定义一些在下文会用到的辅助函数:

1
2
3
4
5
6
7
8
9
10
#include <chrono>
#include <random>
#include <type_traits>

// 使用当前时间初始化 mt19937
std::mt19937 engine(std::chrono::steady_clock::now().time_since_epoch().count());
template<class T, class = std::enable_if_t<std::is_integral_v<T>>> // 要求 T 必须为整型
T rand(T l, T r) {
return std::uniform_int_distribution(l, r)(engine);
}

洗牌!

假设我们有一个数组,现在我们想要随机打乱它,并要求最后每一种排列出现的可能性相等,该怎么操作呢?

我一说,有些人啪就站起来了,说可以用 random_shuffle,很快啊!但是 random_shuffle 内部用的是 rand(),依然会有上面谈及的问题;C++11 推荐使用的是 std::shuffle(begin, end, engine),使用随机数引擎来打乱数组。不过我们今天并不会浮于表面,我们来想想如何实现。

首先你有一个绝对能保证随机的算法:枚举出 $n!$ 种不同的排列,然后随机选择一种,但是这样的时间和空间复杂度都是 $O(n!)$ 的。有人会想到用康托展开,但是当 $n!$ 超出最大的整型范围时我们就无法储存了。于是我们引入一个洗牌算法:Fisher-Yates 洗牌法,又名 Knuth 洗牌法,最初由 Ronald Fisher 和 Frank Yates 在 1938 年提出,并被 Knuth 老头子在他的 The Art of Computer Programming 中提到并传播。它的主要思想只有下面两行:

1
2
for (int i = n - 1; i; --i)
std::swap(arr[i], arr[rand(0, i)]);

没错,就这么简单。现在我们来证明一下为什么这个东西是均匀随机的。

我们考虑归纳证明。$i=1$ 的情况很显然,毕竟只有一种排列。然后我们上面的算法可以写成如下的递归形式:

1
2
3
4
5
void shuffle(int len) {
if (len == 1) return;
std::swap(arr[len - 1], arr[rand(0, len - 1)]);
shuffle(len - 1);
}

然后我们现在假设这个算法能对长度为 $n$ 的数组均匀随机打乱:也就是说,对于任意长度为 $n$ 的排列 $\pi$,我们知道 $\mathbf{P}(\textrm{arr}{[1,n]}=\pi)=\frac1{n!}$。现在我们要证明它能对长度为 $n+1$ 的数组均匀随机打乱,即对于任意长为 $n+1$ 的排列 $\pi’$,有 $\textbf{P}(\textrm{arr}’{[1,n+1]}=\pi’)=\frac1{(n+1)!}$。不妨令 $\pi’’$ 为将 $n+1$ 从 $\pi’$ 中删去后得到的排列,于是我们有 $\textbf{P}(\textrm{arr’}{[1,n+1]}=\pi’)=\textbf{P}(\textrm{arr}{[1,n]}=\pi’’)\times\frac1{n+1}=\frac1{n!}\times\frac1{n+1}=\frac1{(n+1)!}$,由此我们得证。

选择!

又是一个常见的问题:你现在要从 $[1,n]$ 中均匀随机选出 $m$ 个 不重复 的数,其中 $1\le n\le 10^{18}$,$1\le m\le \min\{n,10^6\}$。这时我们上面时间复杂度为 $O(n)$ 的算法纷纷败下阵来。我们有没有时间复杂度更优的做法呢?

实际上是有的。我们把上面那个洗牌算法下标反转一下,得到下面的算法:

1
2
for (int i = 0; i < n; ++i)
std::swap(arr[i], arr[rand(i, n - 1)]);

然后我们发现,由于这个算法在 $i\ge m$ 时对前 $m$ 个元素没有影响,因此就可以转化为:

1
2
for (int i = 0; i < m; ++i)
std::swap(arr[i], arr[rand(i, n - 1)]);

但是我们依旧需要一个长度为 $n$ 的数组。同时现在这个算法是有后效性的,不能像上面一样直接忽略:因为你可能先把 $0$ 号元素换到第 $m$ 位,又把 $m$ 位的 $0$ 号元素换回到第 $1$ 位来。怎么办呢?

我们想到,由于我们只会执行 $m$ 次这样的操作,所以我们不妨用一个 unordered_map 记下除了前 $m$ 个元素的所有元素中,所有发生了 swap 操作元素下标及其值,于是代码就变成了:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include <cstdint>
#include <unordered_map>
#include <vector>

inline std::vector<uint64_t> select(uint64_t n, int m) {
std::unordered_map<uint64_t, uint64_t> rest;
uint64_t tmp[m];
for (uint64_t i = 0; i < m; ++i) tmp[i] = i + 1;
for (uint64_t i = 0; i < m; ++i) {
const uint64_t j = rand<uint64_t>(i, n - 1);
if (j < m) std::swap(tmp[i], tmp[j]);
else if (rest.find(j) == rest.end()) { // 没有找到
rest[j] = tmp[i];
tmp[i] = j + 1;
} else std::swap(tmp[i], rest[j]);
}
return std::vector<uint64_t>(tmp, tmp + m);
}

(这里是为了方便理解,但是为了效率肯定还是要把那个 find 返回的迭代器存下来的,否则会多一次 find 的过程)

这样 unordered_map 中最多只会加入 $m$ 个元素,空间和时间复杂度都是 $O(m)$ 的。

划分!

OI 出题中常常遇到这样的需求:单个数据点范围没有限制,但是所有数据点的总和有限制,这时我们需要把一个大数 $n$ 划分为 $k$ 个有序的部分。其实我们只需要在 $(n+k-1)$ 个点中选出 $(k-1)$ 个点作为分界点,剩下的 $k$ 个区间的长度就是我们选的数了。

由于大数划分的方案和我们这个选点的方案一一对应形成双射,所以可以证得我们的方案生成是均匀随机的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#include <cassert>
#include <type_traits>
#include <vector>

template<class T, class = std::enable_if_t<std::is_integral_v<T>>>
inline std::vector<T> partition(T sum, T count) {
assert(sum >= 0 && count > 0);
const T len = sum + count - 1;
auto ps = choose<T>(0, len - 1, count - 1);
std::sort(ps.begin(), ps.end());
std::vector<T> ret; ret.resize(count);
T last = 0;
for (size_t i = 0; i < ps.size(); ++i) {
ret[i] = ps[i] - last;
last = ps[i] + 1;
}
ret.back() = len - last;
return ret;
}

这样生成的数下界是 $0$,事实上有时候我们还需要要求分出来的每一个数都大于等于一个下界 $l$,实际上也很简单,我们先从 $n$ 个数中划出 $lk$ 的部分预先分配给每个数然后再对剩下的跑刚才的算法即可:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#include <cassert>
#include <type_traits>
#include <vector>

template<class T, class = std::enable_if_t<std::is_integral_v<T>>>
inline std::vector<T> partition(T sum, T count, T min_value = 1) {
if (min_value < 0) min_value = 0;
assert(sum >= 0 && count > 0);
assert(min_value * count <= sum);
const T len = sum + count * (1 - min_value) - 1;
auto ps = choose<T>(0, len - 1, count - 1);
std::sort(ps.begin(), ps.end());
std::vector<T> ret; ret.resize(count);
T last = 0;
for (size_t i = 0; i < ps.size(); ++i) {
ret[i] = ps[i] - last + min_value;
last = ps[i] + 1;
}
ret.back() = len - last + min_value;
return ret;
}

不过带上界的序列均匀随机划分我个人目前还没有找到什么办法… 如果有知道的可以联系我,QQ 洛谷均可,谢谢啦~

括号序列!

我们如下递归定义合法的括号序列:

  • 空串是合法的括号序列
  • 如果 A 是合法的括号序列,那么 (A)A() 也是合法的括号序列

众所周知,长度为 $2n$ 的合法括号序列有 $C_n$(卡特兰数)个。那么我们除了枚举所有 $C_n$ 个括号序列有更好的随机生成合法括号序列的方法吗?

下文在数学公式中用到左右括号的地方都会用下划线表示来与数学中的括号区分。

我们首先来给出几个定义。一个括号序列是平衡的,当且仅当其中的左括号和右括号一样多。为了方便理解,我们先把一个长度为 $2n$ 的平衡括号序列(不一定合法)对应到一条从 $(0,0)$ 走到 $(2n,0)$ 的路径:$\underline($ 左括号对应向左上走 $(1,1)$,$\underline)$ 右括号对应向右下走 $(1,-1)$;并称其为括号路径。例如,下图给出了 $\underline{))()((()()}$ 的括号路径:

示例

然后我们定义一个平衡括号序列 $w$ 的缺陷度 $f(w)$(flaw)为 $k$,当且仅当它所对应的括号路径在 $x$ 轴下的边有 $2k$ 条。很显然每个平衡括号序列的括号路径在 $x$ 轴下的边总有偶数条,因此每个平衡括号序列总有一个缺陷度。例如,上面的路径中有 $6$ 条边在 $x$ 轴下,因此它的缺陷度为 $3$。

我们定义如果 $A$ 是一个括号序列,则 $r(A)$ 为 $A$ 翻转后得到的括号序列。例如 $r(\underline{((()))})=\underline{)))(((}$。我们称一个平衡括号序列 $A$ 不可约,当且仅当不存在两个合法的非空 平衡 括号序列 $B,C$ 满足 $BC=A$($B,C$ 顺次相接与 $A$ 相同)。

引理 1:一个平衡括号序列 $w$ 有且仅有一种分解 $w=w_1w_2w_3\cdots w_k$,其中 $w_i$($1\le i\le k$)为不可约的平衡括号序列。

记 $B_n$ 为所有长度为 $2n$ 的平衡括号序列形成的集合,$B_{n,k}$ 为 $B_n$ 中缺陷度为 $k$ 的平衡括号序列形成的集合。我们现在建立一个从 $B_n$ 到它自身的映射 $\Phi(w)$ 如下(实际上是 Something Comforting 那道题的灵感来源):

我们首先找到一个非空的平衡合法括号序列 $p$ 使得 $p$ 是 $w$ 的前缀且 $p$ 不可约。可以证明这样的 $p$ 一定存在。然后我们把 $w$ 写作 $ps$ 的形式,其中 $s$ 也是一个平衡括号序列。如果 $p$ 不是合法的括号序列,其必定能被写作$\underline)t\underline($ 的形式,其中 $t$ 是可空的平衡括号序列。那么:

$$ \Phi(w)=\begin{cases} p\Phi(s),&p\text{合法}\\ \underline(\Phi(s)\underline)r(t),&p\text{不合法} \end{cases} $$

定理 1:$\Phi(w)$ 建立了从 $B_n$ 到 $B_{n,0}$ 的映射。

也就是说 $\Phi(w)$ 将任意一个平衡括号序列转换为了一个合法的括号序列。考虑归纳证明。首先 $p$ 合法时,$\Phi(w)$ 是否合法等价于 $\Phi(s)$ 是否合法,通过归纳得之;当 $p$ 不合法时,由于 $p$ 不可约,我们得知 $t$ 的缺陷度为 $\frac12|t|$,也就是说 $r(t)$ 的缺陷度为 $0$,它是一个合法的括号序列;又由于 $\underline(\Phi(s)\underline)$ 通过归纳可知是合法的,所以当 $p$ 不合法时 $\Phi(w)$ 也是一个合法的括号序列。我们得证。

定理 2:$\Phi(w)$ 对于任意 $0\le k\le n$ 建立了 $B_{n,k}$ 到 $B_{n,0}$ 的单射。

假设存在两个不相等的平衡括号序列 $w_1$ 和 $w_2$,且 $f(w_1)=f(w_2)=k$,于是我们要证明 $\Phi(w_1)\ne\Phi(w_2)$。沿袭上面的拆分将其写作 $w_1=p_1s_1$,$w_2=p_2s_2$,若 $s_1$ 不合法则 $s_1=\underline)t_1\underline($,$s_2$ 不合法则 $s_2=\underline)t_2\underline($。分类讨论(归纳证明):

  1. $p_1$ 和 $p_2$ 都合法:根据引理 1 我们得知 $\Phi(w_1)\ne\Phi(w_2)$ 等价于 $\Phi(s_1)\ne\Phi(s_2)$,通过归纳得之。

  2. $p_1$ 和 $p_2$ 都不合法:根据引理 1 我们得知 $\Phi(w_1)\ne\Phi(w_2)$ 等价与 $\Phi(t_1)\ne\Phi(t_2)$,通过归纳得之。

  3. $p_1$ 合法而 $p_2$ 不合法:我们要证 $p_1\Phi(s_1)\ne\underline(\Phi(s_2)\underline)r(t_2)$。由于 $p_1$ 和 $\underline(\Phi(s_2)\underline)$ 都不可约,根据引理 1 我们得知 $p_1=\underline(\Phi(s_2)\underline)$ 且 $\Phi(s_1)=r(t_2)$。

    首先 $f(\Phi(s_1))=f(p_1)+f(\Phi(s_1))=f(w_1)=k$,然后由于对于任何平衡括号序列都有 $2f(w)\le |w|$,因此 $k\le \frac12|s_1|$,又有 $|s_1|=|\Phi(s_1)|=|r(s_2)|=|p_2|-2$,即 $k\le \frac12(|p_2|-2)<\frac12|p_2|$。由于 $p_2$ 不合法且不可约,那么 $f(p_2)=\frac12|p_2|$。因为 $f(w_2)=f(p_2)+f(s_2)=k$ 且 $f(s_2)\ge0$,则 $f(p_2)\le k$。那么我们有 $k<f(p_2)\le k$,出现矛盾,假设不成立。

  4. $p_1$ 不合法而 $p_2$ 合法:与第三种情况相同。

综上,我们得证。

定理 3:$\Phi(w)$ 对于任意 $0\le k\le n$ 建立了 $B_{n,k}$ 到 $B_{n,0}$ 的双射。

结合定理 2 和单射的性质我们知道 $|B_{n,k}|\le|B_{n,0}|$,再考虑到根据定义有
$$
\sum_{k=0}^n|B_{n,k}|=|B_n|=\binom{2n}n
$$
由于 $|B_{n,0}|=C_n$,如果存在 $k$ 使得 $|B_{n,k}|<|B_{n,0}|$,那么我们有
$$
\binom{2n}n\le\sum_{k=0}^n|B_{n,k}|<(n+1)C_n=\binom{2n}n
$$
出现矛盾,假设不成立。因此对于任何 $0\le k\le n$ 都有 $|B_{n,k}|=|B_{n,0}|$。结合定理 2,我们得证。

接下来就很愉快啦~我们只需要随机生成一个平衡括号序列(即在 $2n$ 个位置中随机选 $n$ 个设为左括号,其它设为右括号)并对其使用 $\Phi$ 函数即可均匀随机产生一个合法的括号序列。显然算法是线性的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include <algorithm>
#include <string>
#include <map>

inline std::string brackets(size_t n) {
const size_t len = n << 1;
bool arr[len];
std::fill(arr, arr + n, 1);
std::fill(arr + n, arr + len, 0);
std::shuffle(arr, arr + len, engine);
size_t start = 0, end = len;
while (true) {
size_t lef_count = 0, rig_count = 0;
bool ok = true;
for (size_t i = start, j; i < end; ++i) {
++(arr[i]? rig_count: lef_count);
if (lef_count >= rig_count) continue;
for (j = i + 1; j < end; ++j) {
++(arr[j]? rig_count: lef_count);
if (rig_count > lef_count) continue;
// ( ) ) ) ) ( ( ( ) ) ) ( ( (
// i ---S--- j -----T-----

// we swap S and T and then flip S, i and j
const size_t len = j - i - 1;
std::rotate(arr + i + 1, arr + j + 1, arr + end);
std::copy_backward(arr + end - len - 1, arr + end - 1, arr + end);
std::transform(arr + end - len, arr + end, arr + end - len, std::logical_not());
arr[i] = 0; arr[end - len - 1] = 1;
start = i + 1; end = end - len - 1;
ok = false;
break;
}
}
if (ok) break;
}
std::string ret; ret.resize(len);
for (size_t i = 0; i < len; ++i) ret[i] = "()"[arr[i]];
return ret;
}

int main() {
// 验证均匀随机性
std::map<std::string, int> rec;
const int n = 4, T = 100000;
for (int i = 0; i < T; ++i)
++rec[brackets(n)];
for (auto &[str, cnt] : rec)
std::cout
<< str << ": "
<< ((double)cnt / T) << std::endl;
}

下面是在我电脑上运行一次的输出:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
(((()))): 0.07066
((()())): 0.07385
((())()): 0.07146
((()))(): 0.0732
(()(())): 0.07022
(()()()): 0.06976
(()())(): 0.07197
(())(()): 0.07146
(())()(): 0.0715
()((())): 0.07129
()(()()): 0.07212
()(())(): 0.07061
()()(()): 0.07121
()()()(): 0.07069

可以看到还是分布得很均匀的。

二叉树!

注意这里的二叉树是有根的。

定理 4:所有 $n$ 个结点的有根二叉树形成的集合与 $B_{n,0}$ 构成双射。

首先我们直到 $n$ 个结点的有根二叉树共有 $C_n$ 种,因此两个集合大小是相同的。然后我们考虑构造一个映射 $g(tr)$ 来将一个 $n$ 个结点的二叉树映射到一个长度为 $2n$ 的括号序列:
$$
g(tr)=\underline(g(tr.lson)\underline)g(tr.rson)
$$
$|g(tr)|=2n$ 可以显然用归纳法得出。考虑证明这是一个双射。还是归纳证明。根据引理 1 我们得知给定 $g(tr)$ 我们可以唯一确定 $g(tr.lson)$ 的值(因为 $\underline(g(tr.lson)\underline)$ 不可约),那么 $g(tr.rson)$ 的值也随之确定,又由于 根据 $g(tr.lson)$ 和 $g(tr.rson)$ 可以还原出原子树 已经被归纳证明,那么显然 $g(tr)$ 是一个满射。又因为两个集合大小相同,从而我们得证。

那么显然我们只需要实现这个线性的变换即可:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include <cassert>
#include <iostream>
#include <memory>
#include <string>

template<class T>
using Ptr = typename std::unique_ptr<T>;

struct node {
Ptr<node> lson, rson;
node(Ptr<node> lson, Ptr<node> rson):
lson(std::move(lson)), rson(std::move(rson)) {}
};

// 二叉树 -> 括号序列
inline std::string brackets(const Ptr<node> &tr) {
if (!tr) return "";
std::string ret;
ret += '(';
ret += brackets(tr->lson);
ret += ')';
ret += brackets(tr->rson);
return ret;
}

// ... 省略上面的 brackets(size_t) 部分 ...
// 括号序列 -> 二叉树
inline Ptr<node> binary_tree(const std::string &str) {
if (str.empty()) return nullptr;
int i = 0, sum = 0;
do sum += str[i++] == '('? 1: -1; while (sum);
// 此处 [0, i) 是 str 的不可约前缀
// [1, i - 1) 为 g(tr.lson) 的对应部分,[i, n) 是 g(tr.rson) 的对应部分
auto lstr = str.substr(1, i - 2), rstr = str.substr(i); // C++ 的 substr 格式是 (first, count)
return std::make_unique<node>(binary_tree(lstr), binary_tree(rstr));
}
inline Ptr<node> binary_tree(size_t n) { return binary_tree(brackets(n)); }

int main() {
const std::string seq = "()(())()()(()())";
assert(brackets(binary_tree(seq)) == seq);
std::cout << brackets(binary_tree(5)) << std::endl;
}

下面是在我电脑上运行一次的输出:

1
()()(()())

这个不用多说,由于 Prüfer 告诉我们 $n$ 个结点的 有标号 树($n\ge2$)和长度为 $n-2$ 值域为 $[1,n]$ 的序列(Prüfer 序列)构成双射,因此我们只需要随机生成一个 Prüfer 序列再把它转成树即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#include <algorithm>

inline std::vector<std::pair<int, int>> tree(int n) {
assert(n > 0); std::vector<std::pair<int, int>> ret;
if (n == 1) return ret;
std::vector<int> prufer; prufer.resize(n - 2);
for (int &v : prufer) v = rand(1, n);

int deg[n + 1];
std::fill(deg + 1, deg + n + 1, 1);
for (int x : prufer) ++deg[x];
int ptr = 0; while (deg[++ptr] != 1);
int leaf = ptr;
for (int x : prufer) {
ret.emplace_back(x, leaf);
if (--deg[x] == 1 && x < ptr) leaf = x;
else {
while (deg[++ptr] != 1);
leaf = ptr;
}
}
ret.emplace_back(leaf, n);
return ret;
}

int main() {
auto tr = tree(5);
for (auto [x, y] : tr)
std::cout << x << ' ' << y << std::endl;
}

下面是在我电脑上运行一次的输出:

1
2
3
4
1 2
5 3
1 4
1 5

形成的是下面这样的一棵树:

随机树

但是实际造数据时,大多数出题人都使用的是如下造树的方式:

1
2
for (int i = 2; i <= n; ++i)
std::cout << rand(1, i) << ' ' < i << std::endl;

这样显然满足树的性质,但并不均匀随机。可以证得这样生成的树以 $1$ 为根的期望高度是 $O(\log n)$ 的。关于这个证明,我以前由于太 naive 所以写的是假的,所以现在搬一个 EI 的做法。大概的思路是把对 $dep_i$ 求和改为对 $2^{dep_i}$ 求和,使用了一种优于线性的对其最大值进行观测的方式(个人理解)。如下:

定义我们求的树高 $h=\max{d_i}$。令 $t_i=2^{dep_i}$,可得 $2^h=\max{t_i}\le\sum t_i$($t_i\ge0$)。 又由期望的线性性和琴生不等式可知 $2^{\mathbb E[X]}\le\mathbb E[2^X]$,因此 $2^{\mathbb E[h]}\le\mathbb E[2^h]\le\sum\mathbb E[t_i]$。

令 $f_i=\mathbb E[t_i]$,有

$$ \begin{aligned} f_i&=\frac2{i-1}\sum_{j=1}^{i-1}f_j\\ S_i-S_{i-1}&=\frac2{i-1}S_{i-1}\\ S_i&=\frac{i+1}{i-1}S_{i-1}\\ S_i&=\frac{i(i+1)}2 \end{aligned} $$

由此 $\mathbb E[2^h]\le\sum\mathbb E[t_i]=\frac{n(n+1)}2$,则 $\mathbb E[h]\le\log_2\frac{n(n+1)}2=O(\log n)$,我们得证。% EI

图!

要求生成给定点数 $n$ 和边数 $m$ 的无向图。啊这个我会!从 $\binom{2n}n$ 中随机选 $m$ 条加上就好了。但是要求联通呢?

这个 mivik 只会暴力啦。每次随机按照上面的做法生成一张随机图,如果不是联通的就再生成一次,直到图联通为止。这样可以保证是均匀随机的。时间复杂度更低的做法 mivik 就不会了… 还是和上面一样,如果有知道的可以私我。

大部分出题人使用的是下面的做法:先随机生成一棵 $n$ 个结点的有标号树,然后再往上面随机加边。这样并不是均匀随机的,原因很简单:你每张无向图被生成的方案数都是
$$
\text{生成树个数}\cdot \binom{\binom{2n}n}{m-(n-1)}
$$
后面那个是固定的,但前面那个并不固定。例如下图:

两张无向图

虽然两张无向图边数和点数都为 $4$,但第一张有 $3$ 棵生成树,第二张有 $4$ 棵。也就是说,使用上面的生成算法,生成第二张图的概率会更大一些。

The End?

鉴于 mivik 知识有限就只能写到这里了。当然如果有小伙伴知道其它关于随机的有趣知识可以在评论区里面放送~

本文提到的所有算法都已经在我的项目 micrandom.h 中实现,感兴趣的小伙伴可以去看看。

参考资料
  1. Fisher-Yates Shuffle - Wikipedia
  2. Atkinson, M.D. and J.-R. Sack, Generating binary trees at random, Information Processing Letters 41 (1992) 21-23.
作者

Mivik

发布于

2021-01-19

更新于

2021-05-31

许可协议

评论