新的 arc4random_uniform 实现

本月初， Robert Clausecker 替换了 FreeBSD 的 arc4random_uniform(3)。

arc4random_uniform(3) 是 arc4random(3) 之上封装的一个生成一个较小范围伪随机数的函数。 arc4random(3) 采用密码学安全的伪随机数生成一个在 32-bit 范围，即 $[0,\ 2^{32}-1]$ 内均匀分布的伪随机整数， 此处的随机分布是依靠对称加密算法（目前采用的是 Chacha20）中用于实现加密的伪随机置换（Pseudorandom Permutation）来保证的。

伪随机数满足均匀分布是许多应用场景中需要的。举例来说，如果模拟掷骰子，我们一定是希望骰子六个面出现的概率基本相同。 注意到 $2 ^ {32}$ 并不能被 6 整除，直接将 arc4random(3) 取模，或：

arc4random() % 6

得到的并不是我们所希望的 [0, 5] 范围内的均匀分布的随机数，因为 arc4random() 有可能返回 [4294967292, 4294967295] 这几个数，从而使得结果中 0, 1, 2, 3 这四个数出现的概率略微比 4 和 5 高一些。

为了使 0, 1, 2, 3 这几个数字不要出现的概率比其他数字多，我们可以在样本中直接剔除掉开始或结束的这一段产生 0, 1, 2, 3 这几个在样本空间中多出现过一次的数字的原始输出。这可以是 [0, 3]，也可以是 [4294967292, 4294967295]，具体做法是如果 arc4random() 返回了这一区间的数字，则重新生成一个新的伪随机数。

原算法

FreeBSD 此前采用的是 OpenBSD 的算法：

 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

uint32_t
arc4random_uniform(uint32_t upper_bound)
{
	uint32_t r, min;

	if (upper_bound < 2)
		return 0;

	/* 2**32 % x == (2**32 - x) % x */
	min = -upper_bound % upper_bound;

	/*
	 * This could theoretically loop forever but each retry has
	 * p > 0.5 (worst case, usually far better) of selecting a
	 * number inside the range we need, so it should rarely need
	 * to re-roll.
	 */
	for (;;) {
		r = arc4random();
		if (r >= min)
			break;
	}

	return r % upper_bound;
}

即，排除靠近 0 的这一份总数为 min 的导致采样偏差的样本。

上述代码中 min = -upper_bound % upper_bound 利用了 $k$ 为整数、$m$ 为正整数时任意数 $x$ 与 $x \pm km$ 同余，即 $x \equiv x \pm km \pmod{m}$ 的性质，将 $2^{32} \mod upper\_bound$ 改为更容易表示的 $(2^{32}-upper\_bound) \mod upper\_bound$，即 -upper_bound % upper_bound。

但这样做有一个问题：无论如何我们都要进行至少两次除法运算：第一次是算出希望丢弃的最低阈值 min，第二次是返回结果。 更严重的问题是，这两次的除法中的除数均为 upper_bound，这个数未必是二的整数次方幂 ( $2^n$ )。 对于计算机硬件来说，计算除数不是二的整数次方幂的除法比较复杂，而乘法或除数为二的整数次方幂的计算则相对要简单许多， 因此 CPU 在处理这些计算时更为高效，我们希望能够避免发生这类任意除数除法的计算。

新算法

Daniel Lemire 教授 在 2019 年发表了一个新的算法。

考虑我们有一个可以生成 $[0, 2^L-1]$ 范围内均匀分布的伪随机数，希望基于此生成 $[0, s-1]$ 范围内的均匀分布的伪随机数，其中 $0 \leq s \leq 2^L$ 。此处的 $s$ 就是前面算法中的 $upper\_bound$。

新算法具体来说是这样的：

首先，使用 $[0, 2^L-1]$ 范围内均匀分布的伪随机数来生成一个随机数 $x$ 。

为了便于理解，不妨将这个数 $x$ 除以 $2^L$，这样我们就得到一个在 $[0, 1)$ 范围内、间隔为 $1 \over 2^{L}$ 的均匀分布的定点小数。具体而言，由于 $0 \leq x \leq 2^L-1$，因此 ${0 \over 2^L} \leq {x \over 2^L} \leq {{2^L-1} \over {2^L}}$。使用一个包含 $2L$ 个比特的的整数，我们就可以将这个定点数表达为其高 $L$ 位作为其整数部分（全0），而其低 $L$ 位作为其小数部分（$x$）。

定点小数的乘法规则与整数乘法无异。因此，在将长度为 $L$ 的整数 $x$ 扩展到 $2L$ 个比特之后，将 $x$ 乘以 $s$ 我们便可以得到一个整数部分取值范围为 $[0, s-1]$ 的数。从函数角度，这样一来，我们便把所有满足 $0 \leq x \times s < 2^L$ 的 $x$ 对应的输出都映射到 0、所有满足 $2^L \leq x \times s < 2 \times 2^L$ 的 $x$ 对应的输出都映射到 1，或更一般地，将所有满足 $i \times 2^L \leq x \times s < (i+1) \times 2^L$ 的 $x$ 值所对应的输出都映射到 $i$。

