线性同余发生器与伪随机数

本文旨在简单探索线性同余发生器的一些原理和特色，不少思路借鉴于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