Part0:随机数的性质
随机数一般来说符合下面这几个性质.
(马尔科夫性)(1.)它产生时后面那个数与前面的毫无关系.
(不确定性)(2.)给定样本的一部分和随机算法,无法推出样本的剩余部分.
(不可再现性)(3.)其随机样本不可重现.
另外还要说一下统计学伪随机数概念.
统计学伪随机性.统计学伪随机性指的是在给定的随机比特流样本中,1的数量大致等于0的数量,同理,"10""01""00""11"四者数量大致相等.类似的标准被称为统计学随机性.满足这类要求的数字在人类"一眼看上去”是随机的.(摘自百度词条)
实际上这也是在计算机中对伪随机数优劣的概念.
Part1:伪随机数
伪随机数,也就是我们C++中常用的随机数.为什么说它是"伪"随机呢?其实只要稍微详细的了解过C++rand函数的人应该都能懂得这一点.
因为计算机本身并不能够产生和随机数,只能通过产生一组循环节很长的数来"伪造"随机.
C++的rand函数实际上只是根据你提供的种子计算出来的一组循环节很长的数.只要两次提供的种子是一样的,那么rand函数提供的随机数也是一样的.
那么,循环节到底长到什么程度呢?
事实上,rand()函数在LINUX系统下的实现如下：
static unsigned long next=1;
// RAND_MAX assumed to be 32767
int rand()
{
next=next*1103515245+12345;
return ((unsigned)(next/65536)%32768);
}
void srand(unsigned seed)
{
next=seed;
}
通过这个算法我们可以推知,rand函数的循环节应该是在(32768(2^{15}))以内.
因此,在计算机方面,目前来说,如果不借助外部帮助,是无法达到真随机的.
Part2:随机数的优劣判定
在讲随机数算法之前,应该先讲讲随机数优劣的判定.毕竟只有清除了随机数的优劣,我们才能说如何生成优质随机数.
在这里我们就要用到前面说的统计学伪随机性:
统计学伪随机性.统计学伪随机性指的是在给定的随机比特流样本中,1的数量大致等于0的数量,同理,"10""01""00""11"四者数量大致相等.类似的标准被称为统计学随机性.满足这类要求的数字在人类"一眼看上去”是随机的.(摘自百度词条)
结合计算机随机数的特性,我们能够得出以下三项判断随机数优劣的性质:
(1.)随机程度,即随机算法是否足够复杂,随机数之间会不会存在明显联系.
(2.)分布范围,即是否存在随机数在分布区域内大量出现偏大偏小的现象,分布是否足够平均.
(3.)循环长度,即是否会在大量调用时很快地出现循环的情况.
有了这些评判规则,我们就比较好学习优质随机数的生成.
Part3:基于C++rand的优质随机数的生成
我们先来讲一下基于C++rand函数的随机数生成.
1.来回摆动式
这种随机数主要是针对退火算法之类的需要用随机数来修正答案的.既然是修正答案,那么我们希望最好是来回摆动,一正一负.这种随机数的特点便是通过一部分人工处理,将原本的rand函数产生的随机数变成正负交替的.
static int f=3000;
static double del=0.999;// f和del是用来控制随机数幅度不断变小的
static int con=-1;
static int g=1; // 控制正负交替
int rand1()
{
f*=del;
g*=con;
int tmp=f*g*rand();
return tmp;
}
这种随机数的产生引入了退火的思路,当然,你也可以直接使用算法中现成的温度来控制.
2.平均式
这种主要是用于平衡树treap的,特点就是在保证单个数随机的情况下在整体上保证分布比较平均.
int p; // 希望的分布位置
int rand2()
{
int tmp=(p+rand())/2; // 通过取于分布位置的平均数,是产生的数更加靠近期望分布
return tmp;
}
3.多次调用不重复式
当然,如果有人真的需要非常接近真随机的数,也就是多次运行程序也不会出现相同的情况,那就需要用到一定的外部干扰了.
首先是clock函数.上文已经说过,一个程序在不断调用期间.每一次的运行时间都会有细小的变化.我们就可以利用好这个变化.每次调用完后都重置一次随机数种子.
还有一个可能大家都会忽视的方法,计算机本身的误差.众所周知,计算机在做浮点运算时是会产生精度损失的,那么我们也可以利用这个特点辅助``clock调整种子(毕竟程序调用时间相同其实可能性也不小,毕竟clock```只精确到( ext{ms})).
int count;
int rand3()
{
++count;
int t=clock()+1; // 使用当前时间
for(int i=1;i<12121307;++i) // 降速
t+=rand();
t+=clock(); // 降速后扩大时间变化
t*=-1234;
srand(t*count+rand()); // 重置随机数种子
return rand();
}
经过大量实验,笔者发现该函数前三个数出现重复几率相对会比较大(7~9%)建议从第四个开始使用.
上面的代码并没用用精度损失来随机化,因为同一个式子的进度损失值太小,以至于几乎不会发生什么改变,所以并没有使用.
优劣度分析
首先,随机程度方面,虽然之前看过rand()函数代码,可能清楚数字之间的关联.但在实际运用中,这个数字之间的关联还是基本可以忽略的.所以在随机程度方面,rand()函数还是能够勉强通过的.
在平均分布方面,单看代码可能感觉不出来.
那么,笔者就做一个测试:
int data[10007];
int main()
{
for(int i=1;i<=1000000;++i)
{
int tmp=rand()%100000; // 生成一个100000以内的随机数
++data[tmp/10]; // 统计出现次数
}
for(int i=1;i<=1000;++i)
printf("%d
",data[i]);
}
最后结果：
从中我们可以看到,这个分布还是非常平均的.
循环长度...
这个主要就是rand()函数的硬伤了,(32768)这个长度真的挺不够用的.在需要大量调用rand()函数的算法中(比如退火),基本都会把rand()卡出循环.
那有没有既优质又循环节长的算法呢?
Part4:梅森旋转算法(MT19937)
梅森旋转算法是目前产生优质伪随机数的普遍算法,在C++11的标准中,也增加了MT19937的实现,在random库中.
其实,这个算法是由松本真和和西村拓士在1997年开发的一种算法,和梅森没有多大关系.它之所以有这个名字,是因为它有长达(2^{19937}-1)的循环节,这是一个梅森素数.况且,这种算法还能在如此长的循环节下产生均匀的随机数.
那么,MT19937的原理究竟是什么呢?
这个旋转算法实际上是对一个(19937)位的二进制序列作变换.我们知道,一个长度为(n)的二进制序列,它的排列长度最长为(2^n).但是,有时候会因为某些操作不当,导致循环节小于(2^n).而如何将这个序列的排列恰好达到(2^n)个,就是这个旋转算法的精髓.
如果反馈函数的本身(+1)是一个既约多项式,那么它的循环节达到最长,即(2^n-1).
我们拿一个4位寄存器模拟一下:
我们这里使用的反馈函数是(y=x^4+x^2+x+1)(这个不是既约多项式,只是拿来好理解)
这个式子中(x^4),(x^2),(x)的意思就是我们每次对这个二进制序列的从后往前数第4位和第2位做异或运算,然后(x)的意思是我们再拿结果和最后一位做异或运算.把最后的结果放到序列的开头,整个序列后移一位,最后一位舍弃.
1.初始数组({1,0,0,0}).
2.将它的第四位和第二位取出来做异或运算.
3.把刚刚的运算结果和最后一位再做一次运算
4.把最后的运算结果放到第一位,序列后移.最后一位被抛弃.
这就是一次运算,然后这个算法就是不断循环这几步,从而不断伪随机改变这个序列.
因为它所使用的反馈函数(y=x^4+x+1)是既约多项式,所以最后的循环节为(2^4-1=15),运算结果如下：
大家可以看到,这个运算结果包含了(1,2,...,2^4-1)中的所有整数,并且没有循环,同时拥有很好的随机性.
Part5:MT19937的伪代码及C++实现
初始化随机种子:
(indexleftarrow 0)
(MTleftarrow new array with size 624//624 imes 32-31=19937)
(// ext{Above are global variables})
( ext{MT19937_SRAND}(seed):)
(indexleftarrow 0)
(MT[0]leftarrow seed)
(mathbf{for} ileftarrow 1 mathbf{to} 623:)
(quad tleftarrow 1812433253cdot(MT[i-1]oplus(MT[i-1]gg 30))+i//oplus ext{ is the xor operation}, gg ext{ is the right-shift operation})
(quad MT[i]leftarrow t& ext{0xffffffff}// ext{get the last 32 bits})
(//& ext{ is the bit-and operation}, ext{0x means that the number next is a hex number})
梅森旋转:
( ext{MT19937_GENERATE}():)
(mathbf{for} ileftarrow 0 mathbf{to} 623:)
(quad yleftarrow (MT[i]& ext{0x80000000})+(MT[(i+1)mod 624]& ext{0x7fffffff}))
(quad MT[i]leftarrow MT[(i+397)mod 624]oplus(ygg 1))
(quadmathbf{if} y&1=1:)
(quadquad MT[i]leftarrow MT[i]oplus 2567483615)
生成随机数:
( ext{MT19937_RAND}():)
(mathbf{if} index=0:)
(quad ext{MT19937_GENERATE}())
(yleftarrow MT[index])
(yleftarrow yoplus (ygg 11))
(yleftarrow yoplus ((yll 7)&2636928640))
(yleftarrow yoplus ((yll 15)& 4022730752))
(yleftarrow yoplus (ygg 18))
(indexleftarrow (index+1)mod 624)
(mathbf{return} y)
C++实现:
int index;
int MT[624];
// 设置随机数种子
inline void srand(int seed)
{
index=0;
MT[0]=seed;
for(int i=1;i<=623;++i)
{
int t=1812433253*(MT[i-1]^(MT[i-1]>>30))+i;
MT[i]=t&0xffffffff;
}
}
// 梅森旋转
inline void generate()
{
for(int i=0;i<=623;++i)
{
int y=(MT[i]&0x80000000)+(MT[(i+1)%624]&0x7fffffff);
MT[i]=MT[(i+397)%624]^(y>>1);
if(y&1)
MT[i]^=2567483615;
}
}
// 生成随机数
inline int rand()
{
if(index==0)
generate();
int y=MT[index];
y=y^(y>>1);
y=y^((y<<7)&2636928640);
y=y^((y<<15)&4022730752);
y=y^(y>>18);
index=(index+1)%624;
return y;
}