随机数发生器原理

随机数发生器的基本原理：下一个随机数 = 数学函数(前若干个随机数)

评价一个随机数发生器的好坏：

• 看起来是[0,1]之间的随机均匀分布
• 运行快，内存消耗低
• 可以重复，英文是reproducible：支持stream，即支持将一个随机数序列分割为相互无覆盖的子集，每个子集称为一个stream

线性同余发生器(Linear Congruential Generator)

目前主流的随机数发生器之一，1951年由Lehmer提出。从名字可以看出公式，即一个线性函数再取余：

第i个随机数 = (a*第i-1个随机数 + c)(mod m)

随机数的取值范围是[0,m-1]，将其除以m，就可以得到[0,1)之间的随机数。

• 随机数的取值只能是[0,m-1]中离散的整数值，为了让其除以m后看起来很连续，m需要取较大的值
• 一旦在第i个随机数和之前的某个随机数相同，那么后面的随机数将和之前的序列相同，出现周期性
• 全周期：希望[0,m-1]的每个值在一个周期内都出现一次，从而更有效地表现为均匀的[0,1)分布

全周期定理：满足以下条件时随机数发生序列具有全周期性质

• m和c互质
• 如果q是质数（只能被1和它本身整除）且能够整除m，则它能够整除a - 1
• 如果4整除m，4整除a-1

综上，我们可以得到线性同余发生器的设计方案：m=2^31 （除了让m很大外，还可以有效利用32位计算器integer overflow的性质，避免取模运算；之所以是31不是32，是因为int变量的最左边是符号位，数值表示范围只有31位；此外满足全周期定律第二条，2的次方不能被质数整除）；c为奇数（满足全周期定律第一条）；令a-1能够被4整除（第三条）

如果令c=0，可以避免加法运算，但是不可能获得全周期（不满足全周期定理的第一条，c能够被任意数整除）。但是如果妥善地选择m和a，可以得到m-1的周期长度。

如果令m=2^b，那么理论已经证明，能够获得的最大周期是2^(b-2)，也就是说只能产生[0,m)之间四分之一的数据，且这四分之一的数据如何分布是未知的，很可能显示出来的效果并不是均匀的。换个思路，我们将m设置为小于2^b的最大质数。凑巧的是，如果b为31，这个最大质数是2^31-1。确定这个m以后，a的选择满足：min(使a^l-1能被m整除的l)=m-1。这样的选择组合下，周期是m-1。目前大多数软件选择将a设置为a1=7^5=16807或a2=630360016。

这种随机数发生器称为PMMLCG（prime modulus multiplicative linear congruential generators）.

用C语言实现的PMMLCG源代码如下，代码由Averill M. Law提供。

文件1：lcgrand.h

float lcgrand(int stream);

void   lcgrandst(long zset, int stream);

long   lcgrandgt(int stream);

文件2: lcgrand.c

#define MODLUS 2147483647

#define MULT1          24112

#define MULT2          26143

static long zrng[] =

{             1,  

  1973272912, 281629770,   20006270,1280689831,2096730329,1933576050,

   913566091, 246780520,1363774876, 604901985,1511192140,1259851944,

   824064364, 150493284, 242708531,   75253171,1964472944,1202299975,

   233217322,1911216000, 726370533, 403498145, 993232223,1103205531,

   762430696,1922803170,1385516923,   76271663, 413682397, 726466604,

   336157058,1432650381,1120463904, 595778810, 877722890,1046574445,

     68911991,2088367019, 748545416, 622401386,2122378830, 640690903,

  1774806513,2132545692,2079249579,   78130110, 852776735,1187867272,

  1351423507,1645973084,1997049139, 922510944,2045512870, 898585771,

   243649545,1004818771, 773686062, 403188473, 372279877,1901633463,

   498067494,2087759558, 493157915, 597104727,1530940798,1814496276,

   536444882,1663153658, 855503735,   67784357,1432404475, 619691088,

   119025595, 880802310, 176192644,1116780070, 277854671,1366580350,

  1142483975,2026948561,1053920743, 786262391,1792203830,1494667770,

  1923011392,1433700034,1244184613,1147297105, 539712780,1545929719,

   190641742,1645390429, 264907697, 620389253,1502074852, 927711160,

   364849192,2049576050, 638580085, 547070247 };

float lcgrand(int stream)

{

      long zi, lowprd, hi31;

      zi       = zrng[stream];

      lowprd = (zi & 65535) * MULT1;

      hi31    = (zi >> 16) * MULT1 + (lowprd >> 16);

      zi       = ((lowprd & 65535) - MODLUS) +

                    ((hi31 & 32767) << 16) + (hi31 >> 15);

      if (zi < 0) zi += MODLUS;

      lowprd = (zi & 65535) * MULT2;

      hi31    = (zi >> 16) * MULT2 + (lowprd >> 16);

      zi       = ((lowprd & 65535) - MODLUS) +

                    ((hi31 & 32767) << 16) + (hi31 >> 15);

      if (zi < 0) zi += MODLUS;

      zrng[stream] = zi;

      return (zi >> 7 | 1) / 16777216.0;

}

void lcgrandst (long zset, int stream)

{

      zrng[stream] = zset;

}

long lcgrandgt (int stream)

{

      return zrng[stream];

}

文件3：test.c

#include

#include "lcgrand.h"

int main(){

      int i = 0;

      for(i=0;i<10;i++)

      {    

            printf("%f\n",lcgrand(0));

      }    

}

文件4：Makefile

CC=gcc

MAKER_SRC=$(wildcard *.c)

MAKER_OBJ=$(patsubst %.c,%.o,$(MAKER_SRC))

MAKER_EXE=test

.PHONY:all

all:$(MAKER_EXE)

$(MAKER_EXE):$(MAKER_OBJ)

      $(CC) -o $@ $^

%.o:%.c

      $(CC) -c -o $@ $<</div>

.PHONY:clean