与原算法类似，这样一来依然会出现一些采样偏差，我们需要消除这些偏差。

下文中，将小于或等于 $x$ 的最大整数记作 $\lfloor x \rfloor$。 将大于或等于 $x$ 的最小整数记作 $\lceil x \rceil$。 将 $x$ 除以 $y$ 的整数部分即 $\lfloor x / y \rfloor$ 记作 $x \div y$。

观察不难发现，对于整数 $a, b, s$，对于满足 $b > a > 0$ 并且 $s > 0$ 的情形，若 $b-a$ 能被 $s$ 整除，则在 $[a, b)$ 范围内必然存在 $(b-a) \div s$ 个 $s$ 的倍数。

令 $a = i \times 2^L + (2^L \mod s)$， $b = (i+1) \times 2^L$，因此 $b-a = 2^L - (2^L \mod s)$， 故 $b-a$ 一定能被 $s$ 整除，因而在 $[a, b)$ 即 $[i \times 2^L + (2^L \mod s), (i+1) \times 2^L)$ 范围内，存在且仅存在 $2^L \div s$ 个 $s$ 的倍数。

由于我们的目标是在每一个 $i \times 2^L \leq x \times s < (i+1) \times 2^L$ 范围内可能的 $x$ 取值映射到 输出 $i$ 上，这个结论告诉我们，只要我们排除了所有位于 $[ i \times 2^L, i \times 2^L + (2^L \mod s))$ 范围内的 $x \times s$ 值，便可以获得一致的样本了。注意到此处大量出现的 $2^L$，我们恰好可以把这一判断写作检测该乘积的小数部分上，因为我们可以通过判断 $x \times s \mod 2^L < 2^L \mod s$ 来得到同样的结论。

不过，计算 $2^L \mod s$ 依然要用到除法。由于 $s > 2^L \mod s$，我们可以用 $s$ 代替 $2^L \mod s$ 先做一轮判断，即，只在 $x \times s \mod 2^L < s$ 时才去继续进一步测试是否 $x \times s \mod 2^L < 2^L \mod s$。

目前 FreeBSD 的 arc4random_uniform 实现如下：

 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

/*-
 * SPDX-License-Identifier: 0BSD
 *
 * Copyright (c) Robert Clausecker <fuz@FreeBSD.org>
 * Based on a publication by Daniel Lemire.
 * Public domain where applicable.
 *
 * Daniel Lemire, "Fast Random Integer Generation in an Interval",
 * Association for Computing Machinery, ACM Trans. Model. Comput. Simul.,
 * no. 1, vol. 29, pp. 1--12, New York, NY, USA, January 2019.
 */

#include <stdint.h>
#include <stdlib.h>

uint32_t
arc4random_uniform(uint32_t upper_bound)
{
        uint64_t product;

        /*
         * The paper uses these variable names:
         *
         * L -- log2(UINT32_MAX+1)
         * s -- upper_bound
         * x -- arc4random() return value
         * m -- product
         * l -- (uint32_t)product
         * t -- threshold
         */

        if (upper_bound <= 1)
                return (0);

        product = upper_bound * (uint64_t)arc4random();

        if ((uint32_t)product < upper_bound) {
                uint32_t threshold;

                /* threshold = (2**32 - upper_bound) % upper_bound */
                threshold = -upper_bound % upper_bound;
                while ((uint32_t)product < threshold)
                        product = upper_bound * (uint64_t)arc4random();
        }

        return (product >> 32);
}

以下增加一些注解来方便理解。

第32-33行：对于 upper_bound <= 1 的情形，本函数定义为返回 [0, upper_bound) 的随机数，无论 upper_bound 为 0 或 1，它只能返回 0，因而直接返回 0。

第35-46行为原论文中的算法实现。

我们希望尽量避免除法，因此首先在第35行生成一次随机数并将其乘以 upper_bound， 之后做一次粗略判断（第37行），若乘积的低 32-bit 部分超过了 upper_bound，则无需计算 threshold （余数严格小于除数即 upper_bound）就已经知道它符合无偏移采样的条件，因此可以直接到第46行返回结果。

对于乘积低 32-bit 部分小于 upper_bound 的情况，此时乘积的低 32-bit 部分有机会小于 threshold， 如果运气不好，后续生成的随机数可能继续出现乘积低 32-bit 部分有机会小于 threshold 的情况， 因此我们一次性地在第41行把这个数计算出来。

第40-41行：由于 $2^L \equiv (2^L - upper\_bound) \mod upper\_bound$，我们可以把 threshold 写作 -upper_bound % upper_bound。此处不直接计算 $2^{32} \mod n$ 可以避免使用64-bit除法。

第42-43行：接下来，对于乘积低32-bit低于 upper_bound 的情况，再生成一次随机数并计算乘积，直到得到一个符合条件的乘积为止。

最终在第46行：返回乘积的高32-bit。

最后，感谢 Daniel Lemire、Frank 和 Tony Finch 关于新算法中代码快速路径部分的讨论和帮助。

• ← 【公告】关于 woodpecker.org.cn 网站已经停止维护的说明
• 数据的长期保存 →