使用 rand5() 构造 rand7() 之最佳实现

问题

给定一个函数 rand5()，它可以随机、等概率地生成 0 到 4 之间的整数。现在希望使用该函数和确定性的计算过程，实现函数 rand7()，它可以随机、等概率地生成 0 到 6 之间的整数。

解法

简化的类似问题

这个问题如果反过来会容易解决得多。如果我们有一个函数 rand7()，要如何实现函数 rand5() 呢？只需要不断调用函数 rand7()，当输出结果小于 5 时，直接返回，否则再次调用，直到生成要求范围内的数字。

fun rand5(): Int {
    while(true) {
        val n = rand7()
        if (n < 5) {
            return n
        }
    }
}

直觉告诉我们，这个方法没有对 0 到 4 以内的任何一个数字有所偏袒，其任意两次输出之间也不存在任何关联性，因此满足随机、等概率的要求。

前 $n - 1$ 次生成 5 以外的数字，而第 $n$ 次生成 0 到 4 以内的特定数字的概率为 $\frac{1}{7} \cdot \left( \frac{2}{7} \right) ^{n - 1}$。
因此，生成 0 到 4 以内的特定数字的总概率为

$$\begin{aligned} P(0..5) &= \frac{1}{7} + \frac{1}{7} \cdot \frac{2}{7} + \frac{1}{7} \cdot \frac{2}{7} \cdot \frac{2}{7} + \cdots \\ &= \frac{1}{7} \cdot \sum_{n=0}^{\infty} \left( \frac{2}{7} \right) ^n \\ &= \frac{1}{5} \end{aligned}$$

确实符合要求。

扩展状态空间

但这里的问题是利用 rand5() 实现 rand7()。如果沿用上面的思路，我们需要设法扩展 rand5() 的值域，同时保证结果随机、等概率。

rand5() 的输出结果总共只有 5 种可能，而 rand7() 有 7 种可能的输出。任何确定性的算法，都不可能使输出的可能状态数大于输出的状态数，若非如此，则必然存在某一个输入状态对应多于一个输出结果，违背了确定性算法“确定的输入状态对应唯一确定的输出状态”的定义。这意味着我们需要调用 rand5() 多于一次，才能满足题目的要求。

假设使用两个 rand5() 的输出结果，其输出的状态数有 $5^2 = 25$ 种，大于我们需要的 7 种。接下来的问题，就是如何将这 25 种输出状态再次作为输入，随机、等概率地映射到目标的 7 种状态。

首先，为这 25 种输出状态安排唯一的编号。容易想到的编号方式如下：

fun numbered(a: Int, b: Int) = a * 5 + b

这个编号正好是 0 到 24 的正整数，每个整数的概率都相等，我们可以把这个函数叫做 rand25()：

fun rand25() = rand5() * 5 + rand5()

这样就回到容易解决的问题上来了，可以写出：

fun rand5(): Int {
    while(true) {
        val n = rand25()
        if (n < 7) {
            return n
        }
    }
}

优化效率

这么看问题是解决了，但似乎效率有点低。我们只利用了 25 种输出状态的 7 种，剩余 18 种结果都被丢弃了。
如何衡量这个算法的效率呢？可以计算，每次调用 rand7() 生成一个结果，平均需要调用多少次 rand5()。调用次数越多，说明算法的效率越低:

$$\begin{aligned} E(\mathrm{number}\ \mathrm{of}\ rand5()\ \mathrm{calls}) &= \sum_{n=1}^{\infty} 2n \cdot \frac{7}{25} \left( \frac{18}{25} \right) ^{n - 1} \\ &= \frac{50}{7} \\ &\approx 7.14 \end{aligned}$$

效率确实不太高！

怎么能提高效率呢？可以尝试尽量利用而非丢弃生成的随机数。例如，下面的方法，能有效利用 25 种输出状态的 21 种，于是平均调用次数减少到 $\frac{50}{21} \approx 2.38$ 次。

fun rand7(): Int {
    while(true) {
        val num = rand25()
        if (num < 21) {
            return num % 7
        }
    }
}

这个算法还有提升的空间吗？
可以观察到，21 到 24 之间的输入，可以映射为包含 4 个状态的随机数生成器；
此外，对小于 21 的输入除以 7 的商，可以映射为包含 3 个状态的随机数生成器。

这些“随机性”都被浪费掉了，可以尝试将其重新引入到随机数生成中，以减少 rand5() 的调用次数。

