这是一个经典老题了，如何用 fair coins 来构造一个 biased coin，以及如何用 fair coins 来构造一个 discrete uniform distribution。基于这两个问题，我又 yy 了另外几个小 follow-up，试图去分析其时间复杂度。

Case 1: $\text{Bernoulli}(p)$

先来一个 warm-up exercise：在只有 fair coins （即 $\text{Bernoulli}(0.5)$）的情况下，给定任意一个 $0 < p < 1$，如何从 $\text{Bernoulli}(p)$ 中 sample 一个元素？

这个很好做：我们只要 sample $[0, 1)$ 之间的一个实数 $x$，并和 $p$ 比较就好了。我们考虑逐位 sample 并比较：每次得到 $x$ 的二进制表示第 $i$ 位，$x_i \sim \text{Bernoulli}(p)$，如果 $x_i \neq p_i$，则我们就可以得到 $x$ 和 $p$ 的相对大小关系了，否则继续下一位。

def sample(p):
    while True:
        x = random.randint(2)
        p_i = int(p * 2)
        if x != p_i:
            return int(x < p_i)
        p = p * 2 - p_i

follow-up： 时间复杂度

上面这个算法有最坏时间复杂度保证吗？

显然没有。

那么存在一个最坏时间复杂度保证的算法吗？

对于几乎所有的 $p$ 都不存在。假设存在一个算法能保证 $n$ 步内结束，那么我们考虑用到的所有 $n$ 个 random bit，一共有 $2^n$ 种组合。对于每种组合我们一定会返回一个 0 或者 1，所以算法一定是从 $\text{Bernoulli}(\frac{q}{2^n})$ 中 sample，其中 $q$ 为整数。而几乎所有的 $p$ 都不能表示成 $\frac{q}{2^n}$ 这样的形式。

那上面这个算法的期望复杂度是多少？

显然，每一步有 $\frac{1}{2}$ 的概率直接结束，所以期望就是 2 步结束。挺意外的一个结果，和 $p$ 无关。

Case 2: $U(n)$

令 $n \in \mathbb{Z}^+$，如何用 fair coins 来得到一个 $U(n)$ 的 sample 呢？这里 $U(n)$ 指 $\{0, 1, \dots, n - 1\}$ 上的 uniform distribution。

我们每次 sample $k = \lceil \log_2 n \rceil$ 个 bit，得到一个 $U(2^k)$ 的 sample。如果这个数小于 $n$，那么返回即可，否则重新开始。

def sample(n):
    k = int(math.ceil(math.log2(n)))
    while True:
        x = 0
        for i in range(k):
            x = x * 2 + random.randint(2)
        if x < n:
            return x

follow-up：时间复杂度

上面这个算法有最坏时间复杂度保证吗？

显然没有。

那存在一个有时间复杂度保证的算法吗？

和之前一模一样的 argument，没有。

那这个算法的期望时间复杂度是多少？

每一轮需要 $k$ 个 bit，并且有 $\frac{n}{2^k}$ 的概率结束，所以期望时间复杂度就是 $\frac{k 2^k}{n}$。

follow-up：更优的算法

存在期望复杂度更低的算法吗？

存在。和 warm-up exercise 类似，还是 sample 一个 $[0, 1)$ 之间的实数 $x$，答案即为 $\lfloor xn \rfloor$。我们从高位到低位枚举 $x$ 二进制表示，然后用一个区间 $[l, r)$ 表示 $x$ 可能取值，一旦发现 $\lfloor ln \rfloor = \lfloor r n \rfloor$，那么我们就已经得到了 $\lfloor xn \rfloor$，返回即可。

def sample(n):
    l, r = 0, 1
    while True:
        if random.randint(2):
            l = (l + r) / 2
        else:
            r = (l + r) / 2
        if int(l * n) == int(r * n):
            return int(l * n)

follow-up：避免浮点数操作

上面这个做法需要用到两个浮点数 $l, r$，其精度有限。能避免使用无限精度的数据类型吗？

可以。我们回到一开始的算法，它没有达到最优是因为，如果 $n \leq x < 2^k$，我们只会重来，白白浪费了这里面用到的 randomness。新算法如下：我们维护两个数 $x, m$，表示现在我们有一个 $U(m)$ 的 sample $x$。如果 $m \geq n$，那么我们观察 $x$：

• 如果 $x < n$，那么 $x$ 就是 $U(n)$ 的一个 sample，直接返回即可；
• 否则我们就知道 $x$ 是 $U(n, m)$ 的一个 sample，故 $x' = x - n$ 就是 $U(m - n)$ 的一个 sample。 而每次拿到一个新 bit 后，我们可以得到一个 $U(2m)$ 下的 sample $x' = 2x + \text{Bernoulli}(0.5)$，重复操作直到 $m \geq n$ 为止。

