計算機中的偽隨機數序列

02-24

作為隨機選擇的工具，我沒有發現任何優於骰子的東西
—— Francis Galton，1890 年《自然》

1 什麼是偽隨機數序列

如果從字面意思來看，『偽隨機數序列』就是偽造的隨機數序列。為了準確定義，在這裡假設擲骰子可以得到一個真正的隨機數，它的結果符合 1 ~ 6 上的均勻分布（注 1）。如果多次擲骰子，在這個過程中，可能可以看到以下序列

1, 5, 6, 5, 6, 2, 5, 4, 5, 3, 4, 5, 5, 5, 3, 4, 2, 1, 1, 3

也有可能出現這個序列

3, 6, 6, 6, 2, 3, 1, 5, 2, 2, 4, 5, 6, 3, 1, 3, 3, 5, 5, 6

甚至有可能，出現這樣的序列，雖然概率很小

6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6

在不停擲骰子的過程中，在結果出來之前，永遠也不會知道當前這個數是多少。我們只知道在這此擲骰子的結果中某個數可能出現的概率是多少，並且這個概率永遠不會改變。這個性質也可以稱為，獨立同分布。或者轉換一下概念，稱之為不可預測性；相反的，如果只知道當前數字，也不可能確定上一次擲骰子結果。

而『偽隨機數序列』是使用特定的方法，比如通過對某個數學公式的計算，生成一個數字序列，並且讓這個數字序列在各方面看起來儘可能的像一個真正隨機數序列。比如一個最古老的偽隨機數生成演算法——平方取中法（Middle-square method），它的過程如下

$x_n$ 是一個 4 位十進位，比如 1234
計算 $x_{n}^{2}$ ，值為 1522756，補足 8 位，即 01522756
$x_{n+1}$ 等於 $x_{n}^{2}$ 中間 4 位，即 5227

def middle_square(): global rand #偽隨機數生成演算法的內部狀態，也是輸出值 rand = rand ** 2 rand = (rand // 100) % 10000 #也可以考慮轉成字元串再提取 return randrand = 1234 #初始狀態# 將輸出進一步處理，模 6 加 1，獲得[1, 6]範圍內的整數#[middle_square() % 6 + 1 for i in range(0, 20)]

將會輸出序列

2, 6, 3, 1, 4, 1, 2, 5, 3, 5, 5, 1, 6, 4, 1, 5, 5, 1, 6, 4

單獨看這個序列，和前面擲骰子出來結果比較，看起來還是比較像，每個數字都出現了，看似沒有什麼規律性。那麼我們模擬了一個隨機過程——一個偽隨機數序列。

這裡將之為偽隨機數序列，只因為它和真正的隨機數序列比較而言，它實際上是確定的——在平方取中法中，當第一個 rand 確定時，後續序列的所有值也就隨之確定，在這個過程中沒有任何隨機性。

可以使用一個 4 元組 $(Q, sigma, Sigma, f)$ 定義一個偽隨機數生成演算法（在習慣上，也稱之為偽隨機數生成器，Pseudo Random Number Generator，簡稱 PRNG），其中 $Q$ 是有限狀態集合， $sigma$ 為狀態轉移函數， $Sigma$ 是有限輸出集合， $f$ 是從 $Q$ 到 $Sigma$ 的單射。當需要一個偽隨機數序列時，設定初始狀態 $q_0 in Q$ ，每一次輸出都由如下兩個過程組成

$q_{n} = sigma(q_{n-1})$
$output_n = f(q_n), output_n in Sigma$

對於上述的平方取中法，可以直觀的得到 $Q$ 對應 [0, 9999] 中的所有整數， $sigma$ 即平方取中法的計算過程， $Sigma$ 為 [0, 9999] 中的所有整數，而 $f = q_n$ 。在示例上，初始狀態 $q_0 = 1234$

