浅谈随机数的生成

Posted 2020-11-21 anverking

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),运算结果如下：

[ egin{array} {|c|c|c|c|}a_3&a_2&a_1&a_01&0&0&01&1&0&01&1&1&01&1&1&1&1&1&11&0&1&1&1&0&11&0&1&01&1&0&1&1&1&0&0&1&11&0&0&1&1&0&0&0&1&0&0&0&1-&-&-&-1&0&0&0 end{array} ]

大家可以看到,这个运算结果包含了(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 sramd(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;
 } 

本文完