线性同余发生器的参数如何选取?

我们平时所用的伪随机数生成器(PRNGs)主要有两种:线性同余发生器(Linear Congruence Generator)和反馈位移寄存器法(Feedback Shift Register)。

LCG

线性同余发生器是通过这样的递推函数产生随机序列:

x=(a*x+c)%M (x,a,c,M都是非负整数)

这样产生的随机数序列,一定是有周期的,且小于等于M。在实际应用中,当然希望周期越大越好。我在我的笔记本电脑上测试,生成20亿个随机数只需要20秒左右。如果遇上需要大量随机数的程序,很快就会将这个生成器的"随机性"耗尽。

如果以下三个条件都满足,则可以让该随机数序列的周期等于M

  1. c和M互质。
  2. 对于M的任何一个质因子P,a-1能被P整除。
  3. 如果4是M的因子,则a除以4余1。

(这个定理的证明非常复杂,请参考数论的书)

在计算机系统中,M通常选择是2的整数倍,如\( 2^{48} \),\( 2^{32} \)。

M通常由计算机的字长所限,于是在M给定的情况下,如何选择a和c,则很关键。可以证明,如果M是2的整数次幂,如果a和c满足以下条件,则该序列的周期等于M

  1. a除以4余1。
  2. c除以2余1。

(这个证明过程非常显然,请自推一下)

下面假设通过x=(a*x+c)%M得到的数列为\( \{X_n\} \),那么我们希望它满足[0,M-1]之间的等可能分布。那么它的概率密度函数为:

$$ \begin{equation} \rm P(X_n = i ) = \left\{ \begin{array}{ll} \frac{1}{M}, & i=0,1,\dots,M-1 \\ 0, & other \end{array} \right. \end{equation} $$

于是它的数学期望应该为:

$$ E(X_n) = \sum_{i=0}^{M-1}iP\{X_n=i\} = \frac{M-1}{2} $$

方差应该为:

$$ Var(X_n) = E(X_n^2)-(E(X_n))^2 = \frac{1}{M}\sum _{i=0}^{M-1} \left( i- \frac{M-1}{2} \right) ^{2} = \frac{M^2}{12}-\frac{1}{12}$$

下面计算\( X_n \)的1阶自相关系数,

$$ \rho = \frac{E[(X_n - \mu)(X_{n+1} - \mu)]}{\sigma^2} \approx \frac{1}{a} – \frac{6c}{aM}\left(1-\frac{c}{M}\right)$$

因为我们是在模拟随机数生成器,所以自相关系数应当越小越好,所以a应当尽可能大。c相比于a,应当尽可能小。

FreeBSD中的rand函数有两种实现,一种是WEAK_SEEDING。a=1103515245,c=12345,M= \( 2^{31} \)

FreeBSD的另一种rand实现是a=16807,c=0,M=\( 2^{31} -1 \) 。这个后面讲。

java中的java.util.Random和FreeBSD中的lrand48,采用的都是a=25214903917,c=11,M=\( 2^{48} \)。然后将每次迭代结果的高32位向外部返回。如java中的代码大致为:

nextseed = (oldseed * 25214903917+ 11) & ((1L << 48) - 1); //取余的操作通过二进制位与完成。

MCG

当LCG的c=0,它被称为multiplier congruence generator(MCG)。这种情况下不可能达到满周期,即周期一定小于M。且种子不能为0(也不能是M、2M、3M...),否则得到的序列将是全zero。

对于MCG,若\(M= 2^{L} , (L \geq 4) \),当种子为奇数时,若a除以8余3或者余5,则此MCG的周期可达\( 2^{L-2} \)。

为了研究当M是质数的时MCG的周期,需要引用两个数论中的定义:

\( \varphi(m) \)是Euler Phi函数。特别的,当m是质数的时候,\( \varphi(m) = m-1 \)。

阶(order):\( \{d|a^d\equiv 1 \pmod{m} , d \geq 1 \}\)的最小值称为 a modulo m的阶数。

Primitive Roots: 若a modulo m的阶数等于\( \varphi(m) \),则a称为m的Primitive Root。一般来说,m可以有零个或多个Primitive Root。

在MCG中,若a是M的Primitive Root,则此MCG的周期等于T-1。一个具体的例子,当\(M= 2^{31}-1 \),则M具有以下Primitive Root: 16807、 397204094、 764261123、 630360016。就我目前看过的代码中,选16807的居多。

接下来,具体的计算也很有讲究。

首先要说明的是,所有的MCG都可以转换成c不等于0的LCG来计算。

假设我们要计算\( ax \pmod{2^L-g} \),首先令\( k=[\frac{ax}{2^L}] \),k的计算完全可以通过移位操作。然后计算\( ax \pmod{2^L} + kg \),若计算结果大于\(2^L-g\),则再减去\(2^L-g\)即可。

假如我们要计算\( ax \pmod{2^L} \),那么我们可以利用位与代替mod运算,mod最终要利用除法指令来计算,而位与则和加减法一样高效。

最后摘一段leveldb中的代码:

uint32_t seed_;  
uint32_t Next() {  
   static const uint32_t M = 2147483647L;   // 2^31-1
   static const uint64_t A = 16807;

 // We are computing 
 // seed_ = (seed_ * A) % M, where M = 2^31-1 
   uint64_t product = seed_ * A;
   seed_ = static_cast<uint32_t>((product >> 31) + (product & M));
   if (seed_ > M) {
     seed_ -= M;
   }
   return seed_;
 }

它是借助于64位的整数完成了32位的加法、乘法、求余运算。

这里有个技巧。因为M = 2147483647=127773*16807 + 2836 。我们可以用127773把seed分成高位和低位两部分。

seed = hi* 127773 + lo ;

所以 seed * 16807 = hi* 127773 * 16807 + lo * 16807 = hi * (M - 2836) + lo * 16807

上式最右端对M取余后得到 lo * 16807 - hi * 2836 。 但是这个计算结果有可能是负数,此时让它加上M,就回到[0,M)的范围内了。

但是这个技巧其实并不高明。因为把seed拆成hi和lo两部分的时候,需要计算除法。而在现代CPU上,除法的代价是乘法的10-30倍,很可怕。

注意,multiplier congruence generator(c=0)生成的是[1,M-1]上的均匀分布。所以一般来说会把它减去1后再返回,这样得到的就是[0,M-2]上的均匀分布。别看这是个小东西,leveldb的代码不就犯了这个BUG吗?假如你拿它做随机抽样,结果数组中第一个元素永远也抽不中……嗯……嗯……

本文所描述的随机数生成器不适用于与安全相关的用途,如有需求,请参见 http://erngui.com/rng/

此博客中的热门博文

在windows下使用llvm+clang

少写代码,多读别人写的代码

tensorflow distributed runtime初窥