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