def sample(n):
    x, m = 0, 1
    while True:
        x, m = x * 2 + random.randint(2), m * 2
        if m >= n:
            if x < n:
                return x
            x, m = x - n, m - n

follow-up：时间复杂度分析

上面这个算法期望时间复杂度，$f(n)$，具体为多少？

显然对 $(x, m)$ 整一个高斯消元是可行的，但是这有点复杂了。聪明一点的同学可以考虑 $f_m$ 表示在已经有了一个 $U(m)$ 的 sample 后得到 $U(n)$ 的 sample 的期望复杂度，我们有：

$$f_m = \begin{cases} \frac{m-n}{m} f_{m-n}, & m \geq n \\ f_{2m} + 1 , & m < n\end{cases}$$

这样未知数会少一些……

但是 somehow 我们有更聪明的方法。令 $Y$ 表示几轮之后算法结束，我们先来一个常用技巧：

$$\mathbb{E}[Y] = \sum_{y \geq 0}^\infty \Pr[Y \geq y].$$

于是问题转化为计算 $y$ 步之后算法还没有结束的概率是多少。我们把算法作一点点修改：之前是如果 $x < n$ 那么我们返回 $x$，现在是如果 $x \geq m - n$ 我们返回 $x \bmod n$，在 $x \sim U(m)$ 的情况下这两者是等价的。

def sample(n):
    x, m = 0, 1
    while True:
        x, m = x * 2 + random.randint(2), m * 2
        if m >= n:
            if x >= m - n:
                return x % n
            m = m - n

我们可以观察到，第 $y$ 轮结束之后的 $m$ 刚好就是 $2^y \bmod n$，于是假设有无限精度，再改动一点点点点：

def sample(n):
    x, m = 0, 1
    while True:
        x, m = x * 2 + random.randint(2), m * 2
        if x >= m % n:
            return x % n

然后就可以很清晰的看到，$y$ 轮之后仍然不结束的概率是

$$\Pr[Y \geq y] = \Pr[x < (2^y \bmod n)] = \frac{2^y \bmod n}{2^y}.$$

故算法的期望时间复杂度为

$$f(n) = \mathbb{E}[Y] = \sum_{y=0}^\infty \frac{2^y \bmod n}{2^y}.$$

虽然这个求和有无数项，但是只要注意到 $2^y \bmod n$ 有循环节就好办了。

follow-up：时间复杂度上下界

刚才我们得到了 $f(n)$ 的精确表达式，但是看起来不存在一个对于任何一个 $n$ 均成立的 closed-form。那我们能否得到一个 $f(n)$ 的上下界？

一个显然的，不需要任何推导的下界是：

$$f(n) \geq \log_2 n,$$

因为 $U(n)$ 包含 $\log_2 n$ bit 的信息，所以我们期望至少 $\log_2 n$ 次才能得到一个 sample……

稍微一分析，我们就可以得到一个 $f(n)$ 的更强的下界。令 $k = \lceil \log_2 n\rceil$，而 $[0, k-1]$ 内的 $y$ 均有 $\frac{2^y \bmod n}{2^y} = 1$，故有：

$$f(n) \geq k \geq \log_2 n.$$

那上界呢？我能想到的最好的一个是：我们把 $y$ 分成三组来估计 $\frac{2^y \bmod n}{2^y}$,

• $0 \leq y < k$：很显然，此时 $\frac{2^y \bmod n}{2^y} = 1$，这一组内所有的和为 $k$；
• $y = k$：有 $\frac{2^k \bmod n}{2^k} \leq \frac{2^k - n}{2^k} = 1 - \frac{n}{2^k}$；
• $y > k$：有 $\frac{2^y \bmod n}{2^y} < \frac{n}{2^y}$，这一组内所有的和有上界 $\sum_{y > k}^\infty \frac{n}{2^y} = \frac{n}{2^k}$；

故 $f(n)$ 有上界

$$f(n) = \sum_{y=0}^\infty \frac{2^y \bmod n}{2^y} < k + 1 - \frac{n}{2^k} + \frac{n}{2^k} = k + 1.$$

综上，我对 $f(n)$ 的最好的 bound 就是

$$\log_2 n \leq \lceil \log_2 n \rceil \leq f(n) < \lceil \log_2 n \rceil + 1,$$

前面的 $\log_2 n$ 就是信息论的 lower bound。可以看到，这个算法最多比理论下界多 2 bits，还算是不错的吧……