mt19937 随机数


一下内容整理自网络资源:

我们讲的随机数其实暗指伪随机数。不少朋友可能想到C语言的rand(),可惜这个函数产生的随机数随机性非常差,而且速度很慢,相信几乎不能胜任一般的应用。

古老的LCG(linear congruential generator)代表了最好的伪随机数产生器算法。主要原因是容易理解,容易实现,而且速度快。这种算法数学上基于X(n+1) = (a * X(n) + c) % m这样的公式,其中:

模m, m > 0

系数a, 0 < a < m

增量c, 0 <= c < m

原始值(种子) 0 <= X(0) < m

其中参数c, m, a比较敏感,或者说直接影响了伪随机数产生的质量。

一般而言,高LCG的m是2的指数次幂(一般2^32或者2^64),因为这样取模操作截断最右的32或64位就可以了。多数编译器的库中使用了该理论实现其伪随机数发生器rand()。下面是部分编译器使用的各个参数值:

Source m a c rand() / Random(L)的种子位
Numerical Recipes 2^32 1664525 1013904223  
Borland C/C++ 2^32 22695477 1 位30..16 in rand(), 30..0 in lrand()
glibc (used by GCC) 2^32 1103515245 12345 位30..0
ANSI C: Watcom, Digital Mars, CodeWarrior, IBM VisualAge C/C++ 2^32 1103515245 12345 位30..16
Borland Delphi, Virtual Pascal 2^32 134775813 1 位63..32 of (seed * L)
Microsoft Visual/Quick C/C++ 2^32 214013 2531011 位30..16
Apple CarbonLib 2^31-1 16807 0 见Park–Miller随机数发生器

LCG不能用于随机数要求高的场合,例如不能用于Monte Carlo模拟,不能用于加密应用。
LCG有一些严重的缺陷,例如如果LCG用做N维空间的点坐标,这些点最多位于m1/n超平面上(Marsaglia定理),这是由于产生的相继X(n)值的关联所致。
另外一个问题就是如果m设置为2的指数,产生的低位序列周期远远小于整体。
一般而言,输出序列的基数b中最低n位,bk = m (k是某个整数),最大周期bn.
有些场合LCG有很好的应用,例如内存很紧张的嵌入式中,电子游戏控制台用的小整数,使用高位可以胜任。
LCG的一种C++实现版本如下:

  1. //************************************************************************  
  2. // Copyright (C) 2008 – 2009 Chipset  
  3. //  
  4. // This program is free software: you can redistribute it and/or modify  
  5. // it under the terms of the GNU Affero General Public License as  
  6. // published by the Free Software Foundation, either version 3 of the  
  7. // License, or (at your option) any later version.  
  8. //  
  9. // but WITHOUT ANY WARRANTY; without even the implied warranty of  
  10. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the  
  11. // GNU Affero General Public License for more details.  
  12. //  
  13. // You should have received a copy of the GNU Affero General Public License  
  14. // along with this program. If not, see .  
  15. //************************************************************************  
  16. #ifndef LCRANDOM_HPP_  
  17. #define LCRANDOM_HPP_  
  18. #include   
  19.   
  20. class lcrandom  
  21. {  
  22. public:  
  23. explicit lcrandom(size_t s = 0) : seed(s)  
  24. {  
  25.     if (0 == seed) seed = std::time(0);  
  26.     randomize();  
  27. }  
  28.   
  29. void reset(size_t s)  
  30. {  
  31.     seed = s;  
  32.     if (0 == seed) seed = std::time(0);  
  33.     randomize();  
  34. }  
  35.   
  36. size_t rand()  
  37. {  
  38. //returns a random integer in the range [0, -1UL)  
  39.     randomize();  
  40.     return seed;  
  41. }  
  42.   
  43. double real()  
  44. {  
  45. //returns a random real number in the range [0.0, 1.0)  
  46.     randomize();  
  47.     return (double)(seed) / -1UL;  
  48. }  
  49.   
  50. private:  
  51. size_t seed;  
  52. void randomize() { seed = 1103515245UL * seed + 12345UL; }  
  53. };  
  54.   
  55. class lcrand_help  
  56. {  
  57. static lcrandom r;  
  58. public:  
  59. lcrand_help() {}  
  60. void operator()(size_t s) { r.reset(s); }  
  61. size_t operator()() const { return r.rand(); }  
  62. double operator()(double) { return r.real(); }  
  63. };  
  64. lcrandom lcrand_help:: r;  
  65.   
  66. extern void lcsrand(size_t s) { lcrand_help()(s); }  
  67. extern size_t lcirand() { return lcrand_help()(); }  
  68. extern double lcdrand() { return lcrand_help()(1.0); }  
  69.   
  70. #endif // LCRANDOM_HPP_  

