隨機數庫:生成器

簡介

本庫提供了數個偽隨機數生成器。偽隨機數生成器的性質與其算法和參數都有密切的關係。生成算法用取數值參數的類模板實現,封裝在 boost::random 名字空間中。預先選擇好的一組參數通過用 typedef 進行專門化來給出,置於 boost 名字空間中。

程序執行過程中不應當頻繁構造偽隨機數生成器。原因有二:其一,生成器的構造要求對其內部狀態進行完全的初始化;因此包含大量內部狀態的生成器 (詳下) 的構造往往是昂貴的。其二,生成器的構造需要提供序列的「種子」。取得多個良好的「種子」值往往是困難的。取得種子的方法之一是使用精確度盡可能高 (毫秒甚至納秒) 的系統時間。如果再構建一個生成器之後再構建另一個,那麼它們所使用的種子值之差會很接近一個常數。若系統時間的精確度不夠高,這個差值甚至為零,從而得到兩個行為相同的生成器。對於很多應用,這樣的行為都是不合適的。

注意以下描述的所有偽隨機數生成器都是 CopyConstructible and Assignable。生成器的複製或賦值將會複製其內部狀態,因此原有生成器和副本將會生成完全相同的隨機數序列。這種行為往往不是需要的。特別注意標準庫中的算法,如 std::generate。它們取函子的值作為參數,因此會調用複製構造函數。

下表給出了生成器某些性質的概覽。循環長度是給出了生成器性質的粗略描述;相對速度是性能的量度,值越大表示隨機數生成越快。

生成器 循環長度 內存需求 (近似) 相對速度 (近似) 備註
minstd_rand 231-2 sizeof(int32_t) 40 -
rand48 248-1 sizeof(uint64_t) 80 -
lrand48 (C library) 248-1 - 20 狀態由全局變量保存
ecuyer1988 約 261 2*sizeof(int32_t) 20 -
kreutzer1986 ? 1368*sizeof(uint32_t) 60 -
hellekalek1995 231-1 sizeof(int32_t) 3 good uniform distribution in several dimensions
mt11213b 211213-1 352*sizeof(uint32_t) 100 good uniform distribution in up to 350 dimensions
mt19937 219937-1 625*sizeof(uint32_t) 100 good uniform distribution in up to 623 dimensions
lagged_fibonacci607 約 232000 607*sizeof(double) 150 -
lagged_fibonacci1279 約 267000 1279*sizeof(double) 150 -
lagged_fibonacci2281 約 2120000 2281*sizeof(double) 150 -
lagged_fibonacci3217 約 2170000 3217*sizeof(double) 150 -
lagged_fibonacci4423 約 2230000 4423*sizeof(double) 150 -
lagged_fibonacci9689 約 2510000 9689*sizeof(double) 150 -
lagged_fibonacci19937 約 21050000 19937*sizeof(double) 150 -
lagged_fibonacci23209 約 21200000 23209*sizeof(double) 140 -
lagged_fibonacci44497 約 22300000 44497*sizeof(double) 60 -

可以從上表中觀察到,選擇隨機數生成器時往往要在品質、性能、內存三者之間作出權衡。本庫中提供的大量生成器足以讓程序員選擇最適合自己應用領域的一個。此外,在 Monte Carlo 模擬中使用多個有著根本性不同的隨機數生成器將會提高結果的可信度。

如果你不知如何選擇,不妨先使用 mt19937:它快速、可靠。

注意:以上隨機數生成器不能用於要求使用不確定隨機數的應用領域;關於不確定隨機數生成器,參看 nondet_random.html

對於在 概念文檔 中已經給出定義的成員,在下面的描述中不會詳細說明。

<boost/random.hpp> 中的生成器:概覽

namespace boost {
  namespace random {
    template<class IntType, IntType m>
    class const_mod;
    template<class IntType, IntType a, IntType c, IntType m, IntType val>
    class linear_congruential;
  }
  class rand48;
  typedef random::linear_congruential< /* ... */ > minstd_rand0;
  typedef random::linear_congruential< /* ... */ > minstd_rand;

