rng 效率最高的统一随机整数算法是什么？ [复制]

Question

这不是#11766794 “What is the optimal algorithm for generating an unbiased random integer within a range?” 的副本。它的最高投票答案和接受的答案都与 浮点外推法 有关，这甚至不会产生完全一致的整数。这个问题是在询问如何在给定可使用的统一浮点 rand() 值的情况下快速获得统一随机整数的良好近似；我是在 完美 统一随机整数算法的上下文中提出这个问题的，该算法给出了一个真正的random bit generator（确定性或其他方式；无论哪种方式，这个问题都同样适用）。

我要问的是，仅就使用随机比特的效率而言的理论最优性：给定随机比特流，算法是什么在给定范围内生成完全均匀的随机整数的过程中，哪个消耗了它的最少个位数？

例如，CPython 3.9.0's random.randbelow 至少有一个微不足道的低效率——当调用任意 2 的幂（包括微不足道的范围）时，它会浪费一个随机位：

def randbelow(n):
    "Return a random int in the range [0,n).  Returns 0 if n==0."
    if not n:
        return 0
    k = n.bit_length()  # don't use (n-1) here because n can be 1
    r = getrandbits(k)  # 0 <= r < 2**k
    while r >= n:
        r = getrandbits(k)
    return r

虽然通过将“not n”替换为“n <= 1”和“n.bit_length()”替换为“(n-1).bit_length()”很容易解决这个问题，但一点分析表明它还有很多不足之处：

假设正在生成一个在[0, 4097) 范围内的整数：一半 对getrandbits(13) 的所有调用将超过该值：如果第一个位，说，第二个位高，反正会消耗 11 个更多位，并在看起来不需要的时候丢弃它们。所以看起来这个算法显然不是最优的。

今晚一小时内我能想到的最好的算法是以下算法：

def randbelow(n):
    if n <= 1:
        return 0
    k = (n - 1).bit_length() # this is POPCNT for bignums
    while True:
        r = 0
        for i in reversed(range(k)):
            r |= getrandbits(1) << i
            if r >= n:
                break
        else:
            return r

但是，我不是数学家，仅仅因为我修复了效率低下的问题我可以立即看到，这并不能让我相信我只是在一个下午立即发现了最有效的方法统一整数选择算法。

例如，比特是从量子或大气 RNG 服务购买的；或作为多方协议的一部分，其中每个单独的比特生成都需要多次往返；或者在an embedded device without any hardware RNG support... 不管是什么情况，我只是问一个直接的问题：从真正的随机比特流中生成（完美）均匀随机整数的算法是最有效的关于消耗的随机位？ （或者，如果不确定的话， 目前最好的候选人是什么？）

（我在这些示例中使用了 Python，因为这是我本赛季主要研究的内容，但问题绝不是特定于任何语言的，除了 算法本身 必须推广到数字2⁶⁴以上。）

Answer 1

下面的 Python 以精确的算法实现 arithmetic coding。这在计算上变得非常昂贵，但在期望中实现了熵 + O(1) 位，这基本上是最优的。

from fractions import Fraction
from math import floor, log2
import random


meter = 0


def metered_random_bits():
    global meter
    while True:
        yield bool(random.randrange(2))
        meter += 1


class ArithmeticDecoder:
    def __init__(self, bits):
        self._low = Fraction(0)
        self._width = Fraction(1)
        self._bits = bits

    def randrange(self, n):
        self._low *= n
        self._width *= n
        while True:
            f = floor(self._low)
            if self._low + self._width <= f + 1:
                self._low -= f
                return f
            self._width /= 2
            if next(self._bits):
                self._low += self._width


import collections

if __name__ == "__main__":
    k = 3000
    n = 7
    decoder = ArithmeticDecoder(metered_random_bits())
    print(collections.Counter(decoder.randrange(n) for i in range(k)))
    print("used", meter, "bits")
    print("entropy", k * log2(n), "bits")