以下是重新利用随机性并提升生成效率的版本。经测试，该算法的平均调用次数减少到了 1.63 次。很大的进步！

fun Int.power(x: Int): Int {
    var n = 1
    repeat(x) { n *= this }
    return n
}

fun rand5(): Int = Random.nextInt(5)

val pool = arrayListOf<Int>()
fun rand7() = rand7(pool, rand5)
/**
 * 将 m 位 5 进制数转换为 n 位 7 进制数，需要丢弃的比例为 (5^m - 7^n) / 5^m。
 * 运算限制在 32 位整数范围内，使损失最小的组合为 (m, n) = (11, 9)
 */
val exp5 = 11
val exp7 = 9
val exp7Rem = 8
val bMax = 7.power(exp7)
fun rand7(pool: MutableList<Int>,
          rand5: () -> Int): Int {
    
    // 调用 rand5() 11 次，生成 0 到 5^11 - 1 范围内的随机数
    fun rand5Power11(): Int {
        var r = rand5()
        repeat(exp5 - 1) { r = 5 * r + rand5() }
        return r
    }
    
    // 将给定的数字 r 转换为 nDigits 位的 7 进制数（高位用 0 填充），
    // 并将每一位添加到随机数池 pool
    fun supplyPoolWith(r: Int, nDigits: Int) {
        var rVar = r
        repeat(nDigits) {
            pool.add(rVar % 7)
            rVar /= 7
        }
    }
    
    // 补充随机数池
    fun supplyPool() {
        val r = rand5Power11()
        if (r < bMax) {
            // 在 7^9 范围内的数，转换为 9 位 7 进制数，添加到随机数池
            supplyPoolWith(r, exp7)
        } else {
            // 在 7^9 范围外的数，取其超出的部分，转换为 8 位 7 进制数（5^11 - 7^9 > 7^8）
            val rem = r - bMax
            if (rem < 7.power(exp7Rem)) {
                supplyPoolWith(rem, exp7Rem)
            }
        }
    }
    
    while (pool.isEmpty()) {
        supplyPool()
    }
    return pool.removeFirst()
}

优化算法效率之路

算法的效率极限

将 rand5() 调用 $m$ 次产生的随机数序列，称为 rand5() 的 $m$ 次扩展序列。其状态空间的大小为 $5^m$，其中每个状态都是等概率的。

相似地，目标函数 rand7() 的 $n$ 次扩展序列，具有均匀概率分布的 $7^n$ 大小的状态空间。

对于可作用于 rand5() 的任意确定性算法，假设其接受的源序列长度为 $m$，生成的目标序列长度最大为 $n$，根据之前讨论的关于确定性算法的性质，源序列的状态空间大小大于等于目标序列的状态空间大小：

$$\begin{aligned} 5^m &\geq 7^n \\ \frac{m}{n} &\geq \log_5 7 \approx 1.209 \end{aligned}$$

于是，我们得到算法的效率极限：$\log_5 7$。

构造极限效率的算法：算术编码

我想用一个小故事来引入接下来的主角：

一个外星人访问地球时，对人类的文化和知识产生了浓厚的兴趣。

在他准备离开时，陪同的科学家们不禁表达了留恋之情：“您还未完全体验我们丰富的文化和艺术，为何不延长逗留时间呢？”

外星人回答说：“非常感谢你们的热情款待，但在我短暂的访问期间，我已经将地球上所有的文化知识都记录下来了，它们都储存在这根金属棒上的刻痕中。”

科学家们仔细观察那根金属棒，确实在中间部分发现了一条细微的刻痕，但除此之外看不出有什么不同。
他们惊讶地问：“这么一小条刻痕怎么可能储存起人类历史上如此庞大的信息量呢？”

外星人解释说：“其实我们和你们一样，也是使用二进制来存储信息。”
“我把所有要保存的信息转换成二进制数字，并在前面加上一个零和小数点。”
“然后，我把金属棒的一端定义为 0，另一端定义为 1，在代表这个二进制小数的确切位置上刻上一个标记。”
“这样一来，仅用一个数字，我就能保存我所需的全部信息。”

理论上，实数可以存储任意数量的信息，但在物理世界中，由于保存每位数字的精度要求指数上升，这一方法会很快失去可操作性。
假设外星人有普朗克尺度的精度，和可观测宇宙那么长的棍子，那么它可以用上面的方法储存大约 200 比特的信息。