  namespace random {
    template<class DataType, int w, int n, int m, int r, DataType a, int u,
        int s, DataType b, int t, DataType c, int l, IntType val>
    class mersenne_twister;
  }
  typedef random::mersenne_twister< /* ... */ > mt11213b;
  typedef random::mersenne_twister< /* ... */ > mt19937;

  namespace random {
    template<class FloatType, unsigned int  p, unsigned int q>
    class lagged_fibonacci;
  }
  typedef random::lagged_fibonacci< /* ... */ > lagged_fibonacci607;
  typedef random::lagged_fibonacci< /* ... */ > lagged_fibonacci1279;
  typedef random::lagged_fibonacci< /* ... */ > lagged_fibonacci2281;
  typedef random::lagged_fibonacci< /* ... */ > lagged_fibonacci3217;
  typedef random::lagged_fibonacci< /* ... */ > lagged_fibonacci4423;
  typedef random::lagged_fibonacci< /* ... */ > lagged_fibonacci9689;
  typedef random::lagged_fibonacci< /* ... */ > lagged_fibonacci19937;
  typedef random::lagged_fibonacci< /* ... */ > lagged_fibonacci23209;
  typedef random::lagged_fibonacci< /* ... */ > lagged_fibonacci44497;  
} // namespace boost

random::const_mod 類模板

概覽

template<class IntType, IntType m>
class random::const_mod
{
public:
  template<IntType c>
  static IntType add(IntType x);

  template<IntType a>
  static IntType mult(IntType x);

  template<IntType a, IntType c>
  static IntType mult_add(IntType x);

  static IntType invert(IntType x);
private:
  const_mod();         // don't instantiate
};

描述

const_mod 類模板提供了同余運算 (模運算) 函數;這些函數小心設計,確保不會造成溢出。所有的成員函數都是靜態的;const_mod<> 不允許實例化。

模板參數 IntType 為某一整數類型,m 為模數 (modulus)。

注記:以下資料給出了一個方法,允許 m 值很大時仍能快速進行同余乘法:

"A more portable FORTRAN random number generator", Linus Schrage, ACM Transactions on Mathematical Software, Vol. 5, No. 2, June 1979, pp. 132-138

成員函數

template<IntType c> static IntType add(IntType x)

返回: (x+c) mod m

template<IntType a> static IntType mult(IntType x)

返回: (a*x) mod m

template<IntType a, IntType c> static IntType
mult_add(IntType x)

返回: (a*x+c) mod m

static IntType invert(IntType x)

返回: 滿足 (a*i) mod m == 1 的 i
前條件 m 為質數

random::linear_congruential 類模板

概覽

#include <boost/random/linear_congruential.hpp>

template<class IntType, IntType a, IntType c, IntType m, IntType val>
class linear_congruential
{
public:
  typedef IntType result_type;
  static const IntType multiplier = a;
  static const IntType increment = c;
  static const IntType modulus = m;
  static const bool has_fixed_range = true;
  static const result_type min_value;
  static const result_type max_value;
  explicit linear_congruential_fixed(IntType x0 = 1);
  // compiler-generated copy constructor and assignment operator are fine
  void seed(IntType x0);
  IntType operator()();
};

typedef random::linear_congruential<long, 16807L, 0, 2147483647L,
     1043618065L> minstd_rand0;
typedef random::linear_congruential<long, 48271L, 0, 2147483647L,
     399268537L> minstd_rand;

描述

linear_congruential 類模板的實例是 偽隨機數生成器 的模型。以下資料給出了線性同余 (linear congruential) 偽隨機數生成器的描述:

"Mathematical methods in large-scale computing units", D. H. Lehmer, Proc. 2nd Symposium on Large-Scale Digital Calculating Machines, Harvard University Press, 1951, pp. 141-146
令 x(n) 為偽隨機數生成器的返回值序列。對於線性同餘生成器,有 x(n+1) := (a * x(n) + c) mod m。生成器的參數是 x(0), a, c, m。模板參數 IntType 為某一整數類型,a, c, m 須在其可表示範圍內。a, c 應小於 m。

注意: 此生成器的品質很大程度上取決於參數的選擇。用戶代碼應該使用已經選擇好參數的生成器,如 minstd_rand
對參數 a, c, m 的每一組選擇,都附帶定義了一些類別,從而 static 成員不會違反單一定義原則。

成員

