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