# 從模除偏差談亂數分布
> 資料整理： [jserv](https://wiki.csie.ncku.edu.tw/User/jserv)
〈[解讀計算機編碼](https://hackmd.io/@sysprog/binary-representation)〉提到「計算機以模算數 (modular arithmetic) 進行加法和乘法運算……由於電腦的組成是離散系統，可儲存和操作的位元數量有限，因此能夠表達的數值也會有限制」，而模算數在若干應用情境中，可能會導致非預期的偏差，本文藉此出發，聚焦討論 bounded integer generation 中的亂數分布問題。
本文先說明 `rand() % n` 為何會造成 modulo bias，並介紹常見的去偏差方法與實務 API、探討 `rand()` 的侷限，包含執行緒安全性、作業系統 CSPRNG 介面與硬體亂數產生器，及現代 PRNG 的演進，再來從資訊理論角度討論如何用最少隨機位元產生指定範圍的亂數，隨後繼續討論分布轉換、洗牌演算法、蒙地卡羅方法與方案選擇。
## 如何產生 0, 1, 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}$）。
較小的數字出現機率較高，偏差由此而生。不過，當除數 `n` 遠小於 `RAND_MAX` 時，偏差微乎其微；反之，`n` 愈接近 `RAND_MAX`，偏差愈顯著。例如當 `RAND_MAX = 32767` 且 `n = 3` 時，出現機率如下：
* 0 : 10923 次 $\to$ 33.3344%
* 1 : 10923 次 $\to$ 33.3344%
* 2 : 10922 次 $\to$ 33.3313%
偏差僅約 0.003 個百分點，實務上幾乎察覺不到。然而，這個「幾乎察覺不到」的結論僅在特定條件下成立：`RAND_MAX` 遠大於 `n`。當 `n` 接近 `RAND_MAX` 時，偏差會急遽放大。更關鍵的是，即使單次取樣的偏差微小，在需要大量取樣的情境（如蒙地卡羅模擬、密碼學金鑰產生、線上抽獎系統）中，系統性偏差會隨樣本數增加而被統計檢定（如卡方檢定）偵測出來，或被攻擊者利用。因此，後續討論的重點不在於偏差的絕對大小，而在於理解偏差的結構性成因，以及如何從根本上消除它。
Microsoft Windows 上 `RAND_MAX` 定義為 `0x7FFF`，GNU/Linux 一般定義為 `0x7FFFFFFF`。若在 Windows 上對 $[0, 32767]$ 範圍的 `rand()` 輸出直接取模產生 $[0, 9999]$ 區間的亂數，統計出現頻率可得到下圖：
儘管 `rand() % 10000` 可取得此區間內的亂數，其結果並非均勻分布。為何如此？
集合 $\{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` 的樣本空間，記為 $\Omega_1$。
由此可建立從 $\Omega_0$ 到 $\Omega_1$ 的映射 $f$，對任意 $\omega \in \Omega_0$，$f(\omega) \in \Omega_1$。
但映射 $f$ 並非一對一：多個不同的輸入可映射到同一個輸出。對於輸出值 $x \in \Omega_1$，所有滿足 $f(k) = x$ 的 $k \in \Omega_0$ 所構成的集合稱為 $x$ 的原像 (preimage / inverse image)，記作 $f^{-1}(x) = \{k \in \Omega_0 : k \bmod 10000 = x\}$。
關鍵在於：$|\Omega_0| = 32768$ 無法被 $|\Omega_1| = 10000$ 整除。做帶餘除法：
$$
32768 = 3 \times 10000 + 2768
$$
商數為 3，餘數為 2768。這意味著每個輸出值至少有 3 個原像，而前 2768 個值（$0 \leq x \leq 2767$）各多分到一個，共有 4 個原像。
例如 $x = 2767$ 的原像為：
$$
f^{-1}(2767) = \{2767,\; 12767,\; 22767,\; 32767\}
$$
四個元素皆不超過 32767，因此 $|f^{-1}(2767)| = 4$。
而 $x = 2768$ 的原像為：
$$
f^{-1}(2768) = \{2768,\; 12768,\; 22768\}
$$
下一個候選 $32768$ 已超出 $\Omega_0$，因此 $|f^{-1}(2768)| = 3$。同理，$2768 \leq x \leq 9999$ 的所有值皆僅有 3 個原像。
令隨機變數 $X = \texttt{rand()} \bmod 10000$，各輸出值的機率為其原像大小除以 $|\Omega_0|$：
$$
P(X = x) = \frac{|f^{-1}(x)|}{|\Omega_0|}
$$
因此：
$$
P(X = x) = \frac{4}{32768} > P(X = x') = \frac{3}{32768}, \quad x \leq 2767,\; 2768 \leq x' \leq 9999
$$
因此，`rand() % 10000` 產生的亂數必定不是均勻分布。
問題不在於輸出範圍大小，而在於 $|\Omega_1|$ 能否整除 $|\Omega_0|$。若不能整除，分布必然不均。
舉例來說，當你想產生範圍在 $[min, max]$ 的亂數，若滿足以下條件：
$$
(max - min + 1) \nmid (\text{RAND\_MAX} + 1)
$$
> 此處 $\nmid$ 表示無法整除
則使用以下方法產生的亂數必定不是均勻分布：
```c
rand() % (max - min + 1) + min
```
以模擬擲[骰子](https://dict.revised.moe.edu.tw/dictView.jsp?ID=50695)為例：3 個隨機位元可產生 8 個值 ($2^3 = 8$，$8 \bmod 6 = 2 \neq 0$)。若直接 `% 6`，$110_2 = 6$ 和 $111_2 = 7$ 分別映射到 0 和 1，使 0 和 1 的出現機率為 $\frac{2}{8}$，而 2 到 5 各為 $\frac{1}{8}$。
何時 `rand() % n` 會產生均勻分布？當 `RAND_MAX % n == n - 1` 時。例如 `RAND_MAX = 8`、`n = 3`，此時可保證模除結果（`0`, `1`, `2`）機率相等。
最直覺的做法是持續呼叫 `rand()` 直到取得落在 $[0, n)$ 內的值，但當 `n` 遠小於 `RAND_MAX` 時，接受率 $n / (\text{RAND\_MAX}+1)$ 可能極低。
改良版是先計算 `n` 的最大倍數且不超過 `RAND_MAX + 1`，只拒絕超出此倍數的回傳值：
```c
unsigned int limit = (RAND_MAX + 1u) / n * n;
unsigned int x;
do {
x = (unsigned int) rand();
} while (x >= limit);
x %= n;
```
其中 `limit` 必須是 `n` 的整數倍，如此保留下來的樣本數才能平均分配到 `0` 到 `n - 1`。只要 $1 \leq n \leq \text{RAND\_MAX} + 1$，此法的接受率至少為 $\frac{1}{2}$，平均不到 2 次呼叫即可得到結果。當 $n = \text{RAND\_MAX} + 1$ 時，`limit` 會等於 `RAND_MAX + 1`，因此不需拒絕任何值；若 $n > \text{RAND\_MAX} + 1$，則單次 `rand()` 呼叫無法涵蓋整個範圍，必須改用其他策略。
模除（除法）在 CPU 上相對昂貴。Daniel Lemire 提出的「[近乎無除法](https://lemire.me/blog/2019/06/06/nearly-divisionless-random-integer-generation-on-various-systems/)」(Nearly Divisionless) 方法利用 64 位元乘法取代大部分除法：
```c
uint32_t unbiased_bounded_rand(uint32_t range) {
uint64_t multi = (uint64_t) rand() * (uint64_t) range;
uint32_t low = (uint32_t) multi;
if (low < range) {
uint32_t threshold = -range % range;
while (low < threshold) {
multi = (uint64_t) rand() * (uint64_t) range;
low = (uint32_t) multi;
}
}
return multi >> 32;
}
```
大多數情況下只需一次乘法和位移，僅在極少數邊界情況才退化到除法。此方法要求輸入為完整 32 位元的均勻亂數；程式碼中的 `rand()` 應替換為能產生完整 32 位元的來源（如 `arc4random()` 或自行實作的 PRNG）。若僅有 `RAND_MAX = 0x7FFF`（15 位元）的 `rand()`，需要三次呼叫才能填滿 32 位元，例如 `((uint32_t)(rand() & 0x3) << 30) | ((uint32_t)rand() << 15) | (uint32_t)rand()`。但即便位元數正確，LCG 的連續輸出之間存在線性相關，組合後的品質仍不理想，實務上應直接使用 `arc4random()` 等高品質來源。與前述方案相同，這些做法皆屬 PRNG 層級的去偏差技巧，不提供密碼學安全性。
在高效能運算領域（如物理模擬、遊戲伺服器），開發者常需一次產生大量亂數。Daniel Lemire 也將此技巧擴充至 SIMD (Single Instruction, Multiple Data) 架構，可透過 AVX2 或 AVX-512 進一步提升吞吐量。實際基準會依硬體與編譯器而異，但大方向很穩定：若底層 PRNG 品質相近，Nearly Divisionless 方法通常比 `%` 做範圍映射更快，也明顯快於以 CSPRNG 為底層的 `arc4random_uniform`。
### 使用 `arc4random`
使用 [arc4random](https://man7.org/linux/man-pages/man3/arc4random.3.html) 與 `arc4random_uniform()`，是處理模運算偏差的常見高階解法。此 API 起源於 OpenBSD，但自 glibc 2.36 (2022 年) 起已正式納入 GNU/Linux，因此採用 glibc 的現代 Linux 發行版多半已內建，不再需要額外依賴 `libbsd`。以下為 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;
}
```
這樣可保證結果均勻分布，且幾乎不需要多次重試。`arc4random` 底層使用密碼學安全的串流加密產生完整 32 位元亂數，每次呼叫的成本高於簡單 LCG。對安全敏感的應用而言，`arc4random` 與 `getrandom(2)` 等 CSPRNG 介面是正確的選擇。多數現代作業系統的 CSPRNG 已從早期的 RC4 或 SHA-1 改為 ChaCha20 串流加密演算法：OpenBSD 在 2013 年率先將 `arc4random` 底層從 RC4 換為 ChaCha20，Linux 核心自 4.8 版 (2016 年) 起也將 `/dev/urandom` 的內部產生器改為 ChaCha20。ChaCha20 兼具高吞吐量與常數時間執行特性，後者使其免於時序側通道攻擊，這正是這些 API 兼顧速度與安全性的關鍵。
延伸閱讀：
* [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/)（歷史背景：`arc4random` 名稱中的 "arc4" 源自早期使用的 RC4 演算法，現代實作已全面換為 ChaCha20，與 RC4 的安全弱點無關）
## `rand()` 潛藏的問題
[rand()](https://www.man7.org/linux/man-pages/man3/rand.3.html) 和 [srand()](https://man7.org/linux/man-pages/man3/srand.3.html) 的根本問題在於 API 設計：C 標準允許使用者透過 `srand()` 設定種子來重現亂數序列，卻不允許指定所用的演算法。這二個約束同時存在，使得實作者無法在不破壞既有程式行為的前提下改用更好的 PRNG。相對來說，ECMAScript 規範未指定 `Math.random()` 的演算法，也不提供設定種子的介面，因此引擎實作者可自由替換內部 PRNG 而不影響相容性；例如 [V8](https://v8.dev/blog/math-random) 曾將其實作從 MWC (multiply-with-carry) 改為 `xorshift128+` 的變體。需注意的是，即便是 `xorshift128+` 亦非密碼學安全；在 Web 平台中，若需要無偏差且安全的抽獎或洗牌機制，唯一推薦的做法是使用 Web Crypto API 的 `crypto.getRandomValues()`。此外，JavaScript 的 `Number` 型別僅有 53 位元的安全整數空間 ($2^{53} - 1$)，當現代 Web 應用需要產生 64 位元的隨機識別碼 (如 Twitter Snowflake ID) 時，若以 `Math.random()` 搭配位元運算拼湊，高位元會因浮點精度截斷而遺失，導致碰撞機率大幅上升。正確做法是使用 `BigInt` 搭配 `crypto.getRandomValues(new BigUint64Array(1))` 取得完整 64 位元的密碼學安全亂數。FreeBSD 的 `arc4random` 同樣不暴露種子介面，使其得以在不影響使用者的前提下升級底層演算法。
64 位元隨機識別碼的需求也反映分散式系統的 ID 設計問題。[UUIDv4](https://en.wikipedia.org/wiki/Universally_unique_identifier#Version_4_(random)) 依賴 122 位元亂數，在約 $2^{61.24}$ 筆記錄時碰撞機率才達 50%（生日悖論），對多數應用已足夠；但其完全隨機的特性會導致資料庫 B-Tree 索引的插入點隨機散布，造成頁面分裂與快取失效。因此 Twitter Snowflake、[ULID](https://github.com/ulid/spec) 等方案將時間戳置於高位元，使 ID 具備單調遞增特性，兼顧唯一性與索引效能。這些方案的亂數部分仍須使用 CSPRNG，否則同毫秒內產生的 ID 可能因亂數碰撞而重複。
除 API 設計外，`rand()` 還有幾個實務面的缺陷。C 標準未明確定義 `rand()` 與 `srand()` 的演算法及種子初始化方式，即使用相同種子，也無法保證跨實作、跨標準函式庫版本或跨作業系統產生相同序列。若未呼叫 [`srand()`](https://man7.org/linux/man-pages/man3/srand.3.html)，`rand()` 的行為等同於 `srand(1)`，亦即 `rand()` 只能實作為 PRNG，而非真正的 RNG，且同一實作中無論是否呼叫 `srand()`，演算法都必須一致。常見的錯誤是每次需要亂數時都重新呼叫 `srand(time(NULL))`；`time()` 的解析度僅為秒，同一秒內多次呼叫會重設為相同種子而產生重複序列。正確做法是程式啟動時呼叫一次 `srand()`，之後僅呼叫 `rand()`。
另一個常被忽略的限制是 `srand()` 的種子空間。`srand()` 接受 `unsigned` 型別，在多數平台上僅 32 位元（即使 LP64 的 64 位元架構亦然），可選擇的亂數序列最多 $2^{32}$ 種，遠少於許多 PRNG 本身能支援的序列數（如 C++ `mt19937` 的 $2^{19937}$ 種）。
### 執行緒安全性
`rand()` 維護隱含的全域狀態。C 標準未要求 `rand()` 具備執行緒安全性，[POSIX.1-2024](https://pubs.opengroup.org/onlinepubs/9799919799/functions/V2_chap02.html) 亦未將其列入 thread-safe 函式清單，因此在可攜程式中不應假設多執行緒呼叫 `rand()` 是安全的。實務上，各實作的行為並不一致：glibc 的 [rand(3)](https://man7.org/linux/man-pages/man3/rand.3.html) 以互斥鎖保護內部狀態，標記為 MT-Safe，但這也使其在高度並行環境下成為效能瓶頸；其他實作（如 musl libc）則不加鎖，多執行緒同時呼叫可能直接破壞內部狀態。POSIX 曾提供 `rand_r()` 作為執行緒安全的替代方案，但其狀態空間僅 32 位元，品質不佳，已在 POSIX.1-2024 中移除。現代做法是讓每個執行緒維護各自的隨機狀態 (thread-local state)，或直接使用 CSPRNG 介面（見下節）。
### Linux 與 BSD 的 CSPRNG 介面
當應用需要密碼學安全的亂數時，應使用作業系統提供的 CSPRNG 介面，而非 `rand()`。在 Linux 上，核心維護一個以 ChaCha20 為底層的 CSPRNG（參見 [random(7)](https://man7.org/linux/man-pages/man7/random.7.html)），`/dev/random`、`/dev/urandom` 與 `getrandom(2)` 皆為此 CSPRNG 的前端介面。歷史上 `/dev/random` 會在熵池估計值不足時阻塞，`/dev/urandom` 則不會；Linux 核心自 5.6 版移除 blocking pool 後，兩者在 CSPRNG 初始化完成後的行為已無實質差異。需注意的是，在系統啟動早期（CSPRNG 尚未完成初始化時），`/dev/urandom` 仍可能回傳低熵輸出；`getrandom(2)` 預設會阻塞直到 CSPRNG 初始化完成，因此是現代應用的首選介面。在 BSD 系統上，`arc4random` 系列函式提供類似的保證，且不需要開發者管理種子或檔案描述子。
上述 CSPRNG 介面產生的是密碼學安全偽亂數，而非真亂數 (true randomness)；後者依賴實體熵源，核心 CSPRNG 僅在初始化與重新播種時取用這些熵源。真亂數與 CSPRNG 的層級區別，詳見後文「亂數的三個安全層級」。
### 硬體亂數產生器與 CPU 指令
現代 CPU 內建硬體亂數產生器，可作為作業系統 CSPRNG 的熵源之一。x86 架構提供 `RDRAND`（從硬體 CSPRNG 取得處理後的亂數）與 `RDSEED`（從實體熵源取得未經 CSPRNG 處理的種子），Arm v8.5-A 則提供 `RNDR` 與 `RNDRRS` 指令。Linux 核心會將這些硬體來源混入 entropy pool，但並非唯一的熵源——核心同時收集中斷時序、磁碟延遲等多種環境雜訊，以降低對單一硬體元件的信任依賴。
應用程式不應直接呼叫 `RDRAND` 取代 `getrandom(2)`，原因有二。第一，硬體亂數產生器的內部設計不公開，其正確性僅依賴廠商聲明，無法獨立驗證；Linux 核心透過混合多種熵源來緩解這個信任問題。第二，`RDRAND` 在虛擬化環境中可能被 hypervisor 攔截或替換，而 `getrandom(2)` 的混合機制可提供額外防護。簡言之，`RDRAND`/`RDSEED` 是熵的來源之一，不是應用程式該直接使用的介面。
### 線性同餘產生器的侷限
許多 `rand()` 的實作基於線性同餘產生器 (Linear Congruential Generator, LCG)，其遞推關係 (recurrence relation) 為：
$$
y_{i+1} = (a \cdot y_i + c) \bmod m
$$
其中 $y_0$ 為種子、$a$ 為乘數、$c$ 為增量、$m$ 為模數。此遞推關係以前一個狀態 $y_i$ 經由固定的線性運算產生下一個狀態 $y_{i+1}$，整個序列由初始種子 $y_0$ 唯一決定，故稱為「偽亂數」。依 $c$ 是否為零，分為混合型 LCG ($c \neq 0$) 與乘法型 LCG ($c = 0$，即 Lehmer RNG)。
LCG 能否涵蓋所有 $m$ 個值（即達到最大週期 $m$）取決於參數選取。[Hull-Dobell 定理](https://en.wikipedia.org/wiki/Linear_congruential_generator) (1962) 給出充要條件：混合型 LCG 的週期為 $m$ 若且唯若 (1) $\gcd(c, m) = 1$、(2) $m$ 的每個質因數都整除 $a - 1$、(3) 若 $m$ 是 4 的倍數，則 $a - 1$ 也是 4 的倍數。滿足這三個條件時，LCG 的狀態序列恰好走訪 $\{0, 1, \dots, m-1\}$ 中每個值各一次，這是其作為 PRNG 的理論基礎。然而，「走訪所有值」僅保證一維均勻性，不保證高維分布的品質——Marsaglia 在 1968 年證明 LCG 的連續 $d$ 組輸出落在 $d$ 維空間中至多 $(d!\,m)^{1/d}$ 個超平面上（[Marsaglia 定理](https://en.wikipedia.org/wiki/Marsaglia%27s_theorem)），這正是 LCG 的根本局限。
C 標準 (C11 7.22.2.2 Example) 提供的範例實作即為一個簡單 LCG，其內部狀態以 `unsigned long` 型別儲存，在 32 位元平台上等效於：
$$
y_{n+1} = (1103515245 \times y_n + 12345) \bmod 2^{32}
$$
回傳值取位元 16 至 30（即 `(y / 65536) % 32768`）。在 LP64 平台上 `unsigned long` 為 64 位元，內部狀態的模數相應更大，但輸出仍為 15 位元。Microsoft CRT 的 `rand()` 也使用類似的 LCG，參數為 $a = 214013$，$c = 2531011$。
然而，並非所有平台都使用簡單 LCG。glibc 的 `rand()` 實際上與 [random(3)](https://www.man7.org/linux/man-pages/man3/random.3.html) 共用同一個產生器，預設採用 degree-31 的加法回饋產生器 (additive feedback generator)，其內部維護一個 31 個 `long` 的狀態陣列，品質遠優於簡單 LCG。因此，「LCG 低位元週期短」的警告主要適用於使用 C 標準範例實作或 Microsoft CRT 的平台，不適用於 glibc。
以簡單 LCG 為例，其最低位元週期僅為 2（若 $a$ 和 $c$ 皆為奇數，最低位元在 0 與 1 之間交替），因此在這類實作上不應使用 `rand() % 2` 來產生隨機位元，而應使用高位元，如 `(rand() >> 16) & 1`。
LCG 的共同弱點在於其線性結構：給定足夠的輸出值，攻擊者可推導出後續狀態，因此不適用於密碼學。這一弱點不因加法回饋產生器而消除，因為後者本質上仍是線性遞推。
若攻擊者取得連續三個 LCG 輸出 $y_0, y_1, y_2$ 且已知模數 $m$，當 $\gcd(y_1 - y_0,\, m) = 1$（即差值在模 $m$ 下可逆）時，可直接算出乘數 $a \equiv (y_2 - y_1)(y_1 - y_0)^{-1} \pmod{m}$，再回推增量 $c \equiv y_1 - a \cdot y_0 \pmod{m}$。對於滿足 Hull-Dobell 條件且 $m$ 為 2 的冪的常見 LCG，連續輸出的差值必為奇數，因此模逆元必定存在。即便模數 $m$ 也未知，只要有至少四個連續輸出，令差分 $d_i = y_{i+1} - y_i$，再計算 $g_i = d_{i+1} \cdot d_{i-1} - d_i^2$，由於 $d_i \equiv a \cdot d_{i-1} \pmod{m}$，可推得 $g_i \equiv 0 \pmod{m}$，因此 $m$ 整除所有 $g_i$。對多個 $g_i$ 取 $\gcd$ 後，通常可還原 $m$（或其小倍數，再進一步篩選）。這意味著 LCG 在密碼學情境中形同明文。
### 現代 PRNG 演算法：PCG 與 Xoshiro
隨著統計測試（如 BigCrush）的進步，LCG 甚至 Mersenne Twister (MT19937) 的缺陷已逐漸浮現。MT19937 週期極長 ($2^{19937} - 1$)，但有二個根本問題。第一，其 624 個 32 位元字的內部狀態完全可由 624 個連續輸出還原；MT19937 的輸出函式 (tempering) 是可逆的線性變換，攻擊者只需對 624 個已知輸出求解 GF(2) 上的線性方程組，即可重建完整狀態並預測所有後續輸出。Python 的 `random` 模組即使用 MT19937，在 CTF 競賽中經常成為攻擊目標。第二，MT19937 的狀態空間高達 2.5 KB，且在若干線性複雜度測試中會失敗。
此弱點並非 MT19937 獨有：所有基於線性回饋移位暫存器 (LFSR) 的產生器，輸出位元之間皆存在線性關係，給定足夠的已知輸出即可透過 GF(2) 上的高斯消去法還原回饋多項式。$m$ 位元的 LFSR 最大週期為 $2^m - 1$（當回饋多項式為本原多項式時），統計品質尚可，但不適用於任何安全敏感的情境。
現代非密碼學 PRNG 主要有二個方向：[PCG](https://www.pcg-random.org/) (Permuted Congruential Generator) 結合 LCG 的簡潔性與非線性輸出轉換（如 XSH-RR），通過所有嚴苛的統計測試，速度快且記憶體占用小（僅 8 至 16 位元組）；[Xoshiro / Xoroshiro](https://prng.di.unimi.it/) 系列基於 LFSR 變體，為 Julia、Erlang 等語言的預設演算法，具備極高的位元擴充效能，適合高效能模擬。兩者在統計品質上遠超 MT19937，但仍屬非密碼學等級，安全敏感的應用應使用 CSPRNG。
### C++ 的 `<random>` 與現代實踐
C++11 引入 `<random>` 標頭檔，將亂數引擎與分布模型分離，從 API 層面消除以 `%` 做範圍映射的模除偏差，但種子品質與 PRNG 安全性仍需開發者自行把關。
```c
#include <random>
std::random_device rd; // 非確定性種子來源
std::mt19937 gen(rd()); // 32 位元 Mersenne Twister 引擎
std::uniform_int_distribution<> dis(0, 2); // 定義 [0, 2] 均勻分布
int val = dis(gen); // 得到無偏差的 0, 1, 或 2
```
對於浮點亂數，C++ 另提供 `std::generate_canonical<double, 53>(gen)`，專門產生 $[0, 1)$ 區間內具備指定精度位元數的均勻浮點數，避免簡單除法造成的空洞與端點問題（詳見後文「浮點數產生的陷阱」）。現代 C++ 開發應摒棄 `rand()`；若追求更高效率，可考慮 [PCG](https://www.pcg-random.org/) 的 C++ 實作。
### 亂數的三個安全層級
上述討論涉及多種產生器，可依安全性分為三個層級。第一層是統計偽亂數性 (statistical pseudo-randomness)：輸出位元在統計上近似均勻分布，但內部狀態一旦已知即可預測，LCG 與 MT19937 皆屬此類。第二層是密碼學安全偽亂數性 (cryptographically secure PRNG, CSPRNG)：除統計性外，還要求不可預測，亦即給定前 $k$ 個輸出位元，不存在多項式時間演算法能以超過 $\frac{1}{2}$ 的機率預測第 $k+1$ 位元。[Blum Blum Shub](https://en.wikipedia.org/wiki/Blum_Blum_Shub) (BBS) 是經典的理論範例，其遞推式 $x_{i+1} = x_i^2 \bmod n$（其中 $n = pq$，$p \equiv q \equiv 3 \pmod{4}$）的安全性等價於整數分解問題，但實務上速度太慢，通常改用基於 AES-CTR 或 ChaCha20 的串流加密產生器。第三層是真亂數性 (true randomness)：輸出不可重現，依賴實體現象如放射性衰變、熱雜訊或光子到達時間。
混淆層級會造成嚴重的安全漏洞。在 [ECDSA](https://en.wikipedia.org/wiki/Elliptic_Curve_Digital_Signature_Algorithm) 簽章中，每次簽名需要一個不可預測的亂數 $k$。若 $k$ 由 LCG 產生，攻擊者可從二個簽章的 $k$ 值之間的線性關係（$k_2 = a \cdot k_1 + c \pmod{n}$）聯立 ECDSA 方程式，直接還原私鑰。2010 年的 Sony PlayStation 3 事件正是此類攻擊的實例，起因是簽章時使用固定的 $k$ 值。
## 如何用最少的隨機位元產生指定範圍的亂數？
假設存在均勻隨機位元串流 (bitstream) 作為來源，目標是產生 $[0, n)$ 的整數亂數，並盡可能減少消耗的隨機位元數。編碼 $n$ 個等機率結果至少需要 $\log_2 n$ 位元，這是熵的下界。當 $n$ 是 2 的冪時，$\lceil \log_2 n \rceil$ 個位元即足夠，例如 $n = 8$ 只需 3 個位元；當 $n$ 不是 2 的冪時，固定數量的位元無法均勻對應到 $n$ 個結果，必須引入「重試」機制，實際消耗的位元數因此成為隨機變數。以下討論這個權衡的理論基礎與最優演算法。
### 理論背景
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$ 之間。即使是最佳演算法，平均仍可能「浪費」少量位元。
### 最壞情況與無限迴圈的代價
Knuth 與姚期智也指出：任何最優且無偏差的整數產生器，最壞情況下皆可能執行無限次。原因在於，當 $\frac{1}{n}$ 的二進位展開為無限循環小數（$n$ 非 2 的冪），對應的二元樹要不是無限深，要不就在葉節點加入「拒絕」節點。即便平均只用很少的隨機位元，最壞情況仍可能無限重試。
當 $n$ 是 2 的冪時，對應的樹恰好是滿二元樹，無需拒絕節點，也不浪費任何位元。$n = 8$ 只需 3 個位元即可完整編碼；$n = 5$ 則無法以固定位元數均勻編碼，只能藉重試解決。
### Fast Dice Roller 演算法
此最優演算法由 Jérémie Lumbroso 在 2013 年發表於〈[Optimal Discrete Uniform Generation from Coin Flips, and Applications](https://arxiv.org/abs/1304.1916)〉，以「重試事件」保證無偏差。以下是 C 語言改寫的程式碼範例：
```c
#include <stdlib.h>
#include <stdint.h>
/* Returns a random bit (0 or 1).
* This is illustrative code, not a portable implementation.
* On platforms where rand() is a simple LCG, the lowest bits
* have extremely short periods; we extract high bits instead.
* BITS_AVAIL is the number of usable random bits per rand() call,
* i.e. floor(log2(RAND_MAX + 1)). For RAND_MAX = 0x7FFF (32767),
* that is log2(32768) = 15. On GNU/Linux where RAND_MAX = 0x7FFFFFFF,
* set to 31; on a 64-bit PRNG covering [0, 2^64-1], set to 64.
*/
#define BITS_AVAIL 15
int nextbit() {
static uint32_t state = 0;
static int bits_left = 0;
if (bits_left == 0) {
state = rand();
bits_left = BITS_AVAIL;
}
bits_left--;
int bit = (state >> bits_left) & 1;
return bit;
}
/* Generates a random integer in the range [minInclusive, maxExclusive) */
int randomInt(int minInclusive, int maxExclusive)
{
int range = (maxExclusive - minInclusive) - 1;
int x = 1, y = 0;
while (1) {
x = x * 2;
int bit = nextbit();
y = y * 2 + bit;
if (x > range) {
if (y <= range) {
// Accept the generated value
return y + minInclusive;
}
// Rejection: adjust x and y, then retry
x = x - range - 1;
y = y - range - 1;
}
}
}
```
演算法每次擴充一個位元，檢查值是否在可接受範圍內；若不在範圍內，便捨棄並重試。理論上可能無限迴圈，但平均效率接近 $\log_2 n$。須注意，Fast Dice Roller 的保證建立在底層位元來源獨立且無偏差的前提上；若 `nextbit()` 取自品質不佳的 PRNG，整體均勻性也會受影響。
### 最低位元數和理論極限
上述 Fast Dice Roller 屬於最佳演算法，其平均使用位元數介於 $\log_2 n$ 和 $\log_2 n + 2$ 之間。注意此為平均值的上界；在最壞情況下 (即 $n$ 非 2 的冪)，單次產生可能消耗任意多個位元，但每次重試的成功機率至少為 $\frac{1}{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/)
至此可見，不存在同時滿足固定時間、無偏差與最省位元數的解法；這正是亂數產生的根本權衡。上述演算法對機率模擬、蒙地卡羅方法等需要大量高品質亂數的應用尤為重要。
回顧上述各方案，解決模除偏差大致有三個層次。對一般開發而言，直接使用 `arc4random_uniform` 或 C++ `<random>` 即可，不必自行處理偏差。若需要控制底層 PRNG 或追求效能，Lemire 的近乎無除法方法與拒絕取樣都是成熟的選擇。對位元消耗有嚴格預算的情境（如資訊理論研究或硬體亂數產生器設計），則可考慮 Fast Dice Roller 或批次處理技巧。大多數開發者不需要離開第一層。
## 從均勻亂數到分布轉換與應用
前文已回答「如何無偏差地產生 bounded integer」。均勻亂數是幾乎所有隨機演算法的基礎，以下延伸討論如何將均勻亂數轉換為其他分布、如何驗證品質，以及洗牌與蒙地卡羅等應用。
許多應用需要服從特定機率分布的亂數，例如常態分布 (normal distribution)、指數分布 (exponential distribution) 或卜瓦松分布 (Poisson distribution)。這些分布通常以均勻分布亂數為基礎，透過數學變換得到。
### 浮點數產生的陷阱
產生 `[0, 1)` 區間的浮點數（如 `double`）並非簡單地將整數亂數除以 `RAND_MAX`。若以 32 位元整數除以 $2^{32}$，輸出僅有 $2^{32}$ 個可能值，但 `double` 的 53 位元有效尾數在靠近 0 的區域能表示更細微的數值。這種除法使浮點數如同數線上等距的「定點」，在高精度區段（如 $10^{-7}$ 以下），定點間距遠大於 `double` 能表達的最小間隔，形成巨大的「空洞」(holes)，無法取樣到更細微的值。此外，`rand() / (double)RAND_MAX` 在 `rand()` 回傳 `RAND_MAX` 時會產生 1.0，而後續的對數運算（如反函數變換法中的 $\ln(1-u)$）會因此導致 $\ln(0)$ 錯誤。正確做法是取 53 位元亂數 $L$，計算 $L \times 2^{-53}$，確保結果嚴格落在 $[0, 1)$ 內且充分利用 `double` 的精度。
### 反函數變換法
反函數變換法 (inverse transform sampling) 是最基本的分布轉換技巧。其原理為：若 $U$ 服從 $[0, 1)$ 上的均勻分布，且 $F$ 為目標分布的累積分布函數 (CDF)，則 $X = F^{-1}(U)$ 即服從目標分布。
以指數分布為例，其 CDF 為 $F(x) = 1 - e^{-\lambda x}$，反函數為 $F^{-1}(u) = -\frac{1}{\lambda} \ln(1 - u)$。因此，只要從均勻分布取得 $u$，代入即可得到服從指數分布的亂數。
此方法的限制在於：並非所有分布的 CDF 都有封閉形式的反函數，例如常態分布。
### 取捨法
取捨法 (acceptance-rejection method) 適用於 CDF 反函數不易求得的情境。其步驟為：
1. 選擇一個容易取樣的提議分布 (proposal distribution) $g(x)$，並找到常數 $c$ 使得對所有 $x$ 皆有 $f(x) \leq c \cdot g(x)$，其中 $f(x)$ 為目標密度函數
2. 從 $g(x)$ 取樣得到 $x$，再從 $[0, 1)$ 均勻分布取樣得到 $u$
3. 若 $u \leq \frac{f(x)}{c \cdot g(x)}$，則接受 $x$；否則拒絕並重複步驟 2
取捨法的效率取決於常數 $c$ 的大小：$c$ 愈接近 1，接受率愈高。這與前文討論的拒絕取樣 (rejection sampling) 解決模除偏差的原理一脈相承，皆透過丟棄不符合條件的樣本來確保分布正確性。
### Box-Muller 變換
針對常態分布，[Box-Muller 變換](https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform)提供一個精巧的解法。給定二個獨立的均勻分布隨機變數 $u_1, u_2 \in (0, 1)$，可產生二個獨立的標準常態分布隨機變數：
$$
x = \sqrt{-2 \ln u_1} \cos(2\pi u_2), \quad y = \sqrt{-2 \ln u_1} \sin(2\pi u_2)
$$
此變換利用二維常態分布在極座標下的特殊結構，將均勻分布的半徑與角度分量分別對應到常態分布的振幅與相位。
### 組合法
組合法 (composition method) 適用於目標分布可表示為多個較簡單分布的加權混合的情境。若目標密度函數 $f(x) = \sum_{j} w_j f_j(x)$，其中 $w_j \geq 0$ 且 $\sum_j w_j = 1$，則先依權重 $w_j$ 隨機選擇一個子分布 $f_j$，再從該子分布取樣即可。例如 [Rayleigh 分布](https://en.wikipedia.org/wiki/Rayleigh_distribution) 可分解為 $\chi$ 分布的特例來取樣。此方法的關鍵在於找到有效的分解方式，使各子分布容易取樣。
### 離散分布的區間分配法
對於離散型非均勻分布，區間分配法 (interval allocation) 是直觀的實作方式。給定機率向量 $(p_0, p_1, \dots, p_{k-1})$，將 $[0, 1)$ 區間依各機率值切割為 $k$ 個子區間，從均勻分布取樣後判斷落入哪個子區間。
為避免浮點精度問題，實務上常將機率放大為整數。例如 $(0.2, 0.6, 0.2)$ 對應到整數區間 $[0, 20)$、$[20, 80)$、$[80, 100)$，再從 $[0, 100)$ 取樣。此方法的時間複雜度為 $O(k)$；二分搜尋可降至 $O(\log k)$；若採用 [alias method](https://en.wikipedia.org/wiki/Alias_method)，經 $O(k)$ 前置處理後，每次取樣僅需 $O(1)$，代價是額外 $O(k)$ 記憶體。
### 亂數品質的統計驗證
最直接的驗證方式是大量取樣後統計各值的出現頻率，與預期分布比對。例如對均勻分布的 $[0, n)$ 亂數取樣 $N$ 次，預期每個值出現約 $N/n$ 次；若某些值系統性偏高或偏低，即可偵測偏差。
更嚴謹的做法是統計檢定。[卡方檢定](https://en.wikipedia.org/wiki/Chi-squared_test) ($\chi^2$ test) 比較觀察與期望頻率的偏離程度，適合離散分布；[Kolmogorov-Smirnov 檢定](https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test) 比較經驗 CDF 與理論 CDF 的最大距離，適合連續分布。Université de Montréal 開發的 [TestU01](http://simul.iro.umontreal.ca/testu01/tu01.html) 提供 SmallCrush、Crush、BigCrush 等不同強度的測試組合，是目前最廣泛採用的 PRNG 測試套件。
通過統計測試僅表示「未偵測到偏差」，不能證明產生器完美無缺；反之，未通過則可確定存在問題。
### 亂數在洗牌與抽樣中的應用
無偏差的均勻亂數是正確洗牌演算法的基石。[Fisher-Yates shuffle](https://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle) 是標準的原地洗牌演算法，有兩種等價的索引寫法：從末端往前走訪時，每次從 $[0, i]$ 均勻選取索引 $j$；從前端往後走訪時，每次從 $[i, n)$ 均勻選取索引 $j$。兩者的數學效果相同，皆能確保 $n!$ 種排列各以 $\frac{1}{n!}$ 的機率出現。以下以 forward loop 為例說明常見錯誤。
常見的實作錯誤是將 Fisher-Yates shuffle 寫成以下形式：
```c
// 錯誤示範
for (int i = 0; i < n; i++) {
int j = rand() % n; // 模除偏差 + 排列空間錯誤
swap(&arr[i], &arr[j]);
}
```
此實作存在二個缺陷。第一，`rand() % n` 帶有本文討論的模除偏差。第二，更根本的問題在於排列空間：正確的 Fisher-Yates 在第 $i$ 步從 $[i, n)$ 中選取索引，整個過程的決策空間恰好是 $n!$；而上述錯誤版本每步皆從 $[0, n)$ 選取，決策空間為 $n^n$。由於 $n^n$ 無法被 $n!$ 整除（對 $n \geq 3$），即使底層亂數完全均勻，某些排列出現的機率仍永遠高於其他排列。以 $n = 3$ 為例，$3^3 = 27$ 種決策路徑需映射到 $3! = 6$ 種排列，而 $27 \bmod 6 = 3 \neq 0$，因此至少有一種排列的機率與其他排列不同。這是將「亂數品質」與「演算法正確性」混為一談的典型反例：即使解決模除偏差，錯誤的索引範圍仍使排列分布不均。
若演算法設計錯誤（如上述的 $n^n$ 空間問題）或底層亂數產生器存在模除偏差，某些排列的出現機率就會高於其他排列。在線上遊戲的轉蛋機制或抽獎系統中，這種偏差可能導致獎品分布不符合設計者預期。有限獎池適合用洗牌實作；無限獎池則需確保每次獨立取樣的均勻性。
### 從亂數到蒙地卡羅方法
高品質亂數的另一重要應用是蒙地卡羅方法 (Monte Carlo method)：透過大量隨機取樣近似數學計算。經典範例是估算 $\pi$——在單位正方形內隨機撒點，落入內切圓的比例趨近於 $\pi / 4$。
蒙地卡羅積分將此概念推廣：若要計算 $\int_a^b f(x)\,dx$，可從 $[a, b]$ 均勻取樣 $N$ 個點 $x_1, \dots, x_N$，則
$$
\int_a^b f(x)\,dx \approx \frac{b - a}{N} \sum_{i=1}^{N} f(x_i)
$$
此估計量的標準誤差以 $O(1/\sqrt{N})$ 的速率收斂，不受被積函數維度影響，因此在高維積分問題上優於傳統數值方法。
為降低變異數、加速收斂，可採用重要性取樣 (importance sampling)：改從更接近被積函數形狀的提議分布 $q(x)$ 取樣，並以權重 $f(x) / q(x)$ 修正偏差。當目標分布難以直接取樣時，[馬可夫鏈蒙地卡羅](https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) (Markov Chain Monte Carlo, MCMC) 方法（如 Metropolis-Hastings、Gibbs 取樣）透過構造以目標分布為穩態分布的馬可夫鏈，從鏈的軌跡收集樣本，是貝氏推論與統計物理中不可或缺的工具。
### 實作建議：三個值得記住的做法
以下三種做法涵蓋絕大多數情境，無需自行處理模除偏差：
1. C (BSD / glibc)：呼叫 `arc4random_uniform(n)` 即可取得 $[0, n)$ 的無偏差亂數，不需管理種子，底層為 CSPRNG。適用於 OpenBSD、FreeBSD、macOS，以及提供 `arc4random` 系列函式的 Linux 環境（見上方可攜性附註）。
2. C++11 以上：使用 `std::uniform_int_distribution` 搭配適當引擎（如 `std::mt19937` 或 PCG），由標準函式庫處理去偏差邏輯。
```c
std::mt19937 gen(std::random_device{}());
std::uniform_int_distribution<int> dist(0, n - 1);
int val = dist(gen);
```
3. 需要極致效能的情境（如蒙地卡羅模擬、遊戲伺服器）：採用 Lemire 的 Nearly Divisionless 方法，搭配 PCG 或 Xoshiro 等高品質 PRNG，以避免除法開銷。
若應用涉及密碼學（金鑰產生、簽章隨機數 $k$、token 產生），上述 PRNG 方案皆不適用，應改用 `getrandom(2)`（Linux）或 `arc4random`（BSD）等 CSPRNG 介面。
## 延伸閱讀
* [機率研究系列](https://www.prochainsci.com/2020/10/1.html)：涵蓋有限/無限獎池機制、LCG 文獻回顧、分布產生演算法、區間分配法實作與驗證、馬可夫鏈機率
* [抽樣與蒙地卡羅方法](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)：亂數產生、蒙地卡羅積分與重要性取樣、MCMC 方法
* [Randomness 101: LavaRand in Production](https://blog.cloudflare.com/randomness-101-lavarand-in-production/)：Cloudflare 以熔岩燈作為實體亂數來源的工程實踐
* [PCG, A Better Random Number Generator](https://www.pcg-random.org/)
* [Xoshiro / Xoroshiro generators](https://prng.di.unimi.it/)
* [Fast Random Integer Generation in C++](https://lemire.me/blog/2019/06/06/nearly-divisionless-random-integer-generation-on-various-systems/)
Sign in via Google
Sign in via Facebook
Sign in via X(Twitter)
Sign in via GitHub
Sign in via Dropbox
Sign in with Wallet
Connect another wallet