随机数生成（一）：均匀分布

引言

许多应用中都需要用到随机数，如物理仿真、统计采样、密码学、博彩等。随机数一般可以通过两种方法得到。一种是基于物理现象由硬件产生。由此得到的随机数，在产生之前是不可预知的，因此，是真正的随机数。另一种是通过计算机算法产生。通过算法产生的随机数在本质是可以预知，但是在统计上，满足一定的随机性要求，因此一般称作“伪随机数”。
伪随机数要比真正的随机数更容易获取，而且在大多数情况下都能满足应用的要求。我们这里只讨论伪随机数的生成算法，即随机数生成器（Random Number Generator, RNG）。以下除特别指出外，我们提到的随机数都是指伪随机数。关于由硬件生成真正随机数的方法可以参见 这里，一个产生真正随机数的例子请参见 这里。
经常地，随机数在统计上需要满足特定的分布，如均匀分布、正态分布、指数分布等等。一个好的随机数生成器应该都够高效地生成满足待定分布的随机数序列。而产生的随机数序列有多‘随机“，又在在多大程度上符合特定分布，则需要严格的理论证明和统计分析。但这里，首先要记住一点，随机序列的 周期越长越好，虽然周期并不是评介随机数算法的唯一重要标准。
文本介绍符合均匀分布的随机数生成的相关算法。（由于相关算法非常多，因此，这里先介绍常用的若干种，有机会再逐渐补充。）有关随机数产生质量的评估标准及其他分布的随机数生成算法，留待将来介绍。

——————————————————————————————

[a, b)间的均匀分布随机数x，满足 f(x)=1/(b-a)。x是离散的情况下，f(x)代表概率分布函数；x是连续的情况下，f(x)代表概率密度函数。

线性同余法（Linear Congruential Generator, LCG）

线性同余法（Linear Congruential Generator, LCG）可能是最为常见的一种随机数生成算法，许多语言的标准库中的随机函数都是采用这种方法实现的。

线性同余法通过以下递归式产生均匀分布的随机同序列{X}：

其中，a、m为正整数，b为非负整数，a、b都小于m；X0为初始值（种子），是小于m的正整数。LCG产生的随机序列的周期最大为m，通常都会比m小。序列具有最大的周期m，当且仅当：

1.     b和m互质；

2.     m所有质因子的乘积能够整除(a-1)；

3.     (a-1)是4的倍数，若m是4的倍数。

LCG算法产生的随机序列对a、b、m的选取非常敏感；不恰当的选择会使产生的随机数质量很差。

标准C语言中rand()函数的参数设定为：m = 2³¹，a = 1103515245，b = 12345。每次递归过程，只截取X最高的16位作为最终的随机数输出，因为高位周期更长，由此产生的随机数质量会更好。关于部分编程语言中随机函数在用LCG实现时的具体参数设置可参见这里。

LCG非常容易实现，生成速度快，只需要很小的内存来维护状态（通常仅需要4个或8个字节）。因而对于一般性的应用，LCG是首选。但LCG随机数的周期短（32位的随机数数，周期最多仅为 2 ³²），随性机一般。因而对于需要高质量随机数的应用，如蒙特-卡罗算法，由LCG往往并不是理想的选择。

进位相乘法（Multiply-with-carry，MWC）

MWC是一种与LCG在形式上类似的算法。它由如下递归式产生：

其中，cn代表第n代的进位，r是间隔。我们要用r个种子来初使化前n个数，另需要一个种子来初始化cr-1（cr-1 < a ）。

要了解MWC的参数选择，需要一定的数论和群论知识，可以参看Marsaglia（MWC的提出者）的原帖说明和这份文档及文献[4][6]， [8]给出了在给定不同模b的条件下，最佳参数选择及证明。

在不增加计算量的情况下，MWC能够给出比LCG长得多的周期。例如，同在模同为2^32的条件下，MWC可以产生周期约为10^318的序列，而LCG仅为2^32（约为10^10）。在适当的参数下，MWC可以产生周期长达10^453甚至更长的随机序列。

这里可以多提一下MWC的计算量问题。对多数CPU而言，32位的整数乘法在运算过程中会扩展到64位，低32作为最后的结果，而高32位正是真正结果除以2^32的商，即我们这里说的进位c。因此，相比LCG，MWC并不需要额外的运算来得到进位，只需取出高32位即可。