explicit linear_congruential(IntType x0 = 1)

效果:構造一 linear_congruential 生成器,令其 x(0) := x0

void seed(IntType x0)

效果:令生成器當前的 x(n) 值為 x0

專門化

專門化 minstd_rand0 來自以下資料:

A pseudo-random number generator for the System/360, P.A. Lewis, A.S. Goodman, J.M. Miller, IBM Systems Journal, Vol. 8, No. 2, 1969, pp. 136-146
以下資料中有 minstd_rand0 以及 minstd_rand 的討論:
"Random Number Generators: Good ones are hard to find", Stephen K. Park and Keith W. Miller, Communications of the ACM, Vol. 31, No. 10, October 1988, pp. 1192-1201

rand48

概覽

#include <boost/random/linear_congruential.hpp>

class rand48 
{
public:
  typedef int32_t result_type;
  static const bool has_fixed_range = true;
  static const int32_t min_value = 0;
  static const int32_t max_value = 0x7fffffff;
  
  explicit rand48(int32_t x0 = 1);
  explicit rand48(uint64_t x0);
  // compiler-generated copy ctor and assignment operator are fine
  void seed(int32_t x0);
  void seed(uint64_t x0);
  int32_t operator()();
};

描述

rand48 類是 偽隨機數生成器的模型。它使用線性同余算法,其中參數 a = 0x5DEECE66D, c = 0xB, m = 2**48。它能產生與某些系統上的 lrand48() 函數相同的結果 (假設沒有調用 lcong48)。

該類中的常量和枚舉定義等使用了 uint64_t;如果系統沒有提供整數類型 uint64_t,該類不可用。

構造函數

rand48(int32_t x0)

效果:構造一 rand48 生成器,令其 x(0) := (x0 << 16) | 0x330e。

rand48(uint64_t x0)

效果: 構造一 rand48 生成器,令其 x(0) := x0

賦種

void seed(int32_t x0)

效果: 令生成器當前的 x(n) 值為 (x0 << 16) | 0x330e。

void seed(uint64_t x0)

效果: 令生成器當前的 x(n) 值為 x0

random::additive_combine類模板

概覽

#include <boost/random/additive_combine.hpp>

template<class MLCG1, class MLCG2, typename MLCG1::result_type val>
class random::additive_combine
{
public:
  typedef MLCG1 first_base;
  typedef MLCG2 second_base;
  typedef typename MLCG1::result_type result_type;
  static const bool has_fixed_range = true;
  static const result_type min_value = 1;
  static const result_type max_value = MLCG1::max_value-1;
  additive_combine();
  additive_combine(typename MLCG1::result_type seed1, 
                   typename MLCG2::result_type seed2);
  result_type operator()();
  bool validation(result_type x) const;
};

typedef random::additive_combine<
    random::linear_congruential<int32_t, 40014, 0, 2147483563, 0>,
    random::linear_congruential<int32_t, 40692, 0, 2147483399, 0>,
  /* unknown */ 0> ecuyer1988;

描述

additive_combine 類模板的實例是 偽隨機數生成器 的模型。它把兩個倍增線性同餘生成器 (即 c = 0 的生成器) 組合起來。以下資料給出了其描述:

"Efficient and Portable Combined Random Number Generators", Pierre L'Ecuyer, Communications of the ACM, Vol. 31, No. 6, June 1988, pp. 742-749, 774
模板參數 MLCG1MLCG2 應為兩個不同的線性同餘生成器,且都有 c = 0。每次調用返回 X(n) := (MLCG1(n) - MLCG2(n)) mod (m1 - 1),其中 m1 為 MLCG1 的模數。

模板參數 valvalidation 檢驗的參考值。

成員

additive_combine()

效果:用兩個基礎生成器的默認構造函數來構造一 additive_combine 生成器。

additive_combine(typename MLCG1::result_type seed1, 
                 typename MLCG2::result_type seed2)

效果:seed1seed2 分別為兩基礎構造函數的參數來構造一 additive_combine 生成器。

專門化

專門化 ecuyer1988 在上面的資料中給出。

random::shuffle_output 類模板

概覽

#include <boost/random/shuffle_output.hpp>

template<class UniformRandomNumberGenerator, int k, 
  typename UniformRandomNumberGenerator::result_type val = 0>
