【基础知识】vc++中的随机数生成算法

前段时间，看网络上流传的靳志辉写的《LDA数学八卦》，里面提到了蒙特卡洛模拟算法，用来抽样模拟特定概率分布的随机样本序列。文中简单提到了，如何生成均匀分布（同余发生器）的随机样本序列以及根据均匀分布能够模拟生成其他主要分布（例如：正态分布）的样本序列。好奇地查了一下现在的均匀分布的模拟方法。

随机数可以用物理模拟，我猜想是用某些满足随机分布条件的物理值来实时地模拟。不过程序中经常用的是递推模拟，即根据当前的随机数的值，计算生成下一个随机数的值，逐渐递推下去，形成的数列满足特定分布的统计性质，从而作为特定分布的随机样本的模拟。但是既然这些数是计算出来的，即是说他们是确定的，所以称之为伪随机数序列。在这个过程中唯一不确定的就是初始的随机数的值（后续的值都是根据他计算出来的），所记载模拟过程中这个值很重要，通常选择计算机的当前时间作为初始值。

生成伪随机数序列最常用的方法是线性取余法，又叫做线性同余发生器。线性同余发生器的递推公式是

X(n+1) = (a * X(n) + c) % m

其中X(n)是第n个随机数，从中计算出第n+1个随机数。当c=0的时候，就退化到了乘法同余发生器。

看了下vc++ 2008中的代码，用的是这种方法来生成随机数，代码如下：

/***
*void srand(seed) - seed the random number generator
*
*Purpose:
*       Seeds the random number generator with the int given.  Adapted from the
*       BASIC random number generator.
*
*Entry:
*       unsigned seed - seed to seed rand # generator with
*
*Exit:
*       None.
*
*Exceptions:
*
*******************************************************************************/

void __cdecl srand (
unsigned int seed
)
{

// 将种子保存到 _getptd()->_holdrand 中
_getptd()->_holdrand = (unsigned long)seed;

}

/***
*int rand() - returns a random number
*
*Purpose:
*       returns a pseudo-random number 0 through 32767.
*
*Entry:
*       None.
*
*Exit:
*       Returns a pseudo-random number 0 through 32767.
*
*Exceptions:
*
*******************************************************************************/

int __cdecl rand (
void
)
{

_ptiddata ptd = _getptd();

// 根据种子和线性同余发生器，产生下一个随机数
return( ((ptd->_holdrand = ptd->_holdrand * 214013L
+ 2531011L) >> 16) & 0x7fff );

}

在上面的代码中，a =214013，c = 2531011，这两个数是根据什么设定的，我还不太清楚。m = 32767 即 2^16，上面的右移16位（'...>>16'）相当于除32767，接下来取低16位的值（'... & 0x7fff'），相当于模后取余，注意 0x7fff 的值低16位都是1，等于32767。这样，随机数的范围就在[0, 2^16]之间。其实因果关系应该反过来，是因为写程序容易、运行效率高，采用了2^16作为m，才导致随机数的范围在[0, 32767]之间。

需要注意的是，这种方法有很多局限性，统计上有很多缺陷，因此在真正做统计模拟的时候，最好自己设计随机数发生器。

从网上找到一个参考文献，介绍的不错，《纯线性同余随机数生成器》，把原文也copy到下面了：

纯线性同余随机数生成器

线性同余随机数生成器介绍：

古老的LCG(linear congruential generator)代表了最好最朴素的伪随机数产生器算法。主要原因是容易理解，容易实现，而且速度快。

LCG 算法数学上基于公式：

X(n+1) = (a * X(n) + c) % m

其中，各系数为：

模m, m > 0

系数a, 0 < a < m

增量c, 0 <= c < m

原始值(种子) 0 <= X(0) < m

其中参数c, m, a比较敏感，或者说直接影响了伪随机数产生的质量。

一般而言，高LCG的m是2的指数次幂(一般2^32或者2^64)，因为这样取模操作截断最右的32或64位就可以了。多数编译器的库中使用了该理论实现其伪随机数发生器rand()。
下面是部分编译器使用的各个参数值：

