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

本文含 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；但这些 digit 不能反过来参与尚未完成的 Fisher–Yates 过程。

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

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. 一次错误的在线化，以及 3 of 6 反例

我最初尝试把“生成排列后回收 Lehmer code”直接在线化：用部分 Fisher–Yates 维护候选数组，每抽出一个新元素，就立刻把它在当前已选集合中的插入秩回收到熵池。

这个方案有一个看似充分的整体计数：

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

右侧的插入秩序列共有 k! 种，所以如果先完整生成有序样本，再统一回收这些秩，证明没有问题。错误发生在“立刻”二字：Fisher–Yates 的候选数组也是随机状态的一部分，回收的秩与数组中剩余元素的排列相关。下一轮再次从池中抽样时，这种相关性会反馈到输出。

combination(3, 6) 已经足以给出精确反例。令初始熵池依次取

(B, M) = (0, 120), (1, 120), ..., (119, 120)

这 120 个状态等概率。均匀算法应当把它们分到 C(6, 3) = 20 个组合上，每个组合恰好出现 6 次。旧实现的计数却是：

出现次数 组合
8 {0,1,4}、{0,3,4}、{1,2,5}、{2,3,4}
4 {0,1,5}、{0,3,5}、{1,2,4}、{2,3,5}
6 其余 12 个组合

更直接地说，整体双射只能证明“全部生成完再回收”正确，不能自动证明逐轮回收正确。在线算法必须证明每一轮刚回收的 digit 都与保留下来、将影响未来输出的状态独立。

6.1 改用 Floyd 集合转移

修正方法是不再维护随机排列的候选数组，而是使用 Floyd 的无放回集合采样过程。仍令

k = min(m, n - m)

并用支持 insert 和 rank_of 的有序统计集合 S 保存已选元素：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
let start = n - k;
let mut selected = OSBTreeSet::new();

for j in start..n {
    let t = gen_range(j + 1);

    if !selected.insert(t) {
        selected.insert(j);
    }

    let rank = selected.rank_of(&t).unwrap();
    let radix = j - start + 1;
    recycle(rank, radix);
}

第 r 轮的 j = n - k + r - 1。轮开始时，S 是 [j] 中的 r - 1 元集合；先均匀生成 t ∈ [j + 1]：

  • 若 t ∉ S，就插入 t；
  • 若 t ∈ S，就改为插入当前的新最大值 j。

两种情况都会得到 [j + 1] 中的 r 元集合。碰撞时新插入的是 j，但回收的必须仍是原始抽样值 t 在新集合中的秩，不能改为回收 j 的秩。这一点正是下一节局部双射成立的关键。

如果 m > n / 2，算法采样的是大小为 n - m 的补集，最后按自然顺序构造其补集。这个确定性变换不会引入新的随机状态或隐藏排列。

7. 严格证明：局部双射保证在线回收无偏

记

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

并对轮数 r = 1, 2, ..., k 做归纳。第 r 轮开始时令

j = n - k + r - 1

归纳不变量是：

  1. S 在 [j] 的所有 r - 1 元子集上均匀；
  2. 给定 S，熵池状态仍然均匀，且与 S 独立。

r = 1 时 S = ∅，结论显然成立。调用 gen_range(j + 1) 后得到均匀的 t ∈ [j + 1]；根据范围抽样的契约，调用后留下的池状态与 (S, t) 独立。接下来只需研究 Floyd 转移