如果需要高质量的伪随机数,内存充足(约2kb),Mersenne twister算法是个不错的选择。Mersenne twister产生随机数的质量几乎超过任何LCG。不过一般Mersenne twister的实现使用LCG产生种子。

Mersenne twister是Makoto Matsumoto (松本)和Takuji Nishimura (西村)于1997年开发的伪随机数产生器,基于有限二进制字段上的矩阵线性再生。可以快速产生高质量的伪随机数,修正了古老随机数产生算法的很多缺陷。 Mersenne twister这个名字来自周期长度通常取Mersenne质数这样一个事实。常见的有两个变种Mersenne Twister MT19937和Mersenne Twister MT19937-64。

Mersenne Twister有很多长处,例如:周期2^19937 – 1对于一般的应用来说,足够大了,序列关联比较小,能通过很多随机性测试。

关于Mersenne Twister比较详细的论述请参阅
http://www.cppblog.com/Chipset/archive/2009/01/19/72330.html

用Mersenne twister算法实现的伪随机数版本非常多。例如boost库中的高质量快速随机数产生器就是用Mersenne twister算法原理编写的。

C++ : mt19937 Example

C++11 introduces several pseudo-random number generators designed to replace the good-oldrand from the C standard library. I’ll show basic usage examples of std::mt19937, which provides a random number generation based on Mersenne Twister algorithm. Using the Mersenne Twister implementation that comes with C++1 has advantage overrand(), among them:

  1. mt19937 has much longer period than that of rand, e.g. it will take its random sequence much longer to repeat itself.
  2. It much better statistical behavior.
  3. Several different random number generator engines can be initiated simultaneously with different seed, compared with the single “global” seedsrand() provides.

The downside is that mt19937 is a bit less straight-forward to use. However, I hope this post will help with this point :-).

We start with a basic example:

#include 
#include 
 
using namespace std;
 
int main()
{
    mt19937 mt_rand(time(0));
 
    cout << mt_rand() << endl;
 
    return 0;
}

Line 7 creates a new std::mt19937 object and seeds it with time(0). This is just like calling srand(time(0)). Subsequent calls to the newly created object,mt_rand will return a random 32-bit unsigned int.

If you want integers in a specific range, instead all kinds of arithmetics, C++11 provides convinient wrappers:

mt19937::result_type seed = time(0);
auto dice_rand = std::bind(std::uniform_int_distribution(1,6),mt19937(seed));
mt19937::result_type seed = time(0);
auto real_rand = std::bind(std::uniform_real_distribution(0,1),mt19937(seed));

In the first example, each call to
dice_rand() will return an integer between 1 and 6 (inclusive). In the second example each call to
real_rand() will return a double in the range [0,1) (including 0, excluding 1). Note that usage of
std::bind requires
#include .

Finally if you want to go C++11 all the way, you can also replace time(0) with a proper call forstd::chrono when seeding the random number generator:

auto seed = chrono::high_resolution_clock::now().time_since_epoch().count();
mt19937 mt_rand(seed);

The above code requires adding #include to the list of includes.

参考文献:

[1] http://blog.csdn.net/hoxily/article/details/44241387

[2] http://www.guyrutenberg.com/2014/05/03/c-mt19937-example/