（
newsuppy,
轉(zhuǎn)載請注明出處）
當(dāng)我們需要在某個任務(wù)中使用隨機(jī)數(shù)，通常我們習(xí)慣于使用標(biāo)準(zhǔn)庫的rand函數(shù)。像這樣：srand(time(0)); // 時間種子
??????? rand() % MAX_RAND
；
標(biāo)準(zhǔn)庫的rand函數(shù)使用線性同余算法，是生成速度相當(dāng)快的一種隨機(jī)數(shù)生成算法。在多數(shù)情況下也確實能滿足我們的要求，但是對于一些特殊目的應(yīng)用這個算法生成的隨機(jī)數(shù)是不行的，比如某些加密算法，蒙特卡羅積分等（在.NET中創(chuàng)建隨機(jī)密碼的加密安全隨機(jī)數(shù)就不能使用Random類的線性同余隨機(jī)數(shù)，而要使用System.Security.Cryptography命名空間中的相關(guān)隨機(jī)數(shù)生成類）。
這個線性同余算法的實現(xiàn)可以在很多書籍中找到。下面我給出一個The C Programming Langurage 2nd中的一個實現(xiàn)，這也是普遍使用的標(biāo)準(zhǔn)庫隨機(jī)數(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;
??
}
?
|
這個實現(xiàn)的問題在于rand函數(shù)return行中的那個32768，在標(biāo)準(zhǔn)庫中這個數(shù)字定義為RAND_MAX宏，在VisualC++和Mingw32編譯器的stdlib.h頭文件（或者cstdlib)中你都可以發(fā)現(xiàn)RAND_MAX的值為32768。也就是說這個算法的隨機(jī)數(shù)分布在0--RAND_MAX中，而在一般編譯器中就是0--32768。假設(shè)你的算法需要的是300000多個的隨機(jī)數(shù)，那么使用rand函數(shù)會產(chǎn)生重負(fù)次數(shù)近30次的隨機(jī)數(shù)!
所以在這里我將簡單介紹一下如何使用Blitz++庫中的隨機(jī)數(shù)生成類。不過在這里我不能夠證明Blitz++隨機(jī)數(shù)算法相對于標(biāo)準(zhǔn)庫有什么優(yōu)越性。Blitz++的源代碼是開放的，你可以完全了解它的隨機(jī)數(shù)算法的實現(xiàn)。
在Blitz++中隨機(jī)數(shù)類組織在ranlib namespace中，使用方法也非常簡單，seed成員函數(shù)種入種子后，再用random成員函數(shù)就可以了。Blitz++中包括了產(chǎn)生各種概率分布情況的隨機(jī)數(shù)，列舉如下：均勻分布(Uniform),正態(tài)分布(Normal)，指數(shù)分布(Exponential),Beta分布，Gamma分布，Χ方分布，F分布，離散均勻分布。具體地可以參考Blitz++文檔的第九章。本文將會演示正態(tài)分布的隨機(jī)數(shù)類。
NormalUnit<>()??
標(biāo)準(zhǔn)正態(tài)分布，
μ
= 0,
σ
= 1;
?
Normal<>(T mean, T standardDeviation)
正態(tài)分布，N(
μ
,
σ
^2
)
，其中mean就是
μ
,standardDeviation
就是
σ
。
Blitz++
中隨機(jī)數(shù)類的seed是共享的，調(diào)用一個類的seed后，其他類也同樣種入了相同的種子。對于種子的來源，Blitz++文檔中有這樣一個注意事項:
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
，先來看一個標(biāo)準(zhǔn)正態(tài)分布的例子，根據(jù)3
σ
法則(正態(tài)分布的隨機(jī)數(shù)幾乎全部落在(
μ
-3
σ
,
μ
+3
σ
)
中),這些隨機(jī)數(shù)應(yīng)該大部分落在(-3, 3)之間。
下面的程序產(chǎn)生10000個正態(tài)分布隨機(jī)數(shù)，然后導(dǎo)入Matlab生成圖形，橫坐標(biāo)是隨機(jī)數(shù)的序數(shù)，縱坐標(biāo)是隨機(jī)數(shù)值。很顯然，大部分點在3
σ區(qū)間內(nèi)。在區(qū)間
(
μ
-
σ
,
μ
+
σ)中數(shù)據(jù)點最密。
|
#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)庫rand()函數(shù)的缺陷以及Blitz++隨機(jī)數(shù)生成的簡介.files/image001.jpg)
2
，
下面是一個服從于
μ= 10，σ= 2的正態(tài)分布，根據(jù)3σ法則，大部分點應(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)庫rand()函數(shù)的缺陷以及Blitz++隨機(jī)數(shù)生成的簡介.files/image002.jpg)
|
?
3，生成前述正態(tài)分布的鐘型曲線（PDF）。
這里產(chǎn)生1M的float隨機(jī)數(shù)，然后乘以10轉(zhuǎn)型為整數(shù)，統(tǒng)計在區(qū)間0-170之間隨機(jī)數(shù)出現(xiàn)的概率，這樣再導(dǎo)入Matlab生成圖形就是一個近似的正態(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ù)生成類應(yīng)該說是值得信任的。
[
參考文獻(xiàn)
]
1，?
概率論與數(shù)理統(tǒng)計（
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.