线性同余发生器的参数如何选取？（以JDK和leveldb的代码为例）

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

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

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，应当尽可能小。

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

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

另外一个常用的选择是，a=1103515245，c=12345，M= 2^{32} 。据说glibc至今依然采用的是这个。

还有一大类是c=0的，被称为multiplier congruence generator，这种情况下不可能达到满周期，即周期一定小于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_;
 }

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