总之，MWC虽然在在理论上稍显复杂，但设计分析上的复杂性在性能上有着惊人的回报。

线性反馈移位寄存器法（Linear feedback shift register, LFSR）

以上介绍的两种方法，都是直接对整数进行操作， 而下面要介绍的LFSR方法的操作对象是位。同LCG相比，LFSR法也只需要占用很少的内存，但却能产生质量更高的随机数。

LFSR是这样一种移位寄存器：它的输入位是它的当前各位状态的线性函数，即Si+1,1=f(Si)，f(·)为关于Si的线性函数。LFSR最常用到了线性函数是对寄存器的某几位进行异或（XOR）操作。

LFSR可以用硬件实现，也可以通过计算机的移位操作用软件实现。有关前者可以参见 这里，关于软件实现我们会在下面介绍。

斐波那契LFSR（Fabonacci LFSR）

F-LFSR利用斐波那契序列生成随机序列，也称作标准LFSR、多对一LFSR或外部异或门LFSR。下图是一个4位的F-LFSR随机数生成示意图。移位寄存器每周期从左向右移一位。各个位从左到右依次编号为b1、b2、b3、b4，则b1是输入位，b4是输出位。我们需要从输入和输出位之外，另外挑选出若干位，连同输出位一起作为反馈位。这里我们选择b3。线性反馈函数选择为各反馈位（包括输出位）从右向左依次作XOR操作（例如，对16位寄存器，如果我们有b2、b4、b5作为反馈位，则最终反馈函数应为b16^b5^b4^b2）。
假如初始位状态为1-0-1-1，则反馈为（1^1=0）作为输入，各位右移，状态变为0-1-0-1。再下一个周期，（1^0=1）作为输入，各位右移，状态变为1-0-1-0。依次类推。
显然各位全为0的状态是不允许出现的，因为它不能转化为其他的状态。因此，n位的LFSR最多可以产生出2^n-1个不同状态，即周期最长可达为2^n-1。

 

4位F-LFSR

选择不同的反馈位会得到不同周期长度的序列。能够取得最长周期的选择并不唯一，但它们需要满足以下条件：

1. 反馈位的个数为偶数，通常2或4个反馈位就已经中够长周期序列了；
2. 反馈位的编号没有公因数。

产生最长周期的反馈位的配置总是成对出现。例如，如果选择[n, A, B, C]作为反馈位，那么，[n, n-C, n-B, n-A]的选择也会产生一个最长的周期序列。
在FPGA硬件上的g实现参见这里。以下代码片断示例了用C语言实现LFSR，代码完成对LFSR的模拟，并得到周期：

#include <stdint.h>
uint16_t lfsr = 0xACE1u; //初始状态
unsigned bit;            //反馈输入
unsigned period = 0;     //生成序列的周期
do {
  /*反馈位： 16、14、13、11*/ 
 bit  = ((lfsr >> 0) ^ (lfsr >> 2) ^ (lfsr >> 3) ^ (lfsr >> 5) ) & 1; //计算反馈输入（s16^s14^s13^s11）  
 lfsr =  (lfsr >> 1) | (bit << 15); //移位，将反馈至输入位
  ++period;
} while(lfsr != 0xACE1u);

伽罗华LFSR（Galois LFSR）

G-LFSR是另一种形式的LFSR，又称为组合LFSR、一对多LFSR或内部内部异或门LFSR。
每个时钟周期，寄存器的各个位都向右移动一位；非反馈位直接移动，而反馈位在于输出位进行异或操作后，再作右移。异或操作的结果是，当输出位为0时，反馈位相当于直接右移，而当输出位为1时，反馈位反转（即1变为0，0变为1）后右移。

 

16位G-LFSR

如果我们采用与F-LFSR相反的编号方式（即最右位输出位编号为1，依次对各位编号），那么当采用相同的配置时，G-LFSR产生的序列同F-LFSR相同。不同于F-LFSR，G-LFSR的各个异或门可以并行执行，因此，不论软件或硬件实现，G-LFSR都更加快速。
一个关于G-LFRS的硬件实现可以参见这里。以下C代码实现了32位LFSR。