class random::shuffle_output
{
public:
  typedef UniformRandomNumberGenerator base_type;
  typedef typename base_type::result_type result_type;
  static const bool has_fixed_range = false;

  shuffle_output();
  template<class T> explicit shuffle_output(T seed);
  explicit shuffle_output(const base_type & rng);
  template<class T> void seed(T s);

  result_type operator()();
  result_type min() const;
  result_type max() const;
  bool validation(result_type) const;
};

描述

shuffle_output 類模板的實例是 偽隨機數生成器 的模型。它把某種均勻隨機數生成器 (通常是線性同餘生成器) 的多個輸出混合起來,以獲得較好的統計學性質。據 Donald E. Knuth 在 "The Art of Computer Programming, Vol. 2" 中所載,這一算法在以下資料中給出:

"Improving a poor random number generator", Carter Bays and S.D. Durham, ACM Transactions on Mathematical Software, Vol. 2, 1979, pp. 59-64.
基礎生成器的輸出緩衝在長度為 k 的數組中。每一個輸出 X(n) 另有一個用途:它給出了 X(n+1) 在數組中的位置。使用過的數組元素由基礎生成器產生的新輸出取代。

模板參數為基礎生成器和數組長度 k,其中 k 一般取 100 左右的值。模板參數 valvalidation 檢驗的參考值。

成員

shuffle_output()

效果:用基礎生成器的默認構造函數構造一 shuffle_output 生成器。

複雜度:基礎生成器的 k+1 次調用。

template<class T> explicit shuffle_output(T seed)

效果:seed 為基礎生成器單參數構造函數的參數構造一 shuffle_output 生成器。

複雜度:基礎生成器的 k+1 次調用。

explicit shuffle_output(const base_type & rng)

前條件: 模板參數 UniformRandomNumberGenerator 為一 CopyConstructible 類型。

效果: 通過複製傳遞的生成器構造一 shuffle_output 生成器。

複雜度:基礎生成器的 k+1 次調用。

template<class T> void seed(T s)

效果:seed 為參數調用基礎生成器的 seed 方法,並對內部緩衝數組重新初始化。

複雜度:基礎生成器的 k+1 次調用。

專門化

據 Harry Erwin 所述 (私人 e-mail),以下資料給出了專門化 kreutzer1986 的描述:

"System Simulation: programming Styles and Languages (International Computer Science Series)", Wolfgang Kreutzer, Addison-Wesley, December 1986.

random::inversive_congruential 類模板

概覽

#include <boost/random/inversive_congruential.hpp>

template<class IntType, IntType a, IntType b, IntType p>
class random::inversive_congruential
{
public:
  typedef IntType result_type;
  static const bool has_fixed_range = true;
  static const result_type min_value = (b == 0 ? 1 : 0);
  static const result_type max_value = p-1;
  static const result_type multiplier = a;
  static const result_type increment = b;
  static const result_type modulus = p;
  explicit inversive_congruential(IntType y0 = 1);
  void seed(IntType y0);
  IntType operator()();
};

typedef random::inversive_congruential<int32_t, 9102, 2147483647-36884165, 2147483647> hellekalek1995;

描述

inversive_congruential 類模板的實例是 偽隨機數生成器 的模型。它使用了逆同余算法 (inversive congruential algorithm, ICG)。以下資料給出了這一算法的描述:

"Inversive pseudorandom number generators: concepts, results and links", Peter Hellekalek, In: "Proceedings of the 1995 Winter Simulation Conference", C. Alexopoulos, K. Kang, W.R. Lilegdon, and D. Goldsman (editors), 1995, pp. 255-262. ftp://random.mat.sbg.ac.at/pub/data/wsc95.ps
輸出序列由 x(n+1) = (a*inv(x(n)) - b) (mod p) 定義,其中 x(0), a, b 和質數 p 是生成器的參數。inv(k) 為模 p 整數域中 k 的乘法逆元素。特殊地,inv(0) := 0。

模板參數 IntType 應為一帶符號整數類型,p 在其可表示範圍內。a, b 和 p 是生成器的參數。

