（
newsuppy,
轉(zhuǎn)載請(qǐng)注明出處）
當(dāng)我們需要在某個(gè)任務(wù)中使用隨機(jī)數(shù)，通常我們習(xí)慣于使用標(biāo)準(zhǔn)庫(kù)的rand函數(shù)。像這樣：srand(time(0)); // 時(shí)間種子
??????? rand() % MAX_RAND
；
標(biāo)準(zhǔn)庫(kù)的rand函數(shù)使用線性同余算法，是生成速度相當(dāng)快的一種隨機(jī)數(shù)生成算法。在多數(shù)情況下也確實(shí)能滿足我們的要求，但是對(duì)于一些特殊目的應(yīng)用這個(gè)算法生成的隨機(jī)數(shù)是不行的，比如某些加密算法，蒙特卡羅積分等（在.NET中創(chuàng)建隨機(jī)密碼的加密安全隨機(jī)數(shù)就不能使用Random類(lèi)的線性同余隨機(jī)數(shù)，而要使用System.Security.Cryptography命名空間中的相關(guān)隨機(jī)數(shù)生成類(lèi)）。
這個(gè)線性同余算法的實(shí)現(xiàn)可以在很多書(shū)籍中找到。下面我給出一個(gè)The C Programming Langurage 2nd中的一個(gè)實(shí)現(xiàn)，這也是普遍使用的標(biāo)準(zhǔn)庫(kù)隨機(jī)數(shù)算法的實(shí)現(xiàn)：
|
??
unsigned
long
int next =1;
?
??
/* rand:? return pseudo-random integer on 0..32767 */
?
?
int rand(void)
??
{
?????? next = next *1103515245+12345;
??????
return
(
unsigned
int
)(next/65536)%32768;
??
}
?
??
/* srand:? set seed for rand() */
??
void srand(unsignedint seed)
??
{
?????? next = seed;
??
}
?
|
這個(gè)實(shí)現(xiàn)的問(wèn)題在于rand函數(shù)return行中的那個(gè)32768，在標(biāo)準(zhǔn)庫(kù)中這個(gè)數(shù)字定義為RAND_MAX宏，在VisualC++和Mingw32編譯器的stdlib.h頭文件（或者cstdlib)中你都可以發(fā)現(xiàn)RAND_MAX的值為32768。也就是說(shuō)這個(gè)算法的隨機(jī)數(shù)分布在0--RAND_MAX中，而在一般編譯器中就是0--32768。假設(shè)你的算法需要的是300000多個(gè)的隨機(jī)數(shù)，那么使用rand函數(shù)會(huì)產(chǎn)生重負(fù)次數(shù)近30次的隨機(jī)數(shù)!
所以在這里我將簡(jiǎn)單介紹一下如何使用Blitz++庫(kù)中的隨機(jī)數(shù)生成類(lèi)。不過(guò)在這里我不能夠證明Blitz++隨機(jī)數(shù)算法相對(duì)于標(biāo)準(zhǔn)庫(kù)有什么優(yōu)越性。Blitz++的源代碼是開(kāi)放的，你可以完全了解它的隨機(jī)數(shù)算法的實(shí)現(xiàn)。
在Blitz++中隨機(jī)數(shù)類(lèi)組織在ranlib namespace中，使用方法也非常簡(jiǎn)單，seed成員函數(shù)種入種子后，再用random成員函數(shù)就可以了。Blitz++中包括了產(chǎn)生各種概率分布情況的隨機(jī)數(shù)，列舉如下：均勻分布(Uniform),正態(tài)分布(Normal)，指數(shù)分布(Exponential),Beta分布，Gamma分布，Χ方分布，F分布，離散均勻分布。具體地可以參考Blitz++文檔的第九章。本文將會(huì)演示正態(tài)分布的隨機(jī)數(shù)類(lèi)。
NormalUnit<>()??
標(biāo)準(zhǔn)正態(tài)分布，
μ
= 0,
σ
= 1;
?
Normal<>(T mean, T standardDeviation)
正態(tài)分布，N(
μ
,
σ
^2
)
，其中mean就是
μ
,standardDeviation
就是
σ
。
Blitz++
中隨機(jī)數(shù)類(lèi)的seed是共享的，調(diào)用一個(gè)類(lèi)的seed后，其他類(lèi)也同樣種入了相同的種子。對(duì)于種子的來(lái)源，Blitz++文檔中有這樣一個(gè)注意事項(xiàng):
Note: you may be tempted to seed the random number generator from a static initializer. Don't do it! Due to an oddity of C++, there is no guarantee on the order of static initialization when templates are involved. Hence, you may seed the RNG before its constructor is invoked, in which case your program will crash. If you don't know what a static initializer is, don't worry -- you're safe!
（幸運(yùn)的是我也不太清楚static initializer具體指什么,:-） )
1
，先來(lái)看一個(gè)標(biāo)準(zhǔn)正態(tài)分布的例子，根據(jù)3
σ
法則(正態(tài)分布的隨機(jī)數(shù)幾乎全部落在(
μ
-3
σ
,
μ
+3
σ
)
中),這些隨機(jī)數(shù)應(yīng)該大部分落在(-3, 3)之間。
下面的程序產(chǎn)生10000個(gè)正態(tài)分布隨機(jī)數(shù)，然后導(dǎo)入Matlab生成圖形，橫坐標(biāo)是隨機(jī)數(shù)的序數(shù)，縱坐標(biāo)是隨機(jī)數(shù)值。很顯然，大部分點(diǎn)在3
σ區(qū)間內(nèi)。在區(qū)間
(
μ
-
σ
,
μ
+
σ)中數(shù)據(jù)點(diǎn)最密。
|
#include
<random/uniform.h>
#include
<random/normal.h>
#include
<iostream>
#include
<fstream>
#include
<ctime>
#include
<cstdlib>
using
namespace std;
using
namespace ranlib;
?
int
main()
{
???
const size_t MAXGEN =10000;
??? NormalUnit<float> rnd;
?
??? rnd.seed(time(0));
?
??? ofstream file("c:\\file.txt");
?
???
// generator Normal distribution random number
???
for
(size_t i=0; i<MAXGEN;++i)
???
{
??????? file << rnd.random()<<" ";
???
}
}
|
準(zhǔn)庫(kù)rand()函數(shù)的缺陷以及Blitz++隨機(jī)數(shù)生成的簡(jiǎn)介.files/image001.jpg)
2
，
下面是一個(gè)服從于
μ= 10，σ= 2的正態(tài)分布，根據(jù)3σ法則，大部分點(diǎn)應(yīng)該落在(4,16)之間。
|
#include
<random/uniform.h>
#include
<random/normal.h>
#include
<iostream>
#include
<fstream>
#include
<ctime>
#include
<cstdlib>
using
namespace std;
using
namespace ranlib;
?
int
main()
{
???
const size_t MAXGEN =1000000;
???
const
int mu =10;
???
const
int sigma =2;
???
Normal
<
float
> rnd(mu,sigma);
?
??? rnd.seed(time(0));
?
??? ofstream file("c:\\file.txt");
?
???
// generator Normal distribution random number
???
for
(size_t i=0; i<MAXGEN;++i)
???
{
??
?????file << rnd.random()<<" ";
???
}
}
|
|
準(zhǔn)庫(kù)rand()函數(shù)的缺陷以及Blitz++隨機(jī)數(shù)生成的簡(jiǎn)介.files/image002.jpg)
|
?
3，生成前述正態(tài)分布的鐘型曲線（PDF）。
這里產(chǎn)生1M的float隨機(jī)數(shù)，然后乘以10轉(zhuǎn)型為整數(shù)，統(tǒng)計(jì)在區(qū)間0-170之間隨機(jī)數(shù)出現(xiàn)的概率，這樣再導(dǎo)入Matlab生成圖形就是一個(gè)近似的正態(tài)分布的鐘型曲線了。
?
?
|
#include
<random/uniform.h>
#include
<random/normal.h>
#include
<iostream>
#include
<fstream>
#include
<ctime>
#include
<cstdlib>
using
namespace std;
using
namespace ranlib;
?
int
main()
{
???
const size_t MAXGEN =1000000;
???
const
int mu =10;
???
const
int sigma =2;
???
Normal
<
float
> rnd(mu,sigma);
?
??? rnd.seed(time(0));
?
??? ofstream file("c:\\file.txt");
?
???
float
*rndArray =newfloat[MAXGEN];
???
// generator Normal distribution random number
???
for
(size_t i=0; i<MAXGEN;++i)
???
{
??????? rndArray[i]= rnd.random();
???
}
?
???
int
*rndConvertIntArray =newint[MAXGEN];
???
int multiple =10;
???
// convert float random number to integer
???
for
(size_t i=0; i<MAXGEN;++i)
???
{
??????? rndConvertIntArray[i]=int(rndArray[i]* multiple);
???
}
?
???
const size_t PDFSIZE =(mu +3* sigma)* multiple;
???
const size_t redundancy =10;
???
int
*pdf =newint[PDFSIZE+redundancy];
???
for
(size_t i=0; i<PDFSIZE+redundancy;++i)
???
{
??????? pdf[i]=0;
???
}
?
???
// generator PDF(Probability Distribution Function), Normal distribution is a "bell" shape
???
for
(size_t i=0; i<MAXGEN;++i)
???
{
??
?????
if
(rndConvertIntArray[i]<= PDFSIZE+redundancy-1)
???????
{
???????????
++pdf[rndConvertIntArray[i]];
???????
}
???
}
?
???
for
(size_t i=0; i<PDFSIZE+redundancy;++i)
???
{
??????? file << pdf[i]? <<" ";
???
}
?
???
delete
[] rndArray;
???
delete
[] rndConvertIntArray;
???
delete
[] pdf;
}
|
|
?
|
鐘型曲線，當(dāng)然不是特別光滑
(
J
),
大部分隨機(jī)數(shù)都出現(xiàn)在（
4
，
16
）區(qū)間內(nèi)。
?
總結(jié)，
Blitz++
的隨機(jī)數(shù)生成類(lèi)應(yīng)該說(shuō)是值得信任的。
[
參考文獻(xiàn)
]
1，?
概率論與數(shù)理統(tǒng)計(jì)（
3rd
），
浙江大學(xué)
盛驟
謝式千
潘承毅
編
高等教育出版社
2，?
C
數(shù)值算法（
2nd
），
[
美
] William H. Press, Saul A. Teukolsky, William T. Vetterling,
???????????????????
?????Brian P. Flannery
著
?????????????????????
傅祖蕓
趙梅娜
丁巖石
等譯，傅祖蕓
審校
3
，
Matlab 6.5
的文檔
??? The MathWorks, Inc.