寫在前面
偽隨機數生成算法在計算機科學領域應用廣泛，比如槍擊游戲里子彈命中擾動、數據科學里對樣本進行隨機采樣、密碼設計、仿真領域等等，背后都會用到偽隨機數生成算法。
說隨機，那什么是隨機呢？隨機意味着不可預測，沒有任何規律。談隨機數，一定是在序列當中，單拿出一個數談隨機是沒有意義的。給一個數字序列，如果能在其中發現規律可以預測或以一定概率（大於“猜”的概率）預測接下來的數，那么這個序列就不是隨機的。
在20世紀早期科學工作中就開始需要使用隨機數，為了獲取隨機數，研究人員通過物理方式采集了成千上萬的隨機數，並發布給他人使用，比如RAND公司在1955年發布的《A Million Random Digits with 100,000 Normal Deviates》（百萬亂數表）——亞馬遜美國現在還有賣~。但是，通過物理方式采集“真”隨機數並不高效，實時獲取需要附加額外的隨機數發生裝置，而且獲取速度緩慢、序列不可復現，如果將采集到隨機數全保存下來則需要占用額外的存儲空間，而且數量終究是有限的，於是大家開始尋求生成“偽”隨機數的數學方法。偽隨機數，顧名思義，即看起來是隨機的但實際上不是，在不知其背后生成方式的情況下，生成的序列看上去毫無規律可言。
本文源自個人興趣通過查閱參考文獻整理所得，再加上個人的理解，大部分圖片來自WIKI。
統計學檢驗
如何判斷一個序列是否夠隨機呢？偽隨機數生成算法多種多樣，總要分出個孰好孰差，如何對各自的隨機性進行定量評估呢？主要有兩類方式，其出發點都是試圖定量評估序列中是否隱含某種規律或模式：
-
實證檢驗。給定一個隨機序列而不告知其背后的生成方式，嘗試對觀測到的分布進行擬合，以預測后面的序列,或者查看其中是否具有某些統計規律，比如查看分布是否均勻、連續值的相關性、某些數出現位置的間隔分布是否有規律等等。具體有\(\chi ^2\)檢驗、KS檢驗、Frequency test、Serial test等等。
-
理論檢驗。直接分析生成器的理論性質（已知生成方式），生成器通常需要配置一些參數，不同的參數會影響生成序列的質量，比如考察不同參數對隨機序列周期的影響。
可在下一小節對理論檢驗一窺一二，但本文的重點不在於此，就不詳細展開了，詳細內容可見參考資料。
線性同余法
linear congruential generator（LCG）線性同余法是最早最知名的偽隨機數生成算法之一，曾被廣泛應用，后逐漸被更優秀的算法替代，其通過如下遞推關系定義：
其中，\(X\)為偽隨機序列，
- \(m\)，\(m > 0\)，模數，顯然其也是生成序列的最大周期
- \(a\)，\(0 < a < m\)，乘數
- \(c\)，\(0 \leq c < m\)，增量
- \(X_0\)，\(0 \leq X_0 < m\)，種子點（起始值）
當\(c = 0\)時，被稱為multiplicative congruential generator (MCG)，如果\(c \neq 0\)，被稱為mixed congruential generator。
線性同余法的參數應該被小心選取，否則生成的序列將非常糟糕，比如當\(a = 11, c = 0, m = 8, X_0=1\)時，得到的序列是 3、1、3、1、3……從1960s開始使用IBM的RANDU算法，參數設置為\(a = 65539, c = 0, m = 2^{31}\)，最終也被證明是糟糕的設計，因為\(65539 = 2 ^{16} + 3\)，可進一步推導得
因為相鄰3個數間存在這樣的相關性，若將相鄰3個數作為坐標繪制在3維坐標系里，會得到15個明顯的平面
可見，獲得的序列並不是那么隨機，而且沒有均勻地填充整個空間。線性同余法的參數很重要，一些平台和運行時庫中采用的參數如下
使用遞推關系的方式帶來了可復現的便利——只需要記住種子點就可以復現整個序列，而不需要去存儲整個序列，但是帶來的弊端就是相鄰點之間的相關性，隨意設置參數（像RANDU）可能讓序列直落在幾個稀疏的平面上，通常需要將\(m\)選取的足夠大，同時避開2的整數次冪。
馬特賽特旋轉演算法
Mersenne Twister 馬特賽特旋轉演算法，是1997年提出的偽隨機數生成算法，其修復了以往隨機數生成算法的諸多缺陷，可快速生成高質量的偽隨機數，且經過了廣泛的統計學檢驗，目前在各種編程語言和庫中已普遍存在或作為默認的偽隨機數發生器，被認為是更可靠的偽隨機數發生器。下圖截自python的官方文檔：
Mersenne Twister生成隨機數的過程比線性同余法要復雜得多，圖示化如下：
主要流程有3步是：
- 初始化\(n\)個狀態：根據給定的種子點\(x_0\)，通過移位、異或、乘法、加法等操作生成后續的\(n-1\)個狀態\(x_1\)到\(x_{n-1}\)，bit位數為\(w\)
- 生成偽隨機數：根據當前狀態，通過移位、與、異或操作生成隨機數
- 更新\(n\)個狀態：每生成\(n\)個隨機數后，在生成下一個隨機數前，更新狀態
具體參見偽代碼（來自WIKI），如下：
// Create a length n array to store the state of the generator
int[0..n-1] MT
int index := n+1
const int lower_mask = (1 << r) - 1 // That is, the binary number of r 1's
const int upper_mask = lowest w bits of (not lower_mask)
// Initialize the generator from a seed
function seed_mt(int seed) {
index := n
MT[0] := seed
for i from 1 to (n - 1) { // loop over each element
MT[i] := lowest w bits of (f * (MT[i-1] xor (MT[i-1] >> (w-2))) + i)
}
}
// Extract a tempered value based on MT[index]
// calling twist() every n numbers
function extract_number() {
if index >= n {
if index > n {
error "Generator was never seeded"
// Alternatively, seed with constant value; 5489 is used in reference C code
}
twist()
}
int y := MT[index]
y := y xor ((y >> u) and d)
y := y xor ((y << s) and b)
y := y xor ((y << t) and c)
y := y xor (y >> l)
index := index + 1
return lowest w bits of (y)
}
// Generate the next n values from the series x_i
function twist() {
for i from 0 to (n-1) {
int x := (MT[i] and upper_mask)
+ (MT[(i+1) mod n] and lower_mask)
int xA := x >> 1
if (x mod 2) != 0 { // lowest bit of x is 1
xA := xA xor a
}
MT[i] := MT[(i + m) mod n] xor xA
}
index := 0
}
標准實現32bit版本稱之為MT19937，參數設置如下：
- \((w, n, m, r) = (32, 624, 397, 31)\)
- \(a = \rm 9908B0DF_{16}\)
- \((u, d) = (11, \rm FFFFFFFF_{16})\)
- \((s, b) = (7, \rm 9D2C5680_{16})\)
- \((t, c) = (15, \rm EFC60000_{16})\)
- \(l = 18\)
后記
偽隨機數生成算法有很多，遠不止本文介紹的兩種，還有middle-square method（1946）、Additive Congruential Method、xorshift（2003）、WELL（2006，對Mersenne Twister的改進）等等，本文只是從中選取具有代表性的兩種，可閱讀參考文獻了解更多。
參考
- Random number generation
- Pseudorandom number generator
- Linear congruential generator
- RANDU
- Randomness tests
- Randomness Tests: A Literature Survey
- Testing Pseudo-Random Number Generators
- Validation of Pseudo Random Number Generators through Graphical Analysis
- How to Generate Pseudorandom Numbers
- NMCS4ALL: Random number generators
本文出自本人博客：偽隨機數生成算法