線性同餘發生器與僞隨機數

本文旨在簡單探索線性同餘發生器的一些原理和特色，不少思路借鑑於TAOCP，若是想要深刻的探索這方面的知識，建議直接閱讀原著。

1、公式化定義與線性同餘序列的週期

在離散數據及其應用中，若是

那麼，稱a模m同餘b（或者稱模m時，a等價於b），能夠記爲

 而線性同餘式就能夠這樣表示：

線性同餘發生器與上面的線性同餘式多少有一些關係。

2.1 公式化定義

按照The Art of Computer Programming,Volume 2^[1]中3.21. The Linear Congruential Method的思路，線性同餘發生器（LCG：Linear Congruential Generators  ）能夠採用以下公式化定義：

其中：

模數m和係數a是這個公式中最重要的參數，如何合理的選擇這兩個參數決定了其產生的線性同餘序列（LCS：Linear Congruential Scquence）：<X>質量的優劣(<X>=X₁,X₂,X₃..X_n...)。

常數c能夠爲0，也能夠不爲0。一般，若是c=0，那麼（2）式也稱做乘法線性同餘發生器（MCG：Multiplicative Congruential Generator），若是c非0，（2）式,則稱做混合線性同餘發生器（Mixed Linear Congruential Generator）。

X₀稱做初始值，也就是所謂的種子seed。

2.2 線性同餘序列的週期

如上的線性同餘發生器產生的線性同餘序列必然會存在一個週期P。在TAOCP中，做者以一個練習的方式提出了這個問題（exercise 3.1-6）。如下經過一種簡單直白但不嚴謹的推理來解釋這個問題。

將上述線性同餘公式抽象爲一個函數f（將X_n映射爲X_n+1），這個函數具備自封閉特性。不難發現，實際上存在如下已知條件：

定義兩個集合：S和T。初始狀態下，集合S包含了從0到m-1全部的m個元素，集合T是一個空集。

現模擬產生LCS的過程，以任意值X0爲參數，產生第一個僞隨機數X₁。其值必然屬於集合S，此時將X₁從集合S移動到集合T。

以X₁爲參數，產生第二個僞隨機數X₂=f(X₁)，此時，X₂有可能屬於集合S，也有可能屬於集合T。

1）若是X₂屬於集合S，那麼此時尚未產生一個週期；

2）若是X₂屬於集合T，此處也就是恰好等於X₁，那麼此時一個週期產生，週期P=1。

更通常地，假設在生成X₁到X_i-1過程當中，每一個數都是在集合S中找到的，則每一個數都從集合S中移動到了集合T，此時兩個集合的狀態爲：

 

而後生成X_i時，

1）若是X_i=f(X_i-1)在集合S中，則未能產生一個週期；

2）若是X_i=f(X_i-1)在集合T中，則一個週期產生，此時週期P<=i-1。

固然週期P也有可能等於m，也就是集合S最終爲空集，集合T容納了0到m-1的全部元素，且f(X_m)=X₁。

所以，從以上推理不可貴知，LCS必然存在一個週期P，且P<=m。

進而不難推斷：

1）若是某LCG產生的隨機序列的週期P小於m，則選取不一樣的初始值X₀產生的LCS可能有不一樣的週期。

2）若是其週期P=m，則即便選取不一樣的X₀，產生的這些LCS具備相同的週期且一定爲P。

例如，對於下面發生器：

若是以初始值X₀=12產生隨機序列爲：

7,6,9,0,7,6,9,0,7,6,9,0,7,6,9,0......

若是以初始值X₀=13產生隨機序列則爲：

8,3,8,3,8,3,8,3,8,3,8,3,8,3......

可是，對於以下的發生器：

不管選取何值爲種子產生的隨機序列，其週期都是16。

 2、關於參數選擇

構造一個表現良好的線性同餘發生器並不是易事，不但考慮其產生的隨機序列的週期、隨機分佈特色，還要考慮計算的效率。

在參數的選擇方面，最關鍵的就是modulus和multiplier的選擇。

3.1 modulus選擇

modulus應該儘量大，這樣纔有可能產生較長的週期。若是計算機的字長爲w位，TAOCP中推薦，應該取m=2^w，或者m=2^w+1，或者m=2^w-1，也能夠取小於2^w的最大的素數，在論文TABLES OF LINEAR CONGRUENTIAL GENERATORS OF DIFFERENT SIZES AND GOOD LATTICE STRUCTURE^[2] 中就採用了這種m值選取的方式。

一般而言，若是取m=2^w，則利用位運算每每會使計算過程更加方便和高效，可是會存在一個問題：產生的隨機序列中各元素的低比特位的隨機性並非很好。

簡單的解釋是，當m=2^w時，對於一個s位的整數Z，Z模m的結果實際上就是Z的比特位中右邊的s-w位的結果。

做者在原文中用了比較形式化的描述來講明這個問題：

假設d是m的一個因子，q爲某一整數，令Y_n知足以下關係：

而後進行以下變換：

不難發現由公式3-1-1實際就是一個線性同餘發生器，產生的隨機序列也具備週期，可是其週期小於等於d。