從上述定義來看，當 $q_0$ 確定後，後續所有輸出序列 $output$ 即被確定。事實上，對任何圖靈機模型，或者說任何確定性有限自動機模型，其表達能力是不包括任何真正意義上的隨機——當起始狀態、輸入和所有狀態轉移都被確定時，後續所有過程都將被確定。所以，在不依賴外部物理環境，而單純使用邏輯計算，只能生成這種實質沒有任何隨機性的『偽隨機數序列』（注 2）

值得慶幸的是，我們可以定義良好的 4 元組 PRNG，使得生成序列在某些方面和真正的隨機數序列有一定相似性。而這些相似性，達到一定程度時，就足夠特定目的的使用。某些方面包括

A1. 生成的序列每段之間大概率（統計意義上）不一樣

A2. 各項統計指標，和真正隨機數序列儘可能相似，比如自相關檢測

A3. 即使知道生成序列中的某一段數字，也難以確定之前的數、當前內部狀態 $q_n$ 以及之後的數

A4. 即使知道當前 PRNG 的當前內部狀態 $q_n$ ，也難以確定出之前的狀態（如 $q_{n-1}$ ）和生成的數（如 $output_{n-1}$ ）

A1 - A4 是一個非常概念化的說法，更準確的定義可以參考這篇文章。這 4 條，從上到下，表現得越良好，演算法適用性越廣。一般而言，對於一般的娛樂活動，或者模擬計算（比如蒙特卡洛方法），在一定程度上滿足 A1、A2 條即可。對於加密安全或者高要求的賭博活動，第 A3、A4 條不可或缺，稱之為 CSPRNG（Cryptographically Secure Pseudo Random Number Generator）。

在本文中，並不准備定量分析 A1 - A4 條，對於他們的具體分析過於形式化。綜合 A1 - A4、和實踐中的一些要求，提出以下 3 個更弱化並直觀的要求

B1. 循環周期很長（由 A1 導出）。比如 $frac{1}{7} = 0.142857142857142857...$ ，這個小數點後面的輸出序列的周期為 6，而我們希望 PRNG 的輸出序列具備很長的周期，比如 $2^{32}$ ，甚至更大

B2. 輸出數字序列不相關（由 A2 導出）。相關度越高體現出的隨機性也越弱

B3. 高效，可以很快的給出一個偽隨機數序列。計算機內之所以沒有內置一個擲骰子裝置，是因為太慢了^*^（注 3）

A3 和 A4 涉及的領域在本文中將被忽略，有興趣可以參考 CSPRNG - Wikipedia

2 計數器

def counter(): global state state += 1 return statestate = 0#[counter() for i in range(0, 20)]#[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]

這是一個簡單的計數器實現，可以設置初始值。顯而易見，計數器生成序列，周期非常大，在圖靈機模型中可以任意大（但有限）；在現實計算機中，可以達到該計算機的表示能力上界（主要由存儲能力決定）。它滿足 B1，也只滿足 B1。計數器生成遞增序列和真正隨機數序列差距太大，比如強自相關，我們很難將其作為一個 PRNG 用於實踐。

3 平方取中法

在第 1 節中已經簡單描述過這一方法的過程

def middle_square(): global rand #偽隨機數生成演算法的內部狀態，也是輸出值 rand = rand ** 2 rand = (rand // 100) % 10000 #也可以考慮轉成字元串再提取 return rand

在上述示例中，輸出和內部狀態都是由 4 位整數來表達。這個演算法是馮諾依曼在 19 世紀 40 年代提出了，其本身更多還是玩具性質。這個演算法有幾個嚴重缺陷，讓其很難作為一個 PRNG 來使用

接近 20% 的初始值會讓序列收斂到 0
其餘 80% 的初始值同樣會陷入其他短周期循環中

比如 $q_0 = 1234$ ，其第 56 次迭代即收斂到 0

