標籤:

C++ 如何生成大隨機數?

由於C++rand()函數的最大值只有0x7FFF,寫了一種方法類似於這樣的:

如果我想從0-999999中隨機一個數,最好的方法當然是直接隨機了,但是做不到,那麼這樣做,從0-999隨機一個數,乘以1000,再加上0-999中隨機一個數,這樣做,跟直接隨機的效果是一樣的嗎?


因為你並不一定知道原來的 rand() 是怎麼實現的,所以這個問題並沒有一定的回答。在現實中,簡單地說,可以這樣拼合,但要注意你的應用。

如果原來的 rand() 函數實現得很理想,就是獨立同均勻分布從 0 到 RAND_MAX,那你當然可以通過把兩個隨機數拼起來得到一個大的隨機數。這種理想情況很容易算出新的隨機變數也是獨立同均勻分布的。

不過事實上我們知道 C 運行時的 rand() 函數通常都是偽隨機數發生器,我們還知道大部分就是線性同餘演算法。更確切地說,最大值 RAND_MAX == 0x7FFF 的是 VC 的運行時,我們甚至可以直接看到實現:

int __cdecl rand (
void
)
{
_ptiddata ptd = _getptd();

return( ((ptd-&>_holdrand = ptd-&>_holdrand * 214013L
+ 2531011L) &>&> 16) 0x7fff );
}

就是一個線性同餘隨機發生器

r_{n+1}= (214013 	imes r_n + 2531011) mod 2^{32}

輸出是取除最高位的高 15 位的結果(32 位整數先右移 16 位,再取低 15 位)。下面的討論如果涉及具體範圍,我們就假定 RAND_MAX 是上面的 32767。

這種偽隨機數兩個拼起來,統計性質多少會所下降——至少我們知道重複周期會變短為原來的一半對吧。不過,經過適當選擇迭代常數的線性同餘演算法,其高維聯合分布的性質仍然是可接受的——這本來也是選擇偽隨機數演算法參數的一個重要指標。VC 這裡使用的庫函數,實際就是 Intel 推薦的一個具體實現,簡單分析可以看:

Fast Random Number Generator on the Intel速 Pentium速 4 Processor

從上面的二維直方圖可以看出,這個函數的二維分布也可以通過直方圖檢驗,因而像你說的這樣使用沒有大的問題。當然,應該是把兩個 15bit 的隨機數拼合成一個 30bit 的來使用。

在實踐中,現代編譯器所選擇的運行時庫,偽隨機數發生器的質量都還不錯,這與多年前的情況不大一樣。當然,如果是對隨機數要求比較高的場合,還是應該做更完整的測試或者選擇更好的隨機數發生器。

但是等等,只要不溢出,這樣當然可以取得很大的隨機數,但卻並不能取得我們需要的隨機數。因為我們實際中使用的隨機數範圍,通常不會那麼巧合就是 [0, 2^{15}-1] 區間的整數,也不會是 [0, 2^{30}-1] 區間的整數。我們需要的,可能是 [100, 155] 之間的整數,也可能是 [5.0, 27.0] 之間的浮點數。也就是說,rand() 函數得到的隨機數只是固定區間上均勻分布的整數,我們需要把它映射為指定區間上指定分布的指定類型的隨機數。這怎麼辦?

指定分布涉及隨機變數的數學變換,我們暫且放在一邊。如果就是均勻分布,如何?

對於均勻分布的整數,比如說 [100, 155],起點不是 0,終點也不是 RAND_MAX,怎麼辦?起點不是 0 好辦,只要做一個加法平移,[0, RAND_MAX] 區間自然就變換為 [a, RAND_MAX+a] 區間了。但區間長度不同,一個是 155-100+1 = 56,一個是 RAND_MAX+1 = 32768,就需要放縮。

人們對區間放縮想出了一些不同的辦法,最簡單的辦法是取模。比如說,100 + rand() % 55 就是一個可以把 rand() 映射到 [100, 155] 區間上的一個辦法。取模的辦法雖然簡單,但並不準確,因為你不能把 32768 按 56 長度等分(32768 / 56 = 585 余 8),那麼必然最後一些數字比前面的出現概率要小一點——準確點兒說,前面 8 個數出現的概率是 586/32768,後面 48 個數出現的概率是 585/32768,後面的比前面的出現概率小 1/32768。

為了避免這種問題,一種準確的辦法就是,取到 rand() 為 [0, 32767] 區間內的整數後,先判斷它是否超過 585*56=32760,如果超過了就捨去它重新取 rand(),否則才對 rand() 值取模(或者除以 585),得到準確均勻分布的隨機數。——這種捨棄、重取的機制在理論上可能導致演算法不終止,因為有可能每次隨機數發生器總返回需要捨棄的數,但演算法終止的概率隨著步驟數而趨於 1(對於演算法生成的偽隨機數一定在有限步終止),並且只要捨棄的數比保留的數少,就能知道需要重複的步驟數期望小於 2。因而這是一個實用的演算法,如 gcc 使用的 libstdc++ 庫,大約就使用這種方式(見下節)。

上面沒有說清楚的是,如果 RAND_MAX 比需要的區間長度還小,那麼就需要回到最開始的問題了:可以用兩個、三個或更多的 rand() 函數,通過 (RAND_MAX+1) 進位數表示的方式拼合成一個大的隨機數。前面已經說了這個大的隨機數的分布實際就是幾個小隨機數的聯合分布,也是均勻的。於是我們把新問題化歸為一個已知的問題。

