這裡介紹如何使用 C++11 標準中內建的亂數函式庫,產生各種機率分布的隨機亂數。
傳統上在 C/C++ 程式中若要產生亂數,大家最常用的就是標準的 rand
函數,它的用法簡單、快速又方便,雖然其功能比較陽春,生成的亂數品質也比較不好,但因為是標準的函數之一,所以還是非常多人喜歡用。
<random>
這個標準的亂數函式庫,它的功能比較完整,生成的隨機變數品質也會比 rand
函數來得好,以下是簡單的使用教學與範例程式碼。
<random>
概念C++11 的 <random>
的主要功能可分為三大類:
<random>
中的隨機設備(std::random_device
)可以讓我們產生接近隨機的亂數種子,讓生成的亂數無法被預測。<random>
中提供了好幾種不同的亂數產生器,例如 std::mt19937
、std::mt19937_64
等,另外也有最陽春的 std::default_random_engine
。std::normal_distribution
)、指數分布(std::exponential_distribution
)、卜瓦松分布(std::poisson_distribution
)等。將以上三大類的工具組合之後,就可以很輕易的產生各種機率分布的隨機亂數。
這是使用 <random>
亂數函式庫產生 [0, 1)
區間連續型均勻分布(uniform distribution)的範例程式:
#include <iostream> #include <random> /* 亂數函式庫 */ #include <ctime> int main(){ /* 亂數產生器 */ std::default_random_engine generator( time(NULL) ); /* 亂數的機率分布 */ std::uniform_real_distribution<float> unif(0.0, 1.0); /* 產生亂數 */ float x = unif(generator); std::cout << "x = " << x << std::endl; return 0; }
這裡我們使用 std::default_random_engine
這個預設的亂數產生器,以時間的資訊作為亂數種子,然後透過 std::uniform_real_distribution
轉換成 [0, 1)
區間的連續型均勻分布。
這個 hello world 程式只是示範最簡單的亂數產生方式,這樣所產生的亂數品質就跟傳統上的 rand
函數差不多(一樣差),接下來我們會使用隨機設備產生更接近隨機的亂數種子,並且將 std::default_random_engine
換成更好一點的亂數產生器,最後套上各種機率分布的轉換,生成各種不同分布的亂數。
一般的亂數產生器都會需要指定亂數種子,如果亂數種子不是隨機的,那麼所產生出來的亂數也就不是隨機的,嚴格來說就算使用時間當作亂數種子,其所產生的亂數也不是隨機的,這種亂數在某些情況下甚至可以被預測。
隨機設備(std::random_device
)是一個均勻分布(uniform distribution)的整數亂數產生器,可以產生比較接近隨機的亂數(例如從硬體上取得隨機的資料),讓人無法預測,以下是 std::random_device
的一些基本用法。
#include <iostream> #include <random> int main(){ /* 隨機設備 */ std::random_device rd; /* 隨機亂數的範圍 */ std::cout << "Min = " << rd.min() << ", Max = " << rd.max() << std::endl; /* 產生隨機的亂數 */ std::cout << "Random Number = " << rd() << std::endl; /* 隨機設備的熵值 */ std::cout << "Entropy = " << rd.entropy() << std::endl; return 0; }
由於 std::random_device
在產生亂數時會消耗掉系統的熵(entropy),所以當熵耗盡時就會影響到亂數的產生,不適合用來產生大量的亂數。
實務上我們可以使用 std::random_device
所產生的亂數來做為亂數種子,這樣產生的亂數會比較有隨機性:
#include <iostream> #include <random> int main(){ /* 隨機設備 */ std::random_device rd; /* 亂數產生器 */ std::default_random_engine generator( rd() ); /* 亂數的機率分布 */ std::uniform_real_distribution<float> unif(0.0, 1.0); /* 產生亂數 */ float x = unif(generator); std::cout << "x = " << x << std::endl; return 0; }
偽隨機亂數產生器會依照給定的亂數種子,「計算」出一連串類似隨機的亂數,<random>
中提供好幾種不同演算法的亂數產生器:
亂數產生器 | 說明 |
---|---|
std::minstd_rand0 |
Minimal standard,x = x * 16807 % 2147483647 |
std::minstd_rand |
Minimal standard,x = x * 48271 % 2147483647 |
std::mt19937 |
梅森旋轉演算法 |
std::mt19937_64 |
梅森旋轉演算法 |
std::ranlux24_base |
Subtract with carry |
std::ranlux48_base |
Subtract with carry |
std::ranlux24 |
24-bit RANLUX generator by Martin Lüscher and Fred James, 1994 |
std::ranlux48 |
48-bit RANLUX generator by Martin Lüscher and Fred James, 1994 |
std::knuth_b |
Knuth-B generator |
std::default_random_engine |
預設亂數產生器,通常是 LCG 類型的。 |
各種亂數產生器的用法都類似,只要將前面範例中的 std::default_random_engine
抽換掉即可:
#include <iostream> #include <random> int main(){ std::random_device rd; /* 梅森旋轉演算法 */ std::mt19937 generator( rd() ); std::uniform_real_distribution<float> unif(0.0, 1.0); float x = unif(generator); std::cout << "x = " << x << std::endl; return 0; }
其他的用法以此類推,若要查詢亂數的範圍,也是使用 min
與 max
兩個方法函數:
/* 產生亂數的範圍 */ std::cout << "Min = " << generator.min() << ", Max = " << generator.max() << std::endl;
講究亂數品質的人,請避免使用 std::minstd_rand0
、std::minstd_rand
與 std::default_random_engine
這些太簡略的演算法,可以改用 std::mt19937
這類比較新的演算法。
亂數產生器只能產生固定範圍的整數亂數,若我們需要其他類型的亂數,就要靠機率分佈(probability distribution)的函數將基本的離散型均勻分布轉換成我們需要的分布。
以下是 <random>
中所提供的機率分布函數。
機率分布函數的使用也是非常單純,將前面範例中的 std::uniform_real_distribution
置換成自己需要的分布即可,另外再指定機率分布的參數即可:
#include <iostream> #include <random> int main(){ std::random_device rd; std::mt19937 generator( rd() ); /* 標準常態分布 */ std::normal_distribution<double> norm(0.0, 1.0); /* 產生標準常態分布的亂數 */ double x = norm(generator); std::cout << "x = " << x << std::endl; return 0; }
常態分布需要的參數是平均值與標準差,這裡我們產生的是標準常態分布(平均值為 0
、標準差為 1
)的的亂數。
以下是一個產生常態分布亂數,並畫出簡易直方圖的範例:
#include <iostream> #include <random> int main(){ std::random_device rd; std::mt19937 generator( rd() ); /* 常態分布 */ std::normal_distribution<double> norm(5.0, 2.0); /* 常態分布亂數直方圖 */ const int n = 10000; const int nstars = 100; int p[10]={}; for (int i=0; i < n; ++i) { double x = norm(generator); if ((x >= 0.0) && (x < 10.0)) ++p[int(x)]; } for (int i=0; i<10; ++i) { std::cout << i << ": "; std::cout << std::string(p[i]*nstars/n, '*') << std::endl; } return 0; }
執行後的輸出會類似這樣:
0: * 1: **** 2: ********* 3: *************** 4: ******************* 5: ******************* 6: ************** 7: ******** 8: **** 9: *
只要將分布函數替換後,稍微修改一下分布的參數,就可以畫出各種分布的直方圖,以下是指數分布的範例:
#include <iostream> #include <random> int main(){ std::random_device rd; std::mt19937 generator( rd() ); /* 指數分布 */ std::exponential_distribution<double> exp(0.3); /* 指數分布亂數直方圖 */ const int n = 10000; const int nstars = 100; int p[10]={}; for (int i=0; i < n; ++i) { double x = exp(generator); if ((x >= 0.0) && (x < 10.0)) ++p[int(x)]; } for (int i=0; i<10; ++i) { std::cout << i << ": "; std::cout << std::string(p[i]*nstars/n, '*') << std::endl; } return 0; }
執行後的輸出為:
0: ************************** 1: ******************* 2: ************* 3: ********** 4: ******** 5: ***** 6: **** 7: *** 8: ** 9: *
參考資料:GeeksforGeeks、GeeksforGeeks、StackOverflow、dreamincode、PCG