比如 $q_0 = 9999$ ，其第 3 次迭代即開始了周期為 4 的循環，循環序列為 1600, 5600, 3600, 9600

對於周期很短的序列，其在 B1 的表現就已經非常差。

實際上，可以定義不同長度的 $Q$ 和 $Sigma$ ，比如 $Q$ 的元素是 6 位整數，而 $Sigma$ 的元素是 3 位整數

def middle_square(): global rand #偽隨機數生成演算法的內部狀態，也是輸出值 rand = rand ** 2 rand = (rand // 1000) % 100000 #也可以考慮轉成字元串再提取 return rand % 1000

這裡存在同樣的問題，只是陷入周期循環更慢一些，周期相對而言可能較長而已。這其中的原因有本身 $Q$ 的大小限制了周期的上界，也有其實現內在固有的一些性質。比如對這個實例，一旦 $q_{n} < 1000$ ，就會開始向 0 收斂。

整體來看，平方取中法生成序列只有前面一些數可用，具體有多少數字可用由 $Q$ 的位數決定。而當 $Q$ 位數增大時，就需要考慮位數的限制以及計算複雜度。

4 線性同餘生成器，Linear congruential generator

線性同餘生成器（LCG） 4 元組定義

狀態轉移函數 $sigma$ 為

$q_{n + 1} = (aq_{n} + c) space mod space m$

其中 $a 、c、m$ 是固定值

有限狀態集合 $Q$ 為 [0, m - 1] 中的所有整數， $Sigma$ 為 [0, m - 1] 中的所有整數，而 $f = q_n$ 。

LCG 顯然是非常高效的，但是它的問題在於，對不同的 $a 、c、m$ 取值，特性相差很大，比如

def lcg(): a = 31 c = 5 m = 127 global rand rand = (a * rand + c) % m return rand#rand = 1

這個 LCG 函數對任意 $q_0$ 都將輸出一個周期為 63 的序列，除 $q_0 = 21$ 以外，此時所有輸出都將是 21。很容易發現，LCG 的 $m$ 取值決定了其輸出序列周期的上界。為了直觀感受周期性對 PRNG 模擬的重要性，可以使用一段程序（注 4），根據不同 PRNG 生成偽隨機二值圖

使用 Python Random.random 生成 127 * 127

使用 a = 31，c = 5, m = 127 的 LCG 生成 127 * 127

可以看出，前一幅圖的隨機性表現的更好。而後者，有極強的規律性，這主要是由其短周期導致的。

值得慶幸的是，我們可以需要選擇一個很大的 $m$ 值，同時小心翼翼的選擇特定 $a 、c$ 的值（注 5），以產生一個周期很長並可以達到 $m$ 上界的輸出系列。

比如

def LCG(): a = 69069 c = 1 m = 4194304 # 1 << 22 global rand rand = (a * rand + c) % m return rand#rand = 1

其輸出序列周期是 $m = 4194304 = 2^{22}$

我們來看下這個 LCG 生成的偽隨機二值圖

使用 a = 69069，c = 5, m = 4194304 的 LCG 生成 127 * 127

看上去也挺不錯的。我們很容易證明一個輸出序列周期為 $m$ 的 LCG，其輸出是均勻分布在 [0, m - 1] 上。因為在一個周期內，同一個數只會出現一次，而總共只有 $m$ 個數，所以 [0, m - 1] 範圍內的每一個整數各出現一次。但如果我們把範圍縮小，會發現在周期的不同位置，不同大小的數出現分布有顯著差別。現在我們可以繼續看下大尺寸圖像的效果

使用 a = 69069，c = 5, m = 4194304 的 LCG 生成 1024 * 1024

仔細看圖，會發現有明顯的條帶，這是為什麼——在圖大小為 1024 * 1024 = 1048576 小於該 LCG 的周期 4194304 的情況下。事實上，這種規律性並不是由周期太短（B1）而導致的，而是序列相關（B2）導致——當 m 為 2 的冪時，LCG 輸出值很容易保持高位不變，而連續在低位變化。在這個前提下，當我們把它歸一化，比如除以 $m$ ，容易得到一串連續的 0 或者 1。一個更加著名的例子是 RANDU - Wikipedia

