从随机源到熵池：组合抽样中的随机性回收

本文含 AI 生成内容。前些天看到《杀戮尖塔 2》出了伪随机数生成相关性的问题，就想到把这个算法让 AI 整理一下发出来。其实代码主要是我年初的时候写的，但是应该没什么实用价值，因为真随机数实际上并不稀缺。

我有一天看到了一个类似这样的抽奖问题：

只有一个 8-bit 随机数生成器，怎样高效地从 30000 个人中抽出 20000 个中奖者？

第一反应通常是做随机抽样：多读几个随机字节，取模，拒绝采样，直到抽够为止。这个思路可以保证结果无偏，但它忽略了一个更底层的成本：随机性本身。

从 30000 人中抽出 20000 人，与抽出 10000 个未中奖者等价。仅输出这个组合，理论上就需要

log2 C(30000, 20000) = log2 C(30000, 10000) ≈ 27541.20 bit

随机源一次只能给 1 个字节，因此信息论下界是

ceil(27541.20 / 8) = 3443 bytes

如果某个实现消耗了两倍甚至三倍的随机字节，它仍然可以是正确的均匀抽样；但从熵效率看，它离最优还很远。

这篇文章讨论的问题是：拒绝采样里没有变成输出的随机状态，是否一定要丢掉？答案是否定的。我们可以维护一个熵池（entropy pool），把已经读入但尚未被输出消耗的随机状态保留下来。进一步地，在组合抽样中，如果中间过程生成了一个有序样本，而最终只需要无序集合，那么顺序里包含的信息也可以回收。

下面从最简单的范围抽样开始，逐步推到最终的组合抽样算法。

1. 把随机性当作资源

先看一个子问题：如何用 8-bit 随机源生成一个 [0, 30000) 上的均匀整数？

直接做法是读两个字节，得到 [0, 65536) 上的整数。如果它小于 30000 就接受，否则重来。接受率为

30000 / 65536 ≈ 45.78%

平均每产生一个目标整数，需要消耗

16 / 0.4578 ≈ 34.95 bit

而一个 [0, 30000) 均匀整数本身只包含

log2 30000 ≈ 14.87 bit

的信息量。

稍微改进一下，可以只读 15 bit，因为 2^15 = 32768 已经覆盖 30000。此时拒绝率降为

(32768 - 30000) / 32768 ≈ 8.45%

这已经好很多，但仍然会丢掉一段状态空间。关键在于，这段状态空间本身也是均匀随机的。它没有成为这一次的输出，不代表它没有价值。

本文采用的视角是：随机采样算法不仅在“产生结果”，也在“搬运信息”。理想情况下，每一 bit 随机性要么进入输出，要么留在可复用状态中，而不是被控制流悄悄消耗掉。

2. 熵池：保存尚未使用的均匀状态

约定

[n] = {0, 1, ..., n - 1}

熵池维护一对整数

(B, M)

不变量是：

B 在 [M] 上均匀分布

可以把 M 理解为当前池子的状态空间大小，B 是其中的状态编号。池子里的账面熵为

H_pool = log2 M

后续操作只要能解释成状态空间之间的双射，就可以精确保持均匀性。

2.1 注入随机字节

每次从随机源读入一个字节 U ∈ [256]，将池状态更新为

B' = 256B + U
M' = 256M

这是 [M] × [256] 到 [256M] 的双射。因此，如果旧的 B 在 [M] 上均匀，U 独立且均匀，那么新的 B' 在 [M'] 上也均匀。

这一步只是在池中追加 8 bit 新熵：

log2 M' - log2 M = 8 bit

没有损失。

2.2 从熵池生成 [n] 上的均匀整数

假设当前状态空间大小为

M = qn + r,  0 <= r < n

其中

q = floor(M / n)
r = M mod n

如果 B < qn，说明当前状态落在一段可以被整齐切成 n 份的区域中。此时定义

X  = B mod n
B' = floor(B / n)
M' = q

映射

