# 從模除偏差談亂數分布
> 資料整理: [jserv](https://wiki.csie.ncku.edu.tw/User/jserv)
〈[解讀計算機編碼](https://hackmd.io/@sysprog/binary-representation)〉提到「計算機以模算數 (modular arithmetic) 進行加法和乘法運算...由於電腦的組成是離散系統，可儲存和操作的位元數量有限，因此能夠表達的數值也會有限制」，而模算數在若干應用場景中，可能會導致非預期的偏差，本文藉此出發，談亂數分布的議題。
## 如何產生介於 0 到 2 之間的隨機數
[rand](https://man7.org/linux/man-pages/man3/rand.3.html) 是個偽隨機數生成器 (PRNG)，會在 0 到 `RAND_MAX`（定義於 `stdlib.h`）之間選取一個自然數。倘若你想產生一個介於 0 到 2 的隨機數，假設 `RAND_MAX` 是 10，那麼你可能會使用 `rand() % 3`，但這無法產生機率相等的 0, 1, 2 數值。乍聽之下，似乎很不直覺，我們考慮以下表格：
| `rand()` 的輸出值 | `rand() % 3` |
| --- | --- |
| 0 | 0 |
| 1 | 1 |
| 2 | 2 |
| 3 | 0 |
| 4 | 1 |
| 5 | 2 |
| 6 | 0 |
| 7 | 1 |
| 8 | 2 |
| 9 | 0 |
| 10| 1 |
可見 `0` 出現 4 次（機率 $\frac{4}{11}$），`1` 也出現 4 次（機率 $\frac{4}{11}$），但 2 只出現 3 次（機率 $\frac{3}{11}$）。
這樣就產生偏差，因為較小的數字出現的機率比較高。不過，這種偏差只有在你要取模的數字接近 `RAND_MAX` 時才會非常明顯。例如當 `RAND_MAX = 32767` 且 `N = 3` 時，出現機率如下：
* 0 : 10923 次 $\to$ 33.3353%
* 1 : 10923 次 $\to$ 33.3353%
* 2 : 10922 次 $\to$ 33.3323%
當 `N` 很小時，這種偏差幾乎不會明顯出現。
值得留意的是，在 Microsoft Windows 平台上，`RAND_MAX` 定義為 `0x7FFF`，而在 GNU/Linux 一般定義為 `0x7FFFFFFF`。若我們希望產生 $[0, 9999]$ 區間內均勻分布的隨機整數，若在 Windows 平台上，以 `rand()` 函式產生範圍在 $[0, 32767]$ 之間的隨機數，並統計每個整數在 $[0, 9999]$ 區間內出現的頻率（即每個數字出現的次數除以總數），再將此頻率繪製成圖表，如下：
我們不難發現，儘管使用 `rand() % 10000 + 0` 可取得此區間內的隨機整數，但這種方式產生的結果實際上並非均勻分布。
為何如此？
集合 $\{0, 1, 2, \dots, 32767\}$ 是 [rand](https://man7.org/linux/man-pages/man3/rand.3.html) 函式產生一次隨機數的樣本空間，記為 $\Omega_0$。於是集合 $\{0, 1, 2, \dots, 9999\}$ 則是使用 `rand() % 10000 + 0` 後所產生的隨機數樣本空間，記為 $\Omega_1$。
由此可建立一個從 $\Omega_0$ 到 $\Omega_1$ 的映射 $f$，對任意 $\omega \in \Omega_0$ 而言，都有 $f(\omega) \in \Omega_1$。
但這個映射 $f$ 並非從集合 $\Omega_0$ 到 $\Omega_1$ 的一對一映射，而是存在如下情況：
$$
f(\omega_0) = f(\omega_1) = f(\omega_2) = f(\omega_3)
$$
例如：
$$
f(2767) = f(12767) = f(22767) = f(32767)
$$
但並非所有情況都存在上述四對一的映射，例如對於 $f(9999)$ 而言，僅存在：
$$
f(9999) = f(19999) = f(29999)
$$
這表示，對於不同的 $f(\omega_0)$ 和 $f(\omega'_0)$ 而言，$\omega_0$ 所對應的取值個數會多於 $\omega'_0$。
由於
$$
p(x = f(\omega)) = \sum_{(10000 \times \omega_0 + 0) \leq 32767} p(x = \omega)
$$
因此有
$$
p(x = f(\omega_0)) > p(x = f(\omega'_0)), \quad \omega_0 \leq 2767, \quad 2768 \leq \omega'_0 \leq 9999
$$
也因此，若只能產生 $[0, 32767]$ 範圍內隨機數的函式，使用 `rand() % 10000 + 0` 產生的隨機數必定不是均勻分布。
你或許想問：限制範圍在 $[0, 32767]$ 不就可解決問題嗎？
若樣本空間 $\Omega_1$ 的元素數量無法整除原始均勻分布樣本空間 $\Omega_0$ 的元素數量，則 $\Omega_1$ 的分布一定不會均勻。
舉例來說，當你想產生範圍在 $[min, max]$ 的隨機數，若滿足以下條件：
$$
(max - min + 1) \nmid (\text{RAND_MAX} - \text{RAND_MIN} + 1)
$$
> 此處 $\nmid$ 表示無法整除
則使用以下方法產生的隨機數必定不是均勻分布：
```c
rand() % (max - min + 1) + min
```
於是，只要使用模運算把一組隨機輸出壓縮到一個較小集合時，若該集合大小無法整除原集合大小，就會導致機率分佈不均。例如：
- 用 3 個隨機位元模除 6 (模擬擲[骰子](https://dict.revised.moe.edu.tw/dictView.jsp?ID=50695)，值為 0 到 5)，而 3 個位元能產生 8 個值
- 若直接 `% 6`，也就是 110 $\to$ 6 % 6 = 0，111 $\to$ 7 % 6 = 1，也就是說，0 和 1 出現機率高。
何時 `rand() % n` 會產生均勻分佈？當 `RAND_MAX % n == n - 1` 時。例如 `RAND_MAX = 8`、`n = 3`，此時可保證模除結果（`0`, `1`, `2`）機率相等。
- [ ] 可能方案 1
持續呼叫 `rand()` 直到產生落在範圍內的數值：
```c
int x;
do {
x = rand();
} while (x >= n);
```
缺點是，若 `n` 很小，成功的機率低，可能要重試很多次。
- [ ] 可能方案 2
找出一個可整除 `n` 的最大範圍，再限制 `rand()` 的回傳值落在此範圍內：
```c
int x;
do {
x = rand();
} while (x >= (RAND_MAX - RAND_MAX % n));
x %= n;
```
這能確保每個結果機率相等，效率較好。
### 使用 `arc4random`
使用 [arc4random](https://man7.org/linux/man-pages/man3/arc4random.3.html) 是處理模運算偏差的標準解法，例如 OpenBSD 提供類似以下實作：
```c
/* Calculate a uniformly distributed random number less than upper_bound
* avoiding "modulo bias".
*
* Uniformity is achieved by generating new random numbers until the one
* returned is outside the range [0, 2**32 % upper_bound). This
* guarantees the selected random number will be inside
* [2**32 % upper_bound, 2**32) which maps back to [0, upper_bound)
* after reduction modulo upper_bound.
*/
uint32_t arc4random_uniform(uint32_t upper_bound)
{
rand_state* z = sget();
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 = __rand32(z);
if (r >= min)
break;
}
return r % upper_bound;
}
```
這樣保證結果均勻分佈，且幾乎不需要多次重試。
延伸閱讀:
* [Use secure random generator on macOS](https://github.com/microsoft/mimalloc/pull/390)
* [Attack of the week: RC4 is kind of broken in TLS](https://blog.cryptographyengineering.com/2013/03/12/attack-of-week-rc4-is-kind-of-broken-in/)
## `rand()` 潛藏的問題
[rand](https://www.man7.org/linux/man-pages/man3/rand.3.html) 和 `srand` 的問題在於：
- 無法指定 `rand` 所用的亂數產生器演算法
- 允許透過 `srand` 來初始化該演算法，使隨機序列可重現
以上二點限制實作層級在改進 `rand` 演算法上的自由度，亦即無法輕易改用密碼學等級的隨機數生成器（RNG），或其他更優秀的 PRNG。相對來說，JavaScript 的 `Math.random` 與 FreeBSD 的 `arc4random` 沒有這個問題，因為它們不允許應用程式提供種子以達成「可重現的隨機性」。
正因如此，[V8](https://en.wikipedia.org/wiki/V8_(JavaScript_engine)) JavaScript 引擎將 `Math.random` 的實作換成 `xorshift128+` 的變體，同時仍維持向下相容性。
不僅如此，尚有以下幾點問題：
- `rand` 與 `srand` 的演算法及其種子初始化方式皆未被 C 標準明確定義，亦即就算用相同的亂數種子，也無法保證在不同實作、標準庫版本或作業系統間得到相同的隨機序列。
- 若未在 `rand` 呼叫前先呼叫 `srand`，則 `rand` 的行為類似於預設呼叫了 `srand(1)`。換言之，`rand` 只能被實作為 PRNG，而非真正的 RNG，且該 PRNG 的演算法在相同實作中，無論有無呼叫 `srand` 都必須一致，因此無法達成真隨機性。
- `srand` 所接受的種子值型別為 `unsigned`，其空間通常僅為 16 或 32 位元（即使在採用 64 位元架構的 C 實作中，`unsigned` 也不一定是 64 位元，例如在 LP64 中，`unsigned` 僅用 32 位元）。這意味著透過 `srand` 所能選擇的隨機序列最多只有 $2^n$ 種（其中 n 為 `unsigned` 的位元數），即使底層的 `rand` 演算法本身能產生更多種序列（例如 C++ 的 `mt19937` 可提供 $2^{19937}$ 種）。
另一缺陷是，C 語言標準中，未規定 `rand` 所產生的數值需符合特定的分佈，包括「均勻分佈」甚至是「近似均勻分佈」。這點與 C++ 的 [std::uniform_int_distribution](https://en.cppreference.com/w/cpp/numeric/random/uniform_int_distribution) 和 [std::uniform_real_distribution](https://en.cppreference.com/w/cpp/numeric/random/uniform_real_distribution) 形成鮮明對比，C++ 明確指定具體的偽隨機數產生演算法，例如 [std::linear_congruential_engine](https://en.cppreference.com/w/cpp/numeric/random/linear_congruential_engine) 和 [mt19937](https://en.wikipedia.org/wiki/Mersenne_Twister)。
這些限制使得 `rand` 和 `srand` 難以用於嚴格產生隨機數場景。
## 如何用最少的空間產生指定範圍的亂數？
假設存在均勻隨機位元串流 (bitstream) 的資料來源，我們希望能用它來產生介於 0 到 n (含) 之間的隨機整數，並用最少的位元數，也就是說，不超過 $\lfloor \log_2{n} \rfloor + 1$ 個位元數，例如當 $n = 8$ 時，期望的演算法最多只用 3 個位元。該如何達成？
### 隨機整數產生器的理論背景
以下討論在平均位元使用數上，開發最佳 (optimal) 的隨機整數產生器，其前提是我們擁有一個「真正」的隨機位元產生器，能夠產生獨立、無偏差的隨機位元。
在 1976 年，[Donald Knuth](https://en.wikipedia.org/wiki/Donald_Knuth) 和[姚期智](https://en.wikipedia.org/wiki/Andrew_Yao)在〈The complexity of nonuniform random number generation〉證明：
任何透過隨機位元產生整數（依指定機率分佈）的演算法，都可表述為一棵二元樹:
- 每個隨機位元決定往左或右走
- 每個樹葉節點代表一個結果
他們進一步指出：任何最優的二元樹演算法，用來產生範圍 $[0, n)$ 的整數，其平均所需位元數會介於 $log_2{n}$ 和 $log_2{n} + 2$ 之間。 亦即，即使是最佳的演算法，在平均意義上仍可能會「浪費」一些位元。
### 最壞情況與無限迴圈的代價
[Donald Knuth](https://en.wikipedia.org/wiki/Donald_Knuth) 和[姚期智](https://en.wikipedia.org/wiki/Andrew_Yao)也指出：任何最優且無偏差的整數產生器，在最壞情況下都可能執行無限次，因為當 $\frac{1}{n}$ 的二進位展開為無限不循環小數（例如 n 不是 2 的冪時），則對應的二元樹：
- 要不是無限深
- 要不就是在葉節點上加入「拒絕」節點（代表需重試）
因此，即便平均只用很少的隨機位元，但仍可能在最壞情況下執行無限次。
而當 `n` 是 2 的冪時，對應的樹剛好就是滿二元樹，不需要拒絕節點，也不會浪費任何位元。例如：
- n = 8 $\to$ 只需要 3 個位元就能完全編碼所有結果
- n = 5 $\to$ 則無法用固定數量的位元完全均勻編碼，只能用重試法解決
### Fast Dice Roller 演算法
這個最優演算法由 Jérémie Lumbroso 在 2013 年發表於〈[Optimal Discrete Uniform Generation from Coin Flips, and Applications](https://arxiv.org/abs/1304.1916)〉，它用「重試事件」來保證無偏差。 以下是用 C 語言改寫的 Fast Dice Roller 程式碼範例：
```c
#include <stdlib.h>
#include <stdint.h>
/* Returns an unbiased random bit (0 or 1) */
int nextbit() { return rand() & 1; }
/* Generates a random integer in the range [minInclusive, maxExclusive) */
int randomInt(int minInclusive, int maxExclusive)
{
int maxInclusive = (maxExclusive - minInclusive) - 1;
int x = 1, y = 0;
while (1) {
x = x * 2;
int bit = nextbit();
y = y * 2 + bit;
if (x > maxInclusive) {
if (y <= maxInclusive) {
// Accept the generated value
return y + minInclusive;
}
// Rejection: adjust x and y, then retry
x = x - maxInclusive - 1;
y = y - maxInclusive - 1;
}
}
}
```
該演算法每次擴展一個位元，然後檢查產生的值是否在可接受的範圍內，若不在，就進行捨棄並重新嘗試。雖然理論上可能無限迴圈，但平均效率接近最小可用位元數，亦即 $log_2{n}$。
### 最低位元數和理論極限
上述 Fast Dice Roller 屬於最佳演算法，其平均使用位元數不會低於 $log_2{n}$，且在最壞情況下最多也只會超出 2 個位元。
若要進一步降低位元浪費，該考慮：
- 批次處理（batching）：一次產生多個亂數
- 隨機性提取（randomness extraction）：回收被拒絕掉的位元
延伸閱讀：
- 論文〈[Optimal Discrete Uniform Generation from Coin Flips, and Applications](https://arxiv.org/abs/1304.1916)〉第 3.1 節（batching）
- 論文〈[Random variate generation using only finitely many unbiased, independently and identically distributed random bits](https://arxiv.org/abs/1502.02539)〉，第 2.3 節（randomness extraction）
- [Faster random integer generation with batching](https://lemire.me/blog/2024/08/17/faster-random-integer-generation-with-batching/)
- [Nearly Divisionless Random Integer Generation On Various Systems](https://lemire.me/blog/2019/06/06/nearly-divisionless-random-integer-generation-on-various-systems/)
至此，我們知道，不存在一個同時能在固定時間內且又無偏差且最省位元數的解法 —— 這是產生隨機數的權衡。上述演算法對於密碼學或機率模擬等需要高品質隨機數的應用尤為重要。
## 待整理
* [機率研究 (1) - 有限及無限個數機制與應用情景](https://www.prochainsci.com/2020/10/1.html)
* [機率研究 (2) - 亂數產生器文獻探討](https://www.prochainsci.com/2020/10/2-mathematica.html)
* [機率研究 (3) - 機率分佈生成之 - 應用與分佈的意義](https://www.prochainsci.com/2020/11/3.html)
* [機率研究 (4) - 機率分佈生成之 - 倒推分佈生成以及演算法](https://www.prochainsci.com/2020/11/4.html)
* [機率研究 (5) - 指定機率分佈生成之 - 區間分配法實作與分配驗證](https://www.prochainsci.com/2020/11/5.html)
* [機率研究 (6) - Markov Chain Probability](https://www.prochainsci.com/2020/11/6-markov-chain-probability-designing-in.html)
* [抽樣與蒙地卡羅(一)：隨機數生成](https://medium.com/%E6%95%B8%E5%AD%B8-%E4%BA%BA%E5%B7%A5%E6%99%BA%E6%85%A7%E8%88%87%E8%9F%92%E8%9B%87/%E6%8A%BD%E6%A8%A3%E8%88%87%E8%92%99%E5%9C%B0%E5%8D%A1%E7%BE%85-%E4%B8%80-%E9%9A%A8%E6%A9%9F%E6%95%B8%E7%94%9F%E6%88%90-4c4951fc5523)
* [抽樣與蒙地卡羅(二)：蒙地卡羅方法與重要性抽樣](https://medium.com/%E6%95%B8%E5%AD%B8-%E4%BA%BA%E5%B7%A5%E6%99%BA%E6%85%A7%E8%88%87%E8%9F%92%E8%9B%87/%E6%8A%BD%E6%A8%A3%E8%88%87%E8%92%99%E5%9C%B0%E5%8D%A1%E7%BE%85-%E4%BA%8C-%E8%92%99%E5%9C%B0%E5%8D%A1%E7%BE%85%E6%96%B9%E6%B3%95%E8%88%87%E9%87%8D%E8%A6%81%E6%80%A7%E6%8A%BD%E6%A8%A3-7558764151b9)
* [抽樣與蒙地卡羅(三)：馬可夫鏈蒙地卡羅方法](https://medium.com/%E6%95%B8%E5%AD%B8-%E4%BA%BA%E5%B7%A5%E6%99%BA%E6%85%A7%E8%88%87%E8%9F%92%E8%9B%87/%E6%8A%BD%E6%A8%A3%E8%88%87%E8%92%99%E5%9C%B0%E5%8D%A1%E7%BE%85-%E4%B8%89-%E9%A6%AC%E5%8F%AF%E5%A4%AB%E9%8F%88%E8%92%99%E5%9C%B0%E5%8D%A1%E7%BE%85%E6%96%B9%E6%B3%95-1513fe720ca7)
Sign in with Wallet
Connect another wallet