浮點數的問題看起來稍顯複雜。注意到機器浮點數的乘法放縮是準確的(誤差最小),因此 [a, b) 區間的浮點隨機數 x 是可以通過 [0, 1) 區間的浮點隨機數 u 做線性變換得到,即 x = a + u * (b - a)。因此只要能生成 [0, 1) 區間的浮點隨機數 u 就可以了。如果再知道 IEEE 754 標準中浮點數的定義,就知道一個浮點數由符號位、指數、尾數構成。如果符號位和指數均為 0,尾數部分就恰好是 [0, 1) 區間的數。因此一個 [0, 1) 區間的浮點數就是一個尾數位數的整數。例如對於雙精度浮點數,尾數是 52bit,因此可以生成一個 [0, 2^{52}-1] 區間的隨機整數作為浮點數的尾數,就得到 [0, 1) 區間的隨機浮點數。我們就又把問題化歸為已知的問題了。

做為餐後甜點,我們來看看另一個求 [a, b] 區間整數的辦法:先生成一個 [0, 1) 區間的浮點數,然後用它來生成 [a, b+1) 區間的整數,最後用一個 floor 向負無窮大取整。在數學上,這個辦法無懈可擊(注意對於概率,開區間和閉區間並不影響)。Knuth 在 TAOCP 中就是這樣描述隨機區間生成方式的。但在現實中這樣則不夠準確,因為我們知道,所謂 [0, 1) 區間的雙精度浮點數,其實只是一個 52 比特的整數而已,用它計算區間長度小於 2^{52} 的隨機整數時,因為不能等分區間,就會造成不準確;而用它計算區間長度大於2^{52}的隨機整數時,還會因為區間不夠大,會有一些數永遠取不到。求固定區間浮點隨機數時,我們追求的只是誤差最小;而求隨機整數時,我們就有理由追求完全準確了。

另一個求 [a, b] 區間隨機整數的常見用法是寫

x = a + (double)rand() / (RAND_MAX + 1) * b;

這個辦法其實就是前面說的先求 [0, 1) 區間浮點隨機數,再線性映射到整數區間的辦法,只是 RAND_MAX 通常沒有 52bit 之多,因此精確度還不如前面的方法高。——更準確地說,這個方法的精確程度其實和直接取模差不多是一樣的,只是概率較大和較小的數變得分散了一些而已。

=====================================================================

因為隨機數生成器的使用對於大多數缺少相關專業知識的人來說是比較難的,在 C++ 11 標準中,又增加了 & 頭文件,專門用來生成各種分布隨機數,也可以更換其中的隨機數發生器源。如果你使用支持 C++11 標準庫的新版本的 VC、GCC、Clang 等編譯器,就可以參考 & - C++ Reference 來使用相應的函數。

一方面,你可以不必考慮固定整數範圍的基本的隨機數源怎麼用,直接使用 uniform_int_distribution 和 uniform_real_distribution 類來生成指定範圍的隨機整數或浮點數。

另一方面,你也可以在使用 uniform_int_distribution 之類偽隨機數類生成隨機數時,選擇內部狀態空間更大、周期更長的發生器,如 mt19937(有 19937bit 狀態空間)之類。甚至是基於硬體的真隨機數發生器 random_device。


rand() Considered Harmful

一句話總結:對於非 crypto 級別的隨機數,用 uniform_int_distribution + mt19937 。


別用那個不就行了。Boost裡面有MT法,可以直接生成64位整形。


標準庫已經越來越完美了。建議不要再用C Library

std::random_device rd;
std::uniform_int_distribution& dist(0, 9999999);
std::cout&<&

如果再稍微加一些性能方面的考慮,不去在棧生成random_device.可以考慮封裝成singleton.

//random number generator RNG.hpp
//author : AntiMoron
//Date : 2014-11-27
#ifndef __RNG_HPP__
#define __RNG_HPP__

#include &
#include &

class RNG {
public:
const static std::size_t maxRand = std::random_device::max();
static RNG getInstance() {
static RNG instance;
return instance;
}
std::size_t getInteger() noexcept{
return (*dist)(rd);
}
private:
RNG() noexcept{
regen = std::make_shared&(rd());
dist = std::make_shared(std::uniform_int_distribution& (0,maxRand));
}
std::random_device rd;
std::shared_ptr& regen;
std::shared_ptr& &> dist;
};
#define RAND_INT(x) (RNG::getInstance().getInteger() % x)

#endif

徒手擼代碼不知有錯否。目測沒問題。


其實你想自己寫線性同餘法也沒問題。。。


如果要求不那麼苛刻的話,使用下面這個函數就行

int random(int max_range = 1)
{
if (max_range == 1) return rand() &<&< 16 | rand(); return (rand() &<&< 16 | rand()) % max_range; }


可參考numerical recipes一書。有一章專門講這個。


推薦閱讀:

C++ 如何寫一個函數,使得它的返回值是指向該函數自身的指針?
斷言、異常和返回值的選擇問題?
Qt 為什麼在桌面應用(Windows 平台)中不流行呢?
當面試官問我C++ 11新特性的時候,應該怎樣回答?
既然c++模板元編程是圖靈完全的,那麼刷leetcode的時候如果都用模板元編程,能不能全刷成0ms?

TAG:數學 | C | 隨機 |