所以很多現代編程語言標準庫中的 PRNG 並不使用 LCG，比如 Python 的 random 庫中使用的是 Mersenne Twister 。即使使用 LCG 作為 PRNG，也會做適當改進，比如 Java 的 Random 庫中 next 方法使用 64 位整型參與計算，最後狀態保留後 48 位（如果返回 int，則返回 48 位中的高 32 位）

5 更多 PRNG

比如上文提到的 Mersenne Twister，Python 使用它作為 Random 庫的 PRNG，於 1997 年提出。其有巨大的周期 $2^{19937}-1$ （B1），雖然我們已經知道長周期並不一定代表一個好的 PRNG，但是短周期必然不好。 Mersenne Twister 改進了樸素 LCG 中的統計特徵方面的缺陷（A2 / B2），且和 LCG 一樣高效（B3）

使用 Mersenne Twister 生成 1024 * 1024

關於 Mersenne Twister 的更多信息，包括實現，參考

A random number generator (since 1997/10)
Mersenne Twister - Wikipedia

Permuted congruential generator（PCG），其包括一系列變種，統稱為 PCG Family，其 2014 年才被提出。PCG 具有極其良好的特性，可以擁有任意大的周期；良好的統計特性（比 Mersenne Twister 更好，注 6）；高效的計算；同時 PCG 的不可預測性也較前面所說的 PRNG 更好（雖然並未達到 CSPRNG 標準）。PCG Family 應該是目前綜合表現最好的 PRNG 之一。

使用 PCG 生成 1024 * 1024 （注 7）

關於 PCG family 的更多信息，參考 PCG, A Family of Better Random Number Generators

其餘更多 PRNG，參考 List of random number generators - Wikipedia

6 參考閱讀

$pi$ 是一個好的 PRNG 么，類似的還有其他無理數

Is pi a good random number generator?mathoverflow.net

一個關於洗牌演算法問題的小故事

When Random Isnt Random Enough: Lessons from an Online Poker Exploitwww.lauradhamilton.com

War3 錄像使用 PRNG 的確定性，來重演整個偽隨機過程

《魔獸爭霸》的錄像，為什麼長達半小時的錄像大小只有幾百 KB？www.zhihu.com

隨機數生成的 Wiki

Random number generationen.wikipedia.org

注 1：事實上擲骰子得到的數並不會嚴格均勻分布，甚至這個過程是不是嚴格意義上的隨機也有很多不同的意見

注 2：初始狀態的確定可能會引入一些隨機性（隨機種子），如果從不同的角度進行理解，PRNG 的一項功能是將初始狀態的隨機性進行擴散

注 3：計算機內部可能內置其他物理設備，用以搜集諸如外部熱雜訊，或一些電磁現象用以生成隨機數，但是相對而言，還是較慢。可以參考 Hardware random number generator - Wikipedia

注 4: 這段圖片生成代碼如下

#pip install imagefrom PIL import Imagelen = 127 #這個長度根據需要調整，建議根據 prng 周期來選擇img = Image.new(1, (len, len)) pixels = img.load()for i in range(img.size[0]): for j in range(img.size[1]): #將實現 prng 歸一化，或者將 if 條件改下 pixels[i, j] = 0 if prng() > 0.5 else 1img.show()

注 5：關於參數選擇可以參考 Linear congruential generator - Period length

注 6：Mersenne Twister 並沒有通過 TestU01 的所有測試，而 PCG 通過了，TestU01 是一個軟體庫，提供了一系列用於 PRNG 的統計測試

注 7：PCG 實現可以通過 randomstate 庫獲得，randomstate 1.13.3