Source	m	a	c	rand() / Random(L)的种子位
Numerical Recipes	2^32	1664525	1013904223
Borland C/C++	2^32	22695477	1	位30..16 in rand(), 30..0 in lrand()
glibc (used by GCC)	2^32	1103515245	12345	位30..0
ANSI C: Watcom, Digital Mars, CodeWarrior, IBM VisualAge C/C++	2^32	1103515245	12345	位30..16
Borland Delphi, Virtual Pascal	2^32	134775813	1	位63..32 of (seed * L)
Microsoft Visual/Quick C/C++	2^32	214013	2531011	位30..16
Apple CarbonLib	2^31-1	16807	0	见Park–Miller随机数发生器

LCG不能用于随机数要求高的场合，例如不能用于Monte Carlo模拟，不能用于加密应用。

LCG有一些严重的缺陷，例如如果LCG用做N维空间的点坐标，这些点最多位于m1/n超平面上(Marsaglia定理)，这是由于产生的相继X(n)值的关联所致。

另外一个问题就是如果m设置为2的指数，产生的低位序列周期远远小于整体。

一般而言，输出序列的基数b中最低n位，bk = m (k是某个整数)，最大周期bn.

有些场合LCG有很好的应用，例如内存很紧张的嵌入式中，电子游戏控制台用的小整数，使用高位可以胜任。

（2） C语言中伪随机数生成方法：rand(),srand(time(null))的解析

C语言中伪随机数生成算法实际上是采用了"线性同余法"。具体的计算如下：

Xi = (Xi-1 * A + C ) mod M

其中A,C,M都是常数（一般会取质数）。当C=0时，叫做乘同余法。引出一个概念叫seed，它会被作为X0被代入上式中，然后每次调用rand()函数都会用上一次产生的随机值来生成新的随机值。可以看出实际上用rand()函数生成的是一个递推的序列，一切值都来源于最初的 seed。所以当初始的seed取一样的时候，得到的序列都相同。

C语言里面有RAND_MAX这样一个宏，定义了rand()所能得到的随机值的范围。在C里可以看到RAND_MAX被定义成0x7fff，也就是32767。rand()函数里递推式中M的值就是32767。

线性同余法生成的是伪随机数，粗略符合均匀分布。根据中心极限定理，任何分布的噪声，通过反复相加，就可以成为高斯噪声。

函数原型

1、C++标准函数库提供一随机数生成器rand，返回0－RAND_MAX之间均匀分布的伪随机整数。 RAND_MAX必须至少为32767。rand()函数不接受参数，默认以1为种子（即起始值）。随机数生成器总是以相同的种子开始，所以形成的伪随机数列也相同，失去了随机意义。（但这样便于程序调试）

2、C++中另一函数srand（），可以指定不同的数（无符号整数变元）为种子。但是如果种子相同，伪随机数列也相同。一个办法是让用户输入种子，但是仍然不理想。

3、 比较理想的是用变化的数，比如时间来作为随机数生成器的种子。time的值每时每刻都不同。所以种子不同，所以，产生的随机数也不同。

// C++随机函数（VC program）

#include <stdio.h>

#include <iostream>

#include <time.h>

using namespace std;

#define MAX 100

int main(int argc, char* argv[])

{

srand( (unsigned)time( NULL ) );//srand()函数产生一个以当前时间开始的随机种子.应该放在for等循环语句前面 不然要很长时间等待

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

cout<<rand()%MAX<<endl;//MAX为最大值，其随机域为0~MAX-1

return 0;

}

rand()不需要参数，它会返回一个从0到最大随机数的任意整数，最大随机数的大小通常是固定的一个大整数。这样，如果你要产生0~10的10个整数，可以表达为：

int N = rand() % 11;

这样，N的值就是一个0~10的随机数，如果要产生1~10，则是这样：

int N = 1 + rand() % 11;

总结来说，可以表示为：

a + rand() % n

其中的a是起始值，n是整数的范围。

a + rand() % (b-a+1) 就表示ａ～ｂ之间的一个随机数

若要0~1的小数，则可以先取得0~10的整数，然后均除以10即可得到随机到十分位的10个随机小数，若要得到随机到百分位的随机小数，则需要先得到0~100的10个整数，然后均除以100，其它情况依此类推。

通常rand()产生的随机数在每次运行的时候都是与上一次相同的，这是有意这样设计的，是为了便于程序的调试。若要产生每次不同的随机数，可以使用srand( seed )函数进行随机化，随着seed的不同，就能够产生不同的随机数。

如大家所说，还可以包含time.h头文件，然后使用srand(time(0))来使用当前时间使随机数发生器随机化，这样就可以保证每两次运行时可以得到不同的随机数序列(只要两次运行的间隔超过1秒)。