(S, t) -> (S', d)

其中 S' 是本轮结束后的集合，d 是 t 在 S' 中从 0 开始的秩，因此 d ∈ [r]。

这个映射有显式逆映射。给定任意 [j + 1] 中的 r 元集合 S' 和任意 d ∈ [r]：

  1. 将 S' 升序排列，令 t 为从 0 开始编号的第 d 个元素；
  2. 若 j ∈ S'，令 S = S' \ {j}；
  3. 若 j ∉ S'，令 S = S' \ {t}。

验证分三种情形：

  • j ∉ S' 时，t 原先不在 S 中，正向转移插入 t；
  • j ∈ S' 且 t ≠ j 时，t 原先就在 S 中，发生碰撞，正向转移插入 j；
  • t = j 时，j 原先不在 S 中，正向转移直接插入 j。

所以正向与逆向过程逐点互逆，(S, t) <-> (S', d) 是双射。两侧状态数也相等：

C(j, r - 1) · (j + 1)
= C(j + 1, r) · r

(S, t) 在左侧笛卡尔积上均匀，因此 (S', d) 在右侧笛卡尔积上均匀。于是：

  • S' 在所有 r 元子集上均匀；
  • d 在 [r] 上均匀；
  • d 与 S' 独立。

设 gen_range 后留下的池状态为 B ∈ [M]。由于它与 (S', d) 独立，可以执行

B' = B · r + d
M' = M · r

这是 [M] × [r] 到 [Mr] 的双射，所以新池仍然均匀，并且与 S' 独立。归纳不变量由此保持到下一轮。

这次回收也不会导致 u64 乘法溢出。设本轮 gen_range(j + 1) 最后一次成功拆分前的工作模数为 M_work；拆分后池模数为

floor(M_work / (j + 1))

而 r ≤ j + 1，所以回收后的模数不超过 M_work ≤ u64::MAX。

把全部 k 轮的状态数相乘：

∏(j + 1) = n · (n - 1) · ... · (n - k + 1) = P(n, k)
∏r = 1 · 2 · ... · k = k!
P(n, k) = C(n, k) · k!

因此输出集合恰好均匀，算法完整回收了组合不需要的 log2(k!) 顺序信息。与旧方案不同，这里每一轮都有独立成立的局部双射；未来会用到的状态只有集合 S' 和熵池，不存在一个与已回收秩相关的候选排列。

对前面的 combination(3, 6) 做同样的 120 状态穷举，修正后的 20 个组合都恰好出现 6 次；而且对每个固定组合，剩余池值都恰好遍历 0, 1, ..., 5。后一个性质比“输出频数相等”更强，它直接验证了输出集合与剩余池状态独立。

8. 复杂度

Floyd 循环执行

k = min(m, n - m)

轮。每轮调用一次 gen_range，并对 wabi_tree::OSBTreeSet 做常数次插入和秩查询。按这些树操作为 O(log k) 计，随机采样阶段需要

O(k log k) time
O(k) auxiliary space

这里不再有大小为 n 的 Fisher–Yates 候选数组。返回值本身需要保存 m 个元素；走补集分支时，还要顺序扫描 [n] 来构造输出。这是输出物化成本，不是采样算法额外保留的随机状态。

从随机性账本看，第 r 轮抽样的进制依次为

n - k + 1, n - k + 2, ..., n

它们的乘积是 P(n, k)。每轮回收的进制依次为

1, 2, ..., k

乘积是 k!。因此组合层的理想净熵成本为

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

正好达到组合抽样的信息论下界。实际实现中若有极小损失，来源是有限宽熵池在 gen_range 中的接受/拒绝分支，而不是 Floyd 转移或秩回收。

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 等方法。本文的 (B, M) 熵池与这种状态空间回收视角属于同一思想谱系。

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

组合采样这一侧，Floyd 集合转移本身是经典算法。Daniel Ting 的综述 Simple, Optimal Algorithms for Random Sampling Without Replacement 对多种无放回采样方法及其适用场景做了梳理。本文并不把 Floyd 采样本身作为新算法。

本文关注的是两条思路的交点：在 Floyd 第 r 轮中，把转移写成局部双射

(S, t) <-> (S', d)

其中

d = rank of t in S',  d ∈ [r]

再把 d 立即回收到熵池。局部双射既给出组合均匀性的证明，也给出在线回收后集合与池状态独立的证明。

这些工作共同指向同一个原则：随机采样算法不仅要无偏，也应尽量接近信息论下界；而只要随机信息会被在线复用，证明对象就不能只有最终输出，还必须包含所有会影响未来的内部状态。

11. 小结

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

算法可以分成三层：

  1. 熵池层：用 (B, M) 保存尚未被输出消耗的均匀随机状态；
  2. 均匀整数层：用状态空间拆分实现 gen_range(n)，并用 REJECTION_BOUND 控制分支熵损失；
  3. 组合抽样层：用 Floyd 转移直接维护集合，并把原始抽样值在新集合中的秩回收到熵池。

这个过程也留下了一个重要教训：全局双射不等于可以任意在线回收。旧方案对完整有序样本的计数没有错，错的是忽略了 Fisher–Yates 剩余数组中的隐藏排列。修正后的算法不依赖巨大的组合数表，也不需要大小为 n 的候选数组；每一轮都用显式可逆的 Floyd 转移同时证明无偏性和独立性。

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

参考资料