這裏的<Y>序列實際上對應了原線性同餘序列<X>的低位字節，能夠將序列<Y>理解爲是將<X>的低位單獨抽取出來組成的一個序列。例如若是d=2^4，則序列<Y>的週期最大也就是16，對應了序列<X>中各個元素的低4位比特位的週期最大是16，顯然低位的隨機性並非很好。正是因爲這個緣由，有些平臺的實現會丟棄這些隨機性很差的低比特位，截取高比特位以取得一個比較好的隨機性效果。大多數應用場景中，低比特位並不會對最終用途產生影響，所以選取m=2^w基本可以知足要求，實際上不少平臺也確實都是取m=2^w。

若是m取2^w+1或者m=2^w-1，則不會產生上述問題。

3.2 multiplier的選擇

multiplier的選擇推理過程比較複雜一些。通常來說，應該使LCS的週期儘可能長（最長爲m），而後只使用一個週期內的元素，可是週期長的序列可能並不具備很好的隨機性。

好比以下的線性同餘發生器：

以初始值X₀=3產生隨機序列的結果：

4,5,6,7,0,1,2,3,4,5,6,7,0,1,2,3,4,5,6,7,0,1,2,3......

可見，以上序列的隨機性表現很糟糕。

在TAOCP中證實了以下定理，按照以下方式來選定係數a能夠產生最大週期爲m的LCS。

以上定理代表當c不等0時（c與m互質固然就不可能等於0），有可能產生週期爲m的LCS。

另外一方面，當c=0時，也即：

是否也有可能產生週期爲m的LCS呢？答案是必然不可能。

 一個簡單的反證法：

若是c=0時，產生了週期爲m的LCS，那麼0必然在這個序列中，可是若是0在序列中，必然會致使LCS退化成全0的序列，所以原命題必然不成立。

從另外一個角度考慮：

考慮d是m的一個因子，且X_n是d的倍數，也即X_n=kd，其中k爲某整數。因而有以下推導：

由此可知，X_n+1也是d的倍數，同理，後面的X_n+2，X_n+2也是d的倍數，如今這樣的序列也不是一個號的序列。因此若是要求序列中各個元素分別與m造成互質關係，所以這樣的序列的元素個數實際就是歐拉函數對m求值，也即：

簡單總結幾點便是：

1）模數m應該儘量大，一般至少大於2^30，爲了計算效率，一般會結合計算機的字長考慮選取m的值。

2）若是m選取爲2的冪，也即m=2^w，則選取的a一般應該知足a模8等於5。

3）當參數m和a的選定比較合理時，對於c的選擇約束性不是很強烈，但要保證c與m互質。例如c能夠選擇1或者11。

4）種子seed應該是隨機選取的，能夠將時間戳做爲種子。

由於X_n的低有效位的隨機性表現並非很好，因此在對X_n的隨機特性比較敏感的應用場景中，應該儘可能採用高比特位。事實上，更應該將值X_n/m視爲0到1之間的均勻分佈，而不是直接將X_n視爲0到m-1之間的均勻分佈。因此，若是但願產生0到k-1之間的均勻分佈僞隨機數，應該採用kX_n/M的方式。 在具體構造一個LCG時，不仿參考TABLES OF LINEAR CONGRUENTIAL GENERATORS OF DIFFERENT SIZES AND GOOD LATTICE STRUCTURE^[2]中的表格來選取參數，該文中針對MLCG和LCG給出了若干可供選擇的參數m和a的值：針對MLCG，m取小於2^w的最大質數；針對LCG，m取2^w。

3、JDK中僞隨機數生成

不少平臺都採用以上線性同餘發生器實現了僞隨機數的生成機制，好比下面是常見的平臺實現中對應的參數（from wiki）：

代碼演示以下（from rosettacode）：

package com.demo;

import java.util.stream.IntStream;
import static java.util.stream.IntStream.iterate;
 
public class LinearCongruentialGenerator {
    final static int mask = (1 << 31) - 1;
    
    static IntStream randBSD(int seed) {
        return iterate(seed, s -> (s * 1_103_515_245 + 12_345) & mask).skip(1);
    }
 
    static IntStream randMS(int seed) {
        return iterate(seed, s -> (s * 214_013 + 2_531_011) & mask).skip(1)
                .map(i -> i >> 16);
    }
    public static void main(String[] args) {
        System.out.println("BSD:");
        randBSD(0).limit(10).forEach(System.out::println);
        
        System.out.println("\nMS:");
        randMS(0).limit(10).forEach(System.out::println);
    }
}

View Code

java中java.util.Random類的實現中，發生器的關鍵代碼以下：

    private static final long multiplier = 0x5DEECE66DL;
    private static final long addend = 0xBL;
    private static final long mask = (1L << 48) - 1;


    protected int next(int bits) {
        long oldseed, nextseed;
        AtomicLong seed = this.seed;
        do {
            oldseed = seed.get();
            nextseed = (oldseed * multiplier + addend) & mask;
        } while (!seed.compareAndSet(oldseed, nextseed));
        return (int)(nextseed >>> (48 - bits));//丟棄低比特位，保留高比特位
    }

    public int nextInt() {
        return next(32);
    }

同時能夠看到，上滿實現中特地對比特位進行了截取，丟棄低比特位，保留高比特位。

4、References

一、Donald Knuth,The Art of Computer Programming, Volume 2

二、TABLES OF LINEAR CONGRUENTIAL GENERATORS OF DIFFERENT SIZES AND GOOD LATTICE STRUCTURE 

三、https://en.wikipedia.org/wiki/Linear_congruential_generator

完。

轉載請註明原文出處：http://www.cnblogs.com/qcblog/p/8450427.html