注記:目前的實現使用歐幾里德算法 (Euclidian Algorithm) 計算乘法逆元素。因此,逆同餘生成器比其它生成器慢 10-20 倍 (參見 performance 一節)。然而,上面給出的資料指出這一生成器應當只慢 3 倍左右,因此歐幾里德算法可能不是計算乘法逆元素的最優算法。

成員

inversive_congruential(IntType y0 = 1)

效果:構造一 inversive_congruential 生成器,其初始狀態為 y0

void seed(IntType y0)

效果:令當前狀態為 y0

專門化

上面的資料給出了專門化 hellekalek1995 的描述。

random::mersenne_twister 類模板

概覽

#include <boost/random/mersenne_twister.hpp>

template<class DataType, int w, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, IntType val>
class random::mersenne_twister
{
public:
  typedef DataType result_type;
  static const bool has_fixed_range = true;
  static const result_type min_value;
  static const result_type max_value;
  mersenne_twister();
  explicit mersenne_twister(DataType value);
  template<class Generator> explicit mersenne_twister(Generator & gen);
  // compiler-generated copy ctor and assignment operator are fine
  void seed();
  void seed(DataType value);
  template<class Generator> void seed(Generator & gen);
  result_type operator()();
  bool validation(result_type) const;
};

typedef mersenne_twister<uint32_t,351,175,19,0xccab8ee7,11,7,0x31b6ab00,15,0xffe50000,17, /* unknown */ 0> mt11213b;
typedef mersenne_twister<uint32_t,624,397,31,0x9908b0df,11,7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937;

描述

mersenne_twister 的實例是 偽隨機數生成器 的模型。以下資料給出了其算法的描述:

"Mersenne Twister: A 623-dimensionally equidistributed uniform pseudo-random number generator", Makoto Matsumoto and Takuji Nishimura, ACM Transactions on Modeling and Computer Simulation: Special Issue on Uniform Random Number Generation, Vol. 8, No. 1, January 1998, pp. 3-30.
注記:boost 庫中的變種是重新實現的,不派生自也不使用以上網站提供的 mt19937.c。不過兩個生成器的等價性是經過檢驗的。
2005 年 4 月發現算法的一個 缺陷 後修改了用整數賦種的方法。
此生成器的品質很大程度上取決於參數的選擇。用戶代碼應該使用已經選擇好參數的生成器,如 mt19937
此生成器的狀態數組需要的內存是可觀的。例如 mt11213b 需要約 1408 bytes 內存,而 mt19937 需要約 2496 bytes。

構造函數

mersenne_twister()

效果: 構造一 mersenne_twister 並調用 seed()

explicit mersenne_twister(result_type value)

效果: 構造一 mersenne_twister 並調用 seed(value)

template<class Generator> explicit mersenne_twister(Generator & gen)

效果: 構造一 mersenne_twister 並調用 seed(gen)

注意:當使用左值進行直接初始化時 (如 Gen gen2(gen);),編譯器將選擇這一構造函數,而不是複製構造函數。若要複製另一 mersenne_twister 的狀態,使用 Gen gen2 = gen;;這是複製初始化的語法,將會調用複製構造函數。

賦種

void seed()

效果:調用 seed(result_type(5489))

void seed(result_type value)

效果:令狀態 x(0) 為 v mod 2w。然後迭代地,
令 x(i) 為 (i + 1812433253 * (x(i-1) xor (x(i-1) rshift w-2))) mod 2w for i = 1 .. n-1. x(n) 為 operator() 返回的第一個值。

template<class Generator> void seed(Generator & gen)

效果:mersenne_twister 的當前狀態為調用gen n 次的結果序列。

複雜度:genn 次調用。

注意:當用左值調用 seed 時,函數重載的規則都將選擇這一函數模板,除非參數的類型和 result_type 相同。如果要用其它的整數類型作參數,你應該先把它顯式轉換為 result_type

專門化

專門化 mt11213bmt19937 來自上面引用的資料。

random::lagged_fibonacci 類模板

概覽

#include <boost/random/lagged_fibonacci.hpp>