尽管外星人的信息保存法不具有可操作性，但其使用的编码方式，提供了构造接近极限效率算法的思路。

非常粗略地说，算术编码就是将信源序列映射到 $[0, 1)$ 区间的子区间，然后取子区间内的一点作为码字。算数编码的算法，可以逐位递推地计算目标区间范围，而不需要使用整个序列的输入，从而可以实现连续地编码。

篇幅所限，这里不打算介绍算术编码的详细知识，只看如何应用到构造随机数的问题上。

假设我们有对实数进行储存、计算的设备，我们可以使用 rand5() 生成无限长度的序列，并将其映射到 $[0, 1)$ 区间的一个实数上。
$x_n$ 表示生成的序列的第 $n$ 位数字，实数 $R$ 的构造方法为：

$$R = (0.x_1 x_2 x_3 \cdots )_5 = \sum_{n = 1}^{\infty} 5^{-n} x_n$$

这样的操作相当于从 $[0, 1)$ 区间随机采样一个实数。可以证明，对区间内的任意一个子区间 $[a, b)$，采样的实数落在其中的概率严格等于 $b - a$。（证明方法留给读者作为练习）

接下来，可以从这个实数解码出符合要求的 rand7() 序列，基本上是上面构造方法的逆过程。$y_n$ 表示目标序列的第 $n$ 位数字，$Y_n$ 表示目标序列前 $n$ 位组成的七进制小数：

$$Y_n = (0.y_1 y_2 \cdots y_{n})_7$$

则 $y_n$ 可以由 $R$ 和 $Y_{n-1}$ 求得，$Y_n$ 可以由 $y_n$ 和 $Y_{n-1}$ 求得。令 $Y_0 = 0$，可通过递推方式求得目标序列的每一位：

$$\begin{aligned} y_n &= \lfloor 7^n (R - Y_{n-1}) \rfloor \\ Y_n &= Y_{n - 1} + 7^{-n} y_n \end{aligned}$$

实际的计算机系统不能直接对任意精度的实数进行储存和计算，但观察上面的递推公式，可以发现只要给出与实数 $R$ 足够接近的近似数，就不会影响 $y_n$ 的计算结果。

对我们的例子，近似数 $X_m$ 由 rand7() 的 $m$ 次扩展序列给出：

$$X_m = (0.x_1 x_2 \cdots x_{m})_5$$

得到了 $X_m$，就限制了 $R$ 所在的区间范围是 $[X_m, X_m + 5^{-m})$。
使用该限制的上下界代替 $R$ 到递推公式中，如果两者得到的 $y_n$ 一致，代表截止到第 $n$ 位的精度是准确的，不会影响输出结果；如果不一致，就继续扩展输入序列，进行比较，直到满足条件为止。

由此给出下列 Python 风格伪代码实现：

Yn = 0
Xm = 0
m = 0
while True:
    m += 1
    Xm += rand5() * 5**(-m)
    Xmm = Xm + 5**(-m)
    yn1 = floor(7**n * (Xm - Yn))
    yn2 = floor(7**n * (Xmm - Yn))
    if yn1 == yn2:
        Yn += yn1 / 7**n
        yield yn1 # 输出 rand7 的一位结果

将上述方法应用到实际计算中，仍然需要经过以下优化迭代：

1. 使用整数代替浮点数计算：计算机表示实数使用的浮点数的精度是有限的，随着计算进行，Xm 和 Yn 要求的精度会很快超过浮点数的精度上限。此外，计算机进行浮点数的计算通常慢于整数计算。可以通过改进算法，使用整数代替浮点数计算。
2. 位数过长时进行截断：即使采用整数进行计算，且采用的数学库可以处理任意大的整数（内存空间允许的情况下），整数运算的性能也会随着位数的增长而快速下降。这样的算法虽然保证了最大的随机数利用率，但从计算效率的角度来说不是最优的。通常会在位数增长到特定限值时进行截断，然后从初始状态重新开始过程，以略微损失随机数利用率的代价保证计算效率。

下面给出使用 Kotlin 语言和 java.math.BigInteger 计算的算法，截断上界可以自由调整：

import java.math.BigInteger

private class IntegerOverflowException(): ArithmeticException()
/**
 * 截断上界，乘法运算中成员超过该数值时，抛出 IntegerOverflowException
 */
