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