B -> (X, B') = (B mod n, floor(B / n))

是 [qn] 到 [n] × [q] 的双射。因此在接受分支内部：

  • X 在 [n] 上均匀；
  • B' 在 [q] 上均匀；
  • X 与 B' 独立。

需要注意的是，这并不等于整个操作无损。接受之前池子的账面熵是 log2 M；接受之后，输出携带 log2 n，池子剩余 log2 q。二者之和为

log2 n + log2 q = log2(qn) = log2(M - r)

接受分支上的账面损失为

loss_accept
= log2 M - log2(qn)
= log2(M / (M - r))
= -log2(1 - r / M)

损失来自“这次落在接受区域”这个分支事件本身。程序知道这个事件发生了，但这个事实既没有成为用户输出，也没有被编码回熵池。

2.3 拒绝分支也可以回收

如果 B >= qn，传统拒绝采样通常会丢掉当前随机状态，然后重新读随机数。熵池可以做得更细。

拒绝区域为

{qn, qn + 1, ..., M - 1}

大小为 r。可以把它重新编号为

B' = B - qn
M' = r

这是拒绝区域到 [r] 的双射。也就是说，被拒绝区域里的 offset 没有被丢弃，而是留在了池中。

拒绝分支上的账面损失是

loss_reject
= log2 M - log2 r
= log2(M / r)
= -log2(r / M)

这里损失的仍然不是 offset，而是“进入拒绝分支”这一事实。

令

epsilon = r / M

一次尝试的接受概率为 1 - epsilon，拒绝概率为 epsilon。分支事件带来的期望损失正好是二元熵：

E[loss]
= (1 - epsilon) · [-log2(1 - epsilon)]
  + epsilon · [-log2 epsilon]
= H2(epsilon)

因此，熵池能回收拒绝区域内部的随机状态，但不能自动消除控制流分支本身携带的信息。

2.4 REJECTION_BOUND 的作用

实际实现不会在任意 M 下立刻尝试抽 [n]。我使用了类似下面的逻辑提前扩池：

1
2
3
while (self.m % n) * REJECTION_BOUND >= self.m {
    self.append_random_byte()?;
}

这里 self.m % n 就是 r。这段代码等价于：只要

epsilon = r / M >= 1 / REJECTION_BOUND

就继续读入随机字节，直到拒绝区域比例足够小。

它的意义不只是减少循环次数，更重要的是控制熵损失。开始抽样之前可以保证

epsilon < 1 / REJECTION_BOUND

因此接受分支上的损失满足

loss_accept
= -log2(1 - epsilon)
< -log2(1 - 1 / REJECTION_BOUND)

当 REJECTION_BOUND 较大时，这个量约为

1 / (REJECTION_BOUND · ln 2)

拒绝分支单次损失可能较大，但发生概率也是 epsilon。从期望看，单次尝试的分支损失满足

E[loss] <= H2(1 / REJECTION_BOUND)

如果考虑一次 gen_range(n) 可能经历多次拒绝，接受概率至少为

1 - 1 / REJECTION_BOUND

所以期望尝试次数至多为

1 / (1 - 1 / REJECTION_BOUND)

于是每次 gen_range(n) 的期望分支损失可以粗略控制在

H2(1 / REJECTION_BOUND)
/ (1 - 1 / REJECTION_BOUND)

这个界不追求紧，只说明一个工程事实：REJECTION_BOUND 控制的是随机性损耗，而不仅是拒绝采样的循环次数。

2.5 池操作复杂度

本文实现中的 B 和 M 使用机器整数，例如 u64。在这个设定下，一次池操作可以近似看作 O(1)：

  • 注入字节：乘法和加法；
  • 尝试抽样：除法、取模和比较；
  • 回收 digit：乘法和加法。

空间复杂度为 O(1)。

如果推广到任意精度整数，复杂度应按 log M 的字长计，除法和取模会更重。本文关注的是工程实现，因此后文默认池操作在机器整数上近似常数时间。

3. 从均匀整数到随机排列

有了 gen_range(n)，可以直接用 Fisher–Yates 洗牌生成无重复有序样本：

1
2
3
4
for i in 0..m {
    let r = i + gen_range(n - i);
    c.swap(i, r);
}

第 i 步从剩余的 n - i 个元素中均匀选一个放到位置 i。因此前 m 个元素构成一个长度为 m 的无重复有序样本。任意一个有序样本出现的概率为

1 / (n · (n - 1) · ... · (n - m + 1)) = 1 / P(n, m)

其中

P(n, m) = n! / (n - m)!

如果需要构造完整候选数组 c = [0, 1, ..., n - 1]，初始化为 O(n) 时间和 O(n) 空间；前 m 步洗牌需要 O(m) 次 gen_range 和 swap。

但抽奖要的是组合，不是排列。中奖名单不关心顺序。如果先生成长度为 m 的随机排列，再把顺序扔掉，会额外浪费

log2(m!)

bit 的顺序信息。

对于 m = 20000，这不是一个小常数。因此，组合抽样除了要让 gen_range 接近最优，还要处理有序样本和无序集合之间的这部分熵差。

4. 方案一：生成排列后回收 Lehmer code

一个直接方案是：

  1. 先生成长度为 m 的无重复有序样本；
  2. 忽略顺序，把它作为组合输出；
  3. 将这 m 个元素在集合内部的相对顺序编码成 Lehmer code；
  4. 把 Lehmer code 表示的 [m!] 状态回收到熵池。

这个方案的正确性很直观：

有序样本  <->  无序集合 + 集合内部排列

固定一个集合后，内部排列共有 m! 种。Lehmer code 正好把这 m! 种排列编码为一个混合进制数。因此，均匀有序样本可以分解成：

均匀的 m 元集合 + 独立均匀的 [m!] 状态

从熵账本看：

log2 P(n, m) - log2 C(n, m) = log2(m!)

Lehmer code 回收的正是这部分。

这个方案适合作为概念验证，但工程形态不够轻：

  • 生成有序样本需要 O(n + m) 时间和 O(n) 空间；
  • 计算集合内部相对排列通常需要排序或秩查询，约为 O(m log m)；
  • 如果把 Lehmer code 合成为一个 [m!] 上的大整数，需要大整数运算；
  • 如果改成逐位回收混合进制 digit，实现就会接近后文的最终方案。

所以，它说明了“顺序熵可以回收”，但不是最顺手的实现方式。

5. 方案二：按最小值分布递归

另一个思路是直接在组合空间上做递归。

从 [n] 中选出 m 个元素。如果最小值为 x，剩下的 m - 1 个元素必须从 {x + 1, ..., n - 1} 中选。因此

Pr[min = x] = C(n - x - 1, m - 1) / C(n, m)

其中

0 <= x <= n - m

这给出一个递归算法：先按上述分布抽出最小值 x，输出 x，再在后缀中递归抽出 m - 1 个元素。

关键是第一步不是均匀采样，而是带权采样。权重为

w_x = C(n - x - 1, m - 1)

总权重为

W = sum_x w_x = C(n, m)

5.1 带权采样中的 offset 回收

一般地，若结果 i 的权重为 w_i，总权重为

W = sum_i w_i

可以先从熵池中抽一个均匀整数

Y ∈ [W]

再用前缀和找到唯一的 i，使得

prefix_i <= Y < prefix_i + w_i

令

A = Y - prefix_i

则 A ∈ [w_i]。输出只需要 i，而不需要区间内部的 offset A。因此，条件于输出 i，A 可以回收到熵池：

B' = B · w_i + A
M' = M · w_i

这里 (B, M) 是 gen_range(W) 返回 Y 后留下的池状态。由于 A 在 [w_i] 上均匀，且与原池状态独立，上式仍然保持池状态均匀。

熵账本对应为

log2 W = H(i) + E[log2 w_i]

其中 E[log2 w_i] 就是区间内部 offset 可回收的平均熵。

把它用于最小值递归，就是设置

w_x = C(n - x - 1, m - 1)
W   = C(n, m)

5.2 工程代价

这个方案在数学上直接，但实现成本较高。

对本文的例子，C(30000, 10000) 大约需要 27541 bit 才能表示，远超机器整数范围。每一轮都需要处理组合数权重，通常会引入大整数比较、加减，甚至除法。

粗略看：

  • 全局 unranking 可以在 O(n) 或 O(k log n) 次组合数比较内完成，但每次比较的字长为 Θ(log C(n, k))；
  • 逐轮按最小值带权采样需要 k 轮，每轮都要在组合数权重中定位区间；
  • 空间上至少要保存若干大整数，单个大整数规模就是 Θ(log C(n, k)) bit。

因此，这条路更像“组合数 unranking + 熵回收”。它理论上干净，但不符合本文想要的机器整数实现：每一步只处理局部状态，不维护巨大的组合数表，同时回收有序样本中组合不需要的顺序熵。

6. 最终方案：生成组合时回收插入秩

从 n 个元素中抽 m 个，与抽 n - m 个补集元素等价。因此实际只需抽

k = min(m, n - m)

如果 m > n / 2，就抽未中奖者，最后取补集。

核心代码如下。实现中使用 wabi_tree 提供的 order-statistic set；它类似 BTreeSet，但额外支持 rank_of 这类秩查询。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
let k = min(m, n - m);
let mut selected = OSBTreeSet::new(); // from wabi_tree
let mut c: Vec<u32> = (0..n).collect();

for i in (n - k..n).rev() {
    let r = gen_range(i + 1);
    let t = c[r];
    c.swap(r, i);

    selected.insert(t);

    let j = n - i;                  // 已选元素数量
    let rank = selected.rank_of(t); // rank in [j]
    recycle(rank, j);
}

前半段是 Fisher–Yates 的逆向写法。每轮从当前候选区间 [0, i] 中均匀选出一个元素 t，再把它交换到末尾的已抽取区域。因此，按产生顺序看，t 构成一个长度为 k 的无重复有序样本。

后半段是关键：

1
2
let rank = selected.rank_of(t);
recycle(rank, j);

第 j 轮时，新元素 t 插入当前已选集合后，它的秩一定属于 [j]。这个秩表示：在最终集合的相对顺序中，第 j 个被抽到的元素插在第几个位置。

组合只关心最终集合，不关心这些元素被抽到的先后顺序。因此，插入秩是中间过程生成、但最终输出不需要的信息，可以回收到熵池。

这和 Lehmer code 是同一类分解。Lehmer code 是在完整排列生成之后一次性编码；这里是在生成过程中逐位回收混合进制 digit。

j 从 1 到 k，每轮回收一个 [j] 上的 digit。总状态数为

1 · 2 · ... · k = k!

回收的正是

log2(k!)

这部分顺序熵。

7. 为什么在线回收不会引入偏差

这个算法最容易产生疑问的地方是：刚生成的随机信息被立刻放回熵池，而后续抽样又会使用这个熵池。这会不会形成反馈，破坏组合的均匀性？

不会。可以用一个归纳不变量说明。

第 j 轮结束后，保持以下性质：

  1. 当前已选集合 S_j 在所有大小为 j 的子集上均匀分布；
  2. 熵池状态 (B, M) 仍满足 B 在 [M] 上均匀；
  3. 条件于已经输出的集合，池状态仍然是均匀可用的。若池大小固定，也可表述为池状态与 S_j 独立。

j = 0 时显然成立。

假设第 j - 1 轮结束后不变量成立。第 j 轮调用 gen_range，从剩余元素中均匀选出新元素 T_j。根据熵池范围抽样的性质，这个输出与调用后留下的池状态独立。

加入新元素后：

S_j = S_{j-1} ∪ {T_j}

由于 S_{j-1} 均匀，且 T_j 从剩余元素中均匀选择，S_j 在所有大小为 j 的子集上均匀。

再看新元素在 S_j 中的秩：

R_j = rank(T_j in S_j)

给定最终集合 S_j，最后插入的是其中哪一个元素？有序无重复抽样对元素完全对称，因此 T_j 在 S_j 的 j 个元素中均匀。于是

R_j 在 [j] 上均匀

并且 R_j 与 S_j 独立。

另一方面，gen_range 输出 T_j 后留下的池状态与 T_j 独立；结合归纳假设，它也与之前的 S_{j-1} 独立。因此该池状态与新的 (S_j, R_j) 独立。

于是可以把 R_j 当作一个独立的 j 进制随机 digit 回收到池中：

B' = B · j + R_j
M' = M · j

只要乘法不溢出，这是 [M] × [j] 到 [Mj] 的双射。回收后，池状态仍然均匀，并且仍可用于后续抽样。

从整体看，算法隐含了下面的双射分解：

长度为 k 的有序无重复样本
<->
最终 k 元集合 + 插入秩序列

插入秩序列共有 k! 种状态。算法只是把这个分解在线化了。

8. 复杂度

最终方案的复杂度主要由三部分组成。

首先，候选数组 c 的初始化需要

O(n) time, O(n) space

其次，逆向 Fisher–Yates 循环执行

k = min(m, n - m)

轮。每轮包含一次 gen_range、一次数组访问、一次交换，以及一次插入和秩查询。

最后，rank_of 由 wabi_tree::OSBTreeSet 提供。对这里需要的操作，可以按 O(log k) 估计。

因此，基于 wabi_tree 的整体复杂度为

O(n + k log k) time
O(n + k) space

其中 O(n) 空间来自 Fisher–Yates 候选数组，O(k) 空间来自 order-statistic set。

如果 population 不需要完整 materialize，可以用稀疏交换表替代完整数组，把 Fisher–Yates 的数组空间从 O(n) 降到 O(k)。此时整体空间可以接近 O(k)，时间仍主要由 k 次秩树操作决定，即 O(k log k)，另加稀疏交换表的哈希访问开销。

从随机性消耗看，目标输出熵为

log2 C(n, k)

Fisher–Yates 中间生成的有序样本熵为

log2 P(n, k)

插入秩回收

log2(k!)

因此理想账本为

log2 P(n, k) - log2(k!) = log2 C(n, k)

正好回到组合抽样的信息论下界。

9. 熵账本：30000 人抽 20000 人

回到最初的问题。抽出 20000 个中奖者等价于抽出 10000 个未中奖者，因此

n = 30000
k = 10000

理论下界为

log2 C(30000, 10000) ≈ 27541.197884604342 bit

一次测试结果如下：

consumed_bytes  = 3449
consumed_bits   = 3449 * 8 = 27592
retained_states = 1963181291110000
retained_bits   = log2(retained_states)
                ≈ 50.802114828773 bit

只看输出，效率是

output_efficiency
= log2 C(30000, 10000) / consumed_bits
≈ 27541.197884604342 / 27592
≈ 99.815881%

这个数字低于 100%，主要是因为抽样结束后，熵池里还保留了一部分可复用随机状态。剩余池子有

1963181291110000

种可能状态，对应

log2(1963181291110000) ≈ 50.802114828773 bit

如果把剩余池熵也计入账本，则得到

retained_efficiency
= (输出组合熵 + 剩余池熵) / 消耗随机熵
= (log2 C(30000, 10000) + log2 retained_states) / 27592
≈ 99.9999999979%

对应的账面熵损失约为

loss
= consumed_bits - log2 C(30000, 10000) - log2 retained_states
≈ 27592 - 27541.197884604342 - 50.802114828773
≈ 5.67e-7 bit

这个损失很小，但不是零。原因仍然是 gen_range(n) 中的接受/拒绝分支：只要当前池大小 M 不能被目标范围 n 整除，就会有一小部分分支事件信息既没有进入输出，也没有编码回熵池。

这也解释了 REJECTION_BOUND 的价值。它通过限制

epsilon = (M mod n) / M

避免在拒绝区域比例较大时贸然抽样，从而把每次 gen_range 的分支熵损失控制到很小。即使组合抽样需要大量范围抽样，累计损失仍然可以接近不可见。

在这次测试中，随机源提供的 27592 bit 基本进入了两本账：

  1. log2 C(30000, 10000) bit 成为返回的组合；
  2. log2(1963181291110000) bit 留在熵池中，供后续抽样继续使用。

剩下约 5.67e-7 bit，才是真正的账面熵损失。

10. 与已有工作的关系

这个问题可以放在 randomness recycling 的背景下理解。

Draper 和 Saad 的 Efficient Online Random Sampling via Randomness Recycling 系统研究了在线随机采样中的随机性回收：给定离散随机变量序列，如何复用随机算法已经消耗、但没有成为输出信息的那部分随机性，从而让平均熵成本逼近 Shannon 熵率。论文讨论了均匀采样、逆变换采样、查表采样、alias sampling 和 DDG tree sampling 等方法。

更早的 Knuth–Yao 工作为非均匀随机数生成建立了随机 bit 模型下的复杂度视角。Lumbroso 的 Fast Dice Roller 给出了简洁且近乎最优的离散均匀整数生成算法。Fast Loaded Dice Roller 则将精确离散分布采样推进到更实用的形式。

这些工作共同指向同一个原则：随机采样算法不仅要无偏，也应尽量接近信息论下界。

本文关注的是随机组合抽样中的一个具体构造。组合问题里，一个自然的高熵中间对象是“有序无重复样本”，而最终输出只需要“无序集合”。二者相差

log2(k!)

的顺序熵。

本文的做法是把这部分顺序熵拆成每一轮插入时的秩：

R_j = rank(T_j in S_j),  R_j ∈ [j]

并在生成组合的过程中逐步回收这些秩。

因此，这个算法可以看作 randomness recycling 思想在随机组合抽样上的一个工程化构造：它不是先生成完整排列再事后回收 Lehmer code，而是在线维护“当前集合与熵池独立”的不变量，一边生成集合，一边回收刚暴露出来、但最终组合不需要的顺序 digit。

11. 小结

从 8-bit 随机源出发，这个问题最后落到三个关键词上：双射、条件均匀性和熵账本。

算法可以分成三层：

  1. 熵池层：用 (B, M) 保存尚未被输出消耗的均匀随机状态；
  2. 均匀整数层：用状态空间拆分实现 gen_range(n)，并用 REJECTION_BOUND 控制分支熵损失；
  3. 组合抽样层：用 Fisher–Yates 生成无重复有序样本，同时把每轮新元素的插入秩回收到熵池。

它不依赖巨大的组合数表，也不需要把组合问题硬转成完整排列问题。每一步只做一次普通的均匀选择，然后把组合不需要的顺序信息收回池中。

如果把随机性看作一种资源，这个算法的目标很朴素：每一 bit 要么成为输出，要么继续留在池里，尽量不要凭空消失。

参考资料