C 語言隨機數生成器的實現分析

盧亮發表於2016-12-21
隨機

本文是為了記錄和澄清一個由來已久的關於C語言隨機數生成器的誤解。

目前所看到的所有公開的關於C隨機數生成器的中文資料，都提到經典的線性同餘法( LCG, linear congruential generator)，並認為是預設的實現方法。這個說法並不準確。以GCC為例，GLIBC的確實現了線性同餘法，但是實現的程式碼塊分支在日常使用中不會執行到，線性同餘法為C語言預設隨機數生成器的說法已過時。

本文將以GLIBC原始碼為例，結合掌握的文件，做一回搬運工，總結描述一下GLIBC中隨機數生成器的實現。

線性同餘法

線性同餘法，LCG(linear congruential generator)，是經典的偽隨機數產生器演算法，速度快，容易理解實現。 LCG 演算法數學上基於公式：

X(n+1) = (a * X(n) + c) % m

其中，各系數為：模m, 係數a, 0 < a < m，增量c, 0 <= c < m，原始值(種子) 0 <= X(0) < m 。其中引數c, m, a比較敏感，或者說直接影響了偽隨機數產生的質量。

GLIBC中對LCG的實現，取a = 1103515245, c = 12345, m = 134217728，即

X(n+1) = (1103515245 * X(n) + 12345) & 2147483647

由於LCG計算簡單，極省記憶體，很適合記憶體和計算資源比較緊張的嵌入式環境。但LCG有一個嚴重的缺陷，即產生的偽隨機數強依賴於上一次生成的隨機數，且重複週期等於隨機範圍，不能用於隨機數要求高的場合。

原因是單狀態生成器在每次rand()呼叫時不會生成完全的偽隨機數，實際做的是以偽隨機的順序遍歷(0~2^31)範圍內的數。意味著當獲取到一個為偽隨機數時，在當前週期內不會再獲取到同一個數，只有在經過2^31次rand()呼叫之後，才會獲取這個數（而且只會獲取到這個數）。

線性累加反饋法

線性累加反饋法，即LAFM（linear additive feedback method），以下是GLIBC使用的線性累加反饋法的流程描述。其中，2147483647 = 2^31 – 1，4294967296 = 2^32. 所有變數都是整數。 對於給定的種子常量s, 初始化序列r0…r33通過以下步驟計算：

  1. r(0) = s
  2. r(i) = (16807 * r(i-1) ) % 2147483647 (i = 1…30)
  3. r(i) = r(i-31) (i = 31…33)

注意數乘16807的結果應該由足夠大的整數型別儲存，避免取模操作之前發生值溢位。r(i-1)在乘法操作已經是32位整數，r(i)計算結果確保是0到2147483646之間的正整數, 即使r(i-1)為負數。

從r34開始的偽隨機序列，通過以下的線性反饋迴圈來計算：

4. r(i) = (r(i-3) + r(i-31)) % 4294967296 (i ≥ 34)

忽略掉r0…r343序列，rand()函式輸出的偽隨機數o(i)為：

5. o(i) = r(i+344) >> 1

r(i+344)的個位數字移除，生成31位隨機數o(i)。

以下為模擬步驟1~4的程式碼：

GLIBC的實現

GLIBC實現了以上兩種演算法。LAFM生成器標記為 TYPE1, TYPE2, TYPE_3 和 TYPE4 型別，LCG 生成器標記為 TYPE0。相比LCG，LAFM生成器預先生成有很多初始狀態，消除了LCG生成器的週期性遍歷的屬性，在同一個週期內，可以多次獲取到相同的隨機數。

為了提高隨機數生成的時間和空間效率，在計算偽隨機序列時GLIBC使用指標指向包含前驅隨機值的陣列，寫法與按上述公式步驟直譯的方式有所不同。

具體使用哪個生成器依賴於初始狀態集合，由initstate()函式生成：

總結

LCG生成器在狀態陣列(buf->state)長度為8位元組時才會使用。 狀態陣列長度更大時則會啟用LAFM生成器。通常在使用rand()方法時，會使用srand()設定種子常量，這時狀態陣列預設就是128位元組, 所以實際會啟用LAFM生成器。

最後，以上分析如果有偏差，很可能是作者對資料或程式碼的理解問題，歡迎即時反饋。

參考資料：

  1. 純線性同餘隨機數生成器介紹
    http://www.cnblogs.com/xkfz007/archive/2012/03/27/2420154.html
  2. random/initstate原始碼
    https://github.com/lattera/glibc/blob/master/stdlib/random_r.c
  3. glibc rand function implementation
    http://stackoverflow.com/questions/18634079/glibc-rand-function-implementation
  4. The GLIBC random number generator
    http://www.mscs.dal.ca/~selinger/random/

相關文章