浅谈随机数的生成
Part0:隨機數(shù)的性質(zhì)
隨機數(shù)一般來說符合下面這幾個性質(zhì).
(馬爾科夫性)(1.)它產(chǎn)生時后面那個數(shù)與前面的毫無關(guān)系.
(不確定性)(2.)給定樣本的一部分和隨機算法,無法推出樣本的剩余部分.
(不可再現(xiàn)性)(3.)其隨機樣本不可重現(xiàn).
另外還要說一下統(tǒng)計學偽隨機數(shù)概念.
統(tǒng)計學偽隨機性.統(tǒng)計學偽隨機性指的是在給定的隨機比特流樣本中,1的數(shù)量大致等于0的數(shù)量,同理,"10""01""00""11"四者數(shù)量大致相等.類似的標準被稱為統(tǒng)計學隨機性.滿足這類要求的數(shù)字在人類"一眼看上去”是隨機的.(摘自百度詞條)
實際上這也是在計算機中對偽隨機數(shù)優(yōu)劣的概念.
Part1:偽隨機數(shù)
偽隨機數(shù),也就是我們C++中常用的隨機數(shù).為什么說它是"偽"隨機呢?其實只要稍微詳細的了解過C++rand函數(shù)的人應該都能懂得這一點.
因為計算機本身并不能夠產(chǎn)生和隨機數(shù),只能通過產(chǎn)生一組循環(huán)節(jié)很長的數(shù)來"偽造"隨機.
C++的rand函數(shù)實際上只是根據(jù)你提供的種子計算出來的一組循環(huán)節(jié)很長的數(shù).只要兩次提供的種子是一樣的,那么rand函數(shù)提供的隨機數(shù)也是一樣的.
那么,循環(huán)節(jié)到底長到什么程度呢?
事實上,rand()函數(shù)在LINUX系統(tǒng)下的實現(xiàn)如下：
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函數(shù)的循環(huán)節(jié)應該是在(32768(2^{15}))以內(nèi).
因此,在計算機方面,目前來說,如果不借助外部幫助,是無法達到真隨機的.
Part2:隨機數(shù)的優(yōu)劣判定
在講隨機數(shù)算法之前,應該先講講隨機數(shù)優(yōu)劣的判定.畢竟只有清除了隨機數(shù)的優(yōu)劣,我們才能說如何生成優(yōu)質(zhì)隨機數(shù).
在這里我們就要用到前面說的統(tǒng)計學偽隨機性:
統(tǒng)計學偽隨機性.統(tǒng)計學偽隨機性指的是在給定的隨機比特流樣本中,1的數(shù)量大致等于0的數(shù)量,同理,"10""01""00""11"四者數(shù)量大致相等.類似的標準被稱為統(tǒng)計學隨機性.滿足這類要求的數(shù)字在人類"一眼看上去”是隨機的.(摘自百度詞條)
結(jié)合計算機隨機數(shù)的特性,我們能夠得出以下三項判斷隨機數(shù)優(yōu)劣的性質(zhì):
(1.)隨機程度,即隨機算法是否足夠復雜,隨機數(shù)之間會不會存在明顯聯(lián)系.
(2.)分布范圍,即是否存在隨機數(shù)在分布區(qū)域內(nèi)大量出現(xiàn)偏大偏小的現(xiàn)象,分布是否足夠平均.
(3.)循環(huán)長度,即是否會在大量調(diào)用時很快地出現(xiàn)循環(huán)的情況.
有了這些評判規(guī)則,我們就比較好學習優(yōu)質(zhì)隨機數(shù)的生成.
Part3:基于C++rand的優(yōu)質(zhì)隨機數(shù)的生成
我們先來講一下基于C++rand函數(shù)的隨機數(shù)生成.
1.來回擺動式
這種隨機數(shù)主要是針對退火算法之類的需要用隨機數(shù)來修正答案的.既然是修正答案,那么我們希望最好是來回擺動,一正一負.這種隨機數(shù)的特點便是通過一部分人工處理,將原本的rand函數(shù)產(chǎn)生的隨機數(shù)變成正負交替的.
static int f=3000;
static double del=0.999;// f和del是用來控制隨機數(shù)幅度不斷變小的
static int con=-1;
static int g=1; // 控制正負交替
int rand1()
{
f*=del;
g*=con;
int tmp=f*g*rand();
return tmp;
}
這種隨機數(shù)的產(chǎn)生引入了退火的思路,當然,你也可以直接使用算法中現(xiàn)成的溫度來控制.
2.平均式
這種主要是用于平衡樹treap的,特點就是在保證單個數(shù)隨機的情況下在整體上保證分布比較平均.
int p; // 希望的分布位置
int rand2()
{
int tmp=(p+rand())/2; // 通過取于分布位置的平均數(shù),是產(chǎn)生的數(shù)更加靠近期望分布
return tmp;
}
3.多次調(diào)用不重復式
當然,如果有人真的需要非常接近真隨機的數(shù),也就是多次運行程序也不會出現(xiàn)相同的情況,那就需要用到一定的外部干擾了.
首先是clock函數(shù).上文已經(jīng)說過,一個程序在不斷調(diào)用期間.每一次的運行時間都會有細小的變化.我們就可以利用好這個變化.每次調(diào)用完后都重置一次隨機數(shù)種子.
還有一個可能大家都會忽視的方法,計算機本身的誤差.眾所周知,計算機在做浮點運算時是會產(chǎn)生精度損失的,那么我們也可以利用這個特點輔助``clock調(diào)整種子(畢竟程序調(diào)用時間相同其實可能性也不小,畢竟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()); // 重置隨機數(shù)種子
return rand();
}
經(jīng)過大量實驗,筆者發(fā)現(xiàn)該函數(shù)前三個數(shù)出現(xiàn)重復幾率相對會比較大(7~9%)建議從第四個開始使用.
上面的代碼并沒用用精度損失來隨機化,因為同一個式子的進度損失值太小,以至于幾乎不會發(fā)生什么改變,所以并沒有使用.
優(yōu)劣度分析
首先,隨機程度方面,雖然之前看過rand()函數(shù)代碼,可能清楚數(shù)字之間的關(guān)聯(lián).但在實際運用中,這個數(shù)字之間的關(guān)聯(lián)還是基本可以忽略的.所以在隨機程度方面,rand()函數(shù)還是能夠勉強通過的.
在平均分布方面,單看代碼可能感覺不出來.
那么,筆者就做一個測試:
int data[10007];
int main()
{
for(int i=1;i<=1000000;++i)
{
int tmp=rand()%100000; // 生成一個100000以內(nèi)的隨機數(shù)
++data[tmp/10]; // 統(tǒng)計出現(xiàn)次數(shù)
}
for(int i=1;i<=1000;++i)
printf("%d
",data[i]);
}
最后結(jié)果：
從中我們可以看到,這個分布還是非常平均的.
循環(huán)長度...
這個主要就是rand()函數(shù)的硬傷了,(32768)這個長度真的挺不夠用的.在需要大量調(diào)用rand()函數(shù)的算法中(比如退火),基本都會把rand()卡出循環(huán).
那有沒有既優(yōu)質(zhì)又循環(huán)節(jié)長的算法呢?
Part4:梅森旋轉(zhuǎn)算法(MT19937)
梅森旋轉(zhuǎn)算法是目前產(chǎn)生優(yōu)質(zhì)偽隨機數(shù)的普遍算法,在C++11的標準中,也增加了MT19937的實現(xiàn),在random庫中.
其實,這個算法是由松本真和和西村拓士在1997年開發(fā)的一種算法,和梅森沒有多大關(guān)系.它之所以有這個名字,是因為它有長達(2^{19937}-1)的循環(huán)節(jié),這是一個梅森素數(shù).況且,這種算法還能在如此長的循環(huán)節(jié)下產(chǎn)生均勻的隨機數(shù).
那么,MT19937的原理究竟是什么呢?
這個旋轉(zhuǎn)算法實際上是對一個(19937)位的二進制序列作變換.我們知道,一個長度為(n)的二進制序列,它的排列長度最長為(2^n).但是,有時候會因為某些操作不當,導致循環(huán)節(jié)小于(2^n).而如何將這個序列的排列恰好達到(2^n)個,就是這個旋轉(zhuǎn)算法的精髓.
如果反饋函數(shù)的本身(+1)是一個既約多項式,那么它的循環(huán)節(jié)達到最長,即(2^n-1).
我們拿一個4位寄存器模擬一下:
我們這里使用的反饋函數(shù)是(y=x^4+x^2+x+1)(這個不是既約多項式,只是拿來好理解)
這個式子中(x^4),(x^2),(x)的意思就是我們每次對這個二進制序列的從后往前數(shù)第4位和第2位做異或運算,然后(x)的意思是我們再拿結(jié)果和最后一位做異或運算.把最后的結(jié)果放到序列的開頭,整個序列后移一位,最后一位舍棄.
1.初始數(shù)組({1,0,0,0}).
2.將它的第四位和第二位取出來做異或運算.
3.把剛剛的運算結(jié)果和最后一位再做一次運算
4.把最后的運算結(jié)果放到第一位,序列后移.最后一位被拋棄.
這就是一次運算,然后這個算法就是不斷循環(huán)這幾步,從而不斷偽隨機改變這個序列.
因為它所使用的反饋函數(shù)(y=x^4+x+1)是既約多項式,所以最后的循環(huán)節(jié)為(2^4-1=15),運算結(jié)果如下：
[egin{array}
{|c|c|c|c|}\
a_3&a_2&a_1&a_0\
1&0&0&0\
1&1&0&0\
1&1&1&0\
1&1&1&1\
0&1&1&1\
1&0&1&1\
0&1&0&1\
1&0&1&0\
1&1&0&1\
0&1&1&0\
0&0&1&1\
1&0&0&1\
0&1&0&0\
0&0&1&0\
0&0&0&1\
-&-&-&-\
1&0&0&0
end{array}
]
大家可以看到,這個運算結(jié)果包含了(1,2,...,2^4-1)中的所有整數(shù),并且沒有循環(huán),同時擁有很好的隨機性.
Part5:MT19937的偽代碼及C++實現(xiàn)
初始化隨機種子:
(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})
梅森旋轉(zhuǎn):
( 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)
生成隨機數(shù):
( 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++實現(xiàn):
int index;
int MT[624];
// 設(shè)置隨機數(shù)種子
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;
}
}
// 梅森旋轉(zhuǎn)
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;
}
}
// 生成隨機數(shù)
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;
}
本文完
總結(jié)
- 上一篇： 红烧肉家庭做法（一定要收藏的红烧肉，简单
- 下一篇： 七：策略模式（不同等级会员打折算法）