private val maxBigInteger: BigInteger = (BigInteger.TWO * (BigInteger.ONE + BigInteger.valueOf(Long.MAX_VALUE))).pow(64)
@Throws(IntegerOverflowException::class)
fun BigInteger.timesBounded(other: Int): BigInteger {
    if (this > maxBigInteger) throw IntegerOverflowException()
    return this.times(BigInteger.valueOf(other.toLong()))
}
@Throws(IntegerOverflowException::class)
fun BigInteger.times5(): BigInteger = this.timesBounded(5)
@Throws(IntegerOverflowException::class)
fun BigInteger.times7(): BigInteger = this.timesBounded(7)

class Generator(private val rand5: () -> Int) {
    /** S_{m, n} = (X_m - Y_{n-1}) * 5^m * 7^n ≈ (Y_n - Y_{n-1}) * 5^m * 7^n */
    private lateinit var state: BigInteger
    /** 5^m */
    private lateinit var powerOf5: BigInteger
    /** 7^n */
    private lateinit var powerOf7: BigInteger
    private var m: Int = 0
    private var n: Int = 1

    fun get(): () -> Int {
        // 进行首次初始化
        reset()
        return ::gen
    }

    /**
     * 初始状态，整数溢出时也会重置到该状态
     */
    private fun reset() {
        state = BigInteger.ZERO
        powerOf5 = BigInteger.ONE
        powerOf7 = BigInteger.valueOf(7L)
        m = 0
        n = 1
    }

    /**
     * 拦截溢出异常，重置状态，重新生成
     */
    private fun gen(): Int {
        return try {
            genInner()
        } catch (e: IntegerOverflowException) {
            reset()
            gen()
        }
    }

    @Throws(IntegerOverflowException::class)
    private fun genInner(): Int {
        while (true) {
            // y_{n, low} = (X_m - Y_{n - 1}) * 7^n = S_{m, n} / 5^m
            val low = (state / powerOf5).toInt()
            // y_{n, high} = (X_m - Y_{n - 1} + 1) * 7^n = (S_{m, n} + 7^n) / 5^m
            val high = ((state + powerOf7) / powerOf5).toInt()
            if (low == high) {
                if (low !in 0..6) {
                    println("invalid: $low, $powerOf5, $powerOf7, $state")
                    throw IntegerOverflowException()
                }
                // S_{m, n + 1} / 5^m = 7 * (S_{m, n} / 5^m - y_n)
                state = (state - powerOf5.timesBounded(low)).times7()
                powerOf7 = powerOf7.times7()
                n++
                println("<- $low S_{$m, $n}: " +
                        "5:${(state / powerOf7).toString(5)}, " +
                        "7:[${(state / powerOf5).toString(7)}, " +
                        "${((state + powerOf7) / powerOf5).toString(7)}]"
                )
                return low
            } else {
                // S_{m + 1, n} / 7^n = S{m, n} / 7^n * 5 + x_m
                val rand = rand5()
                state = state.times5() + powerOf7.timesBounded(rand)
                powerOf5 = powerOf5.times5()
                m++
                println("-> $rand S_{$m, $n}: " +
                        "5:${(state / powerOf7).toString(5)}, " +
                        "7:[${(state / powerOf5).toString(7)}, " +
                        "${((state + powerOf7) / powerOf5).toString(7)}]"
                )
            }
        }
    }
}

使用 4096 位无符号整数最大值作为截断上界，平均生成 $4096/log_2 7 \approx 1459$ 个随机数时会遭遇一次截断。
得到的平均调用次数为 $1.210$ 次，非常接近 $log_5 7 \approx 1.209$ 的理论下界。

感兴趣的同学，可以使用下面的驱动代码进行测试：

val nTimes: Int = 100000
var nCount: Int = 0
val map = mutableMapOf<Int, Int>()
repeat(7) {
    map[it] = 0
}
val rand7 = Generator { Random.nextInt(5).apply { nCount++ } }.get()
repeat(nTimes) {
    val value = rand7()
    map[value] = map[value]?.plus(1)?:1
}

println("$nCount / $nTimes = ${nCount / nTimes.toFloat()}")
println(map)

参考链接

1. Expand a random range from 1–5 to 1–7: Adam Rosenfield’s Answer
2. 《信息论基础(第2版)》田宝玉 杨洁 贺志强 许文俊 |微信读书
3. 算数编码的原理及C++实现