#include <stdint.h>
uint32_t lfsr = 1;
unsigned period = 0;
do {
  /* 反馈位: 32 31 29  */
  lfsr = (lfsr >> 1) ^ (-(lfsr & 1u) & 0xD0000001u); 
  ++period;
} while(lfsr != 1u);

Mersenne twister（MT）

MT由日本学者Makoto Matsumoto和Takuji Nishimura提出。MT同LFSR同属于一类称为GFSR（generalizedfeedback shift register）的方法。选择合适的参数，MT方法可以生成周期长达2^19937-1（算法也因为周期是梅森质数而得名）的32位整数，并且在高达623维下满足均等分布[7]。要得到64精度的整数，可以将两个32位的随机数组合。

MT算法需要634个字长的内存空间（在某些硬件下，这可能是一个缺点），运行速度与C语言的rand函数几乎一样快。 

MT算法的具体流程（包括C代码）及参数选择参见[7]，各种语言的实现可以从这里找到：http://en.wikipedia.org/wiki/Mersenne_twister#Implementations_in_various_languages，官方版本：http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html。基于SIMD计算优化的实现：http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index.html 。Python、Ruby、R、PHP、MATLAB中随机函数的默认算法都是MT，而Boost C++ Library、NAG Numerical Library以及GNU的Scientific Library也实现都了用MT算法实现的RNG。GPU上CUDA版本curand（随SDK发布），另一个版本，但不如前者快。 

MT的一个缺点是启动慢，当MT的种子各位含有较多的0时，需要较长时间跳出来，一般用LCG等初始化。有关这个改进可以参考[11]，从http://en.wikipedia.org/wiki/Well_Equidistributed_Long-period_Linear 可以找到改进一些具体实现。

结语

虽然随机数生成是计算机的一项常规任务，但很多情况下并没有引起足够的重视。Knuth曾说过，“查一查你所在机构的每台计算机上安装的函数库，将里面所有的随机数生成函数都用更好的替换掉。尽量不要对你发现的结果太过于吃惊。”确实，许多随机数函数质量不够好，有些甚至可以说是很坏。当然，这并不意味着实现一个更好的版本是一件容易的事，这需要用到“真正的”数学！对于一般性的任务，一个设计良好的LCG方法可能已经足够“随机”了，并且LCG非常高效，不会占用太多的内存空间。对于随机性要求高的应用，MT一般是首要考虑的方法，但需要占用相对较多的内存空间可能是它的一个不足(大多数情况下，这似乎根本不是个问题:)。其他的算法或者更适合于特定的应用，或者更便于在特定的硬件上实现。总之，当我们在使用随机数的时候，我们最好知道自己在做什么。

Refs

[1]Wikipedia:Random Number Generation.

[2]Wikipedia:Linear Congruential Generator.

[3]Wikipedia:Linear Feedback Shift Register.

[4]GeorgeMarsaglia and Arif Zaman: A New Class of Random Number Generators.(1991)

[5]PierreL'Ecuyer: Uniform random number generation.(1994)

[6]RaymondCouture and Pierre L'Ecuyer: Distribution properties of multiply-with-carryrandom number generators.(1997)

[7]MakotoMatsumoto and Takuji Nishimura: Mersenne twister: a 623-dimensionallyequidistributed uniform pseudo-random number generator.(1998)

[8]MarkGoresky: Efficient Multiply-with-Carry Random NumberGenerators with Maximal Period.(2008)

[9]GeorgeMarsaglia: Xorshift RNGs.(2003)

[10]LeeHowes: Efficient Random Number Generation and Application Using CUDA.(GPU Gems3 Chapter 37)

[11]PierreL'Ecuyer, François Panneton and Makoto Matsumoto: Improved Long-PeriodGenerators Based on Linear Recurrences Modulo 2.(2006)

[12]DavidB. Thomas, Lee Howes and Wayne Luk: A Comparison of CPUs, GPUs, FPGAs, andMassively Parallel Processor Arrays for Random Number Generation.(2009)

[13]C. J. A. Bastos-Filho, M. A. C. OliveiraJunior and A. D. Ramos: Impact of the Random Number Generator Quality onParticle Swarm Optimization Algorithm Running on Graphic Processor Units.(2010)