template<class FloatType, unsigned int p, unsigned int q>
class lagged_fibonacci
{
public:
  typedef FloatType result_type;
  static const bool has_fixed_range = false;
  static const unsigned int long_lag = p;
  static const unsigned int short_lag = q;
  result_type min() const { return 0.0; }
  result_type max() const { return 1.0; }
  lagged_fibonacci();
  explicit lagged_fibonacci(uint32_t value);
  template<class Generator>
  explicit lagged_fibonacci(Generator & gen);
  // compiler-generated copy ctor and assignment operator are fine
  void seed(uint32_t value = 331u);
  template<class Generator> void seed(Generator & gen);
  result_type operator()();
  bool validation(result_type x) const;
};

typedef random::lagged_fibonacci<double, 607, 273> lagged_fibonacci607;
typedef random::lagged_fibonacci<double, 1279, 418> lagged_fibonacci1279;
typedef random::lagged_fibonacci<double, 2281, 1252> lagged_fibonacci2281;
typedef random::lagged_fibonacci<double, 3217, 576> lagged_fibonacci3217;
typedef random::lagged_fibonacci<double, 4423, 2098> lagged_fibonacci4423;
typedef random::lagged_fibonacci<double, 9689, 5502> lagged_fibonacci9689;
typedef random::lagged_fibonacci<double, 19937, 9842> lagged_fibonacci19937;
typedef random::lagged_fibonacci<double, 23209, 13470> lagged_fibonacci23209;
typedef random::lagged_fibonacci<double, 44497, 21034> lagged_fibonacci44497;

描述

lagged_fibonacci 類模板的實例是 偽隨機數生成器 的模型。It uses a lagged Fibonacci algorithm with two lags p and q, evaluated in floating-point arithmetic: x(i) = x(i-p) + x(i-q) (mod 1) with p > q. See

"Uniform random number generators for supercomputers", Richard Brent, Proc. of Fifth Australian Supercomputer Conference, Melbourne, Dec. 1992, pp. 704-706.

注意:此生成器的品質很大程度上取決於參數的選擇。用戶代碼應該使用已經選擇好參數的生成器,如 lagged_fibonacci607
此生成器的狀態數組需要的內存是可觀的。例如 lagged_fibonacci607 需要約 4856 bytes 內存,而 lagged_fibonacci44497 需要約 350 KBytes。

構造函數

lagged_fibonacci()

效果:構造一 lagged_fibonacci 生成器並調用 seed()

explicit lagged_fibonacci(uint32_t value)

效果:構造一 lagged_fibonacci 生成器並調用 seed(value)

template<class Generator> explicit lagged_fibonacci(Generator & gen)

效果:構造一 lagged_fibonacci 生成器並調用 seed(gen)

賦種

void seed()

效果:調用 seed(331u)

void seed(uint32_t value)

效果:value 作參數構造一 minstd_rand0,並用它作參數調用 seed

template<class Generator> void seed(Generator & gen)

效果:lagged_fibonacci 的狀態為對 uniform_01<gen, FloatType>p 次調用的結果序列。
複雜度:genp 次調用。

專門化

專門化 lagged_fibonacci607 ... lagged_fibonacci44497 (見上) use well tested lags. (參考資料待添加)

性能

測試程序 random_speed.cpp 用一個緊湊的循環測試 random.hpp 中上述算法的實現的運行速度。測試環境:Pentium Pro 200 MHz with gcc 2.95.2, Linux 2.2.13, glibc 2.1.2。

每次調用耗時 [微秒, usec]
rand48 0.096
rand48 run-time configurable 0.697
lrand48 glibc 2.1.2 0.844
minstd_rand 0.174
ecuyer1988 0.445
kreutzer1986 0.249
hellekalek1995 (inversive) 4.895
mt11213b 0.165
mt19937 0.165
mt19937 original 0.185
lagged_fibonacci607 0.111
lagged_fibonacci4423 0.112
lagged_fibonacci19937 0.113
lagged_fibonacci23209 0.122
lagged_fibonacci44497 0.263

誤差範圍為 +/- 10 納秒 (nsec)。


Valid HTML 4.01 Transitional

Revised 05 December, 2006

Copyright © 2000-2005 Jens Maurer

Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

中文版修訂:2008/6/28

Copyright © 2008 xiaq

在 Boost Software License, Version 1.0 的條款下發佈。(參看文件 LICENSE_1_0.txt 或在線副本 http://www.boost.org/LICENSE_1_0.txt)