随机数算法

最新推荐文章于 2022-07-10 20:16:50 发布

sdfgh2046

最新推荐文章于 2022-07-10 20:16:50 发布

阅读量3.4k

点赞数 1

分类专栏：计算机语言文章标签：算法 generator list algorithm microsoft transformation

计算机语言专栏收录该内容

23 篇文章 0 订阅

订阅专栏

Mersenne twister -- 目前为止最好的随机数算法

Mersenne twister -- From Wikipedia, the free encyclopedia

The Mersenne twister is a pseudorandom number generator developed in 1997 by Makoto Matsumoto ( 松本眞 , Makoto Matsumoto ^? ) and Takuji Nishimura ( 西村拓士 , Takuji Nishimura ^? )^[1] that is based on a matrix linear recurrence over a finite binary field F ₂ . It provides for fast generation of very high-quality pseudorandom numbers, having been designed specifically to rectify many of the flaws found in older algorithms.

Its name derives from the fact that period length is chosen to be a Mersenne prime . There are at least two common variants of the algorithm, differing only in the size of the Mersenne primes used. The newer and more commonly used one is the Mersenne Twister MT19937, with 32-bit word length. There is also a variant with 64-bit word length, MT19937-64, which generates a different sequence.

Application

Unlike Blum Blum Shub , the algorithm in its native form is not suitable for cryptography . Observing a sufficient number of iterates (624 in the case of MT19937) allows one to predict all future iterates.

Another issue is that it can take a long time to turn a non-random initial state into output that passes randomness tests , due to its size. A small lagged Fibonacci generator or linear congruential generator gets started much quicker and usually is used to seed the Mersenne Twister. If only a few numbers are required and standards aren't high it is simpler to use the seed generator. But the Mersenne Twister will still work.

For many other applications, however, the Mersenne twister is quickly becoming the pseudorandom number generator of choice^{[citation needed ]} . Since the library is portable, freely available and quickly generates good quality pseudorandom numbers it is rarely a bad choice.

It is designed with Monte Carlo simulations and other statistical simulations in mind. Researchers primarily want good quality numbers but also benefit from its speed and portability.

The commonly used variant of Mersenne Twister, MT19937 has the following desirable properties:

It was designed to have a period of 2¹⁹⁹³⁷ − 1 (the creators of the algorithm proved this property). In practice, there is little reason to use larger ones, as most applications do not require 2¹⁹⁹³⁷ unique combinations (2¹⁹⁹³⁷ is approximately 4.3 × 10⁶⁰⁰¹ ).
It has a very high order of dimensional equidistribution (see linear congruential generator ). This implies that there is negligible serial correlation between successive values in the output sequence.
It passes numerous tests for statistical randomness, including the Diehard tests . It passes most, but not all, of the even more stringent TestU01 Crush randomness tests.

The Mersenne Twister algorithm has received some criticism in the computer science field, notably by George Marsaglia . These critics claim that while it is good at generating random numbers, it is not very elegant and is overly complex to implement. Marsaglia has provided several examples of random number generators that are less complex yet which he claims provide significantly larger periods. For example, a simple complementary multiply-with-carry generator can have a period 10³³⁰⁰⁰ times as long, be significantly faster, and maintain better or equal randomness.^{[2] [3]}

Algorithmic detail

The Mersenne Twister algorithm is a twisted generalised feedback shift register ^[4] (twisted GFSR, or TGFSR) of rational normal form (TGFSR(R)), with state bit reflection and tempering. It is characterized by the following quantities:

w : word size (in number of bits)
n : degree of recurrence
m : middle word, or the number of parallel sequences, 1 ≤ m ≤ n
r : separation point of one word, or the number of bits of the lower bitmask, 0 ≤ r ≤ w - 1
a : coefficients of the rational normal form twist matrix
b , c : TGFSR(R) tempering bitmasks
s , t : TGFSR(R) tempering bit shifts
u , l : additional Mersenne Twister tempering bit shifts

with the restriction that 2^nw ^− r − 1 is a Mersenne prime. This choice simplifies the primitivity test and k -distribution test that are needed in the parameter search.

For a word x with w bit width, it is expressed as the recurrence relation

with | as the bitwise or and ⊕ as the bitwise exclusive or (XOR), x ^u , x ^l being x with upper and lower bitmasks applied. The twist transformation A is defined in rational normal form

with I_n _− 1 as the (n − 1) × (n − 1) identity matrix (and in contrast to normal matrix multiplication, bitwise XOR replaces addition). The rational normal form has the benefit that it can be efficiently expressed as

where

In order to achieve the 2^nw ^− r − 1 theoretical upper limit of the period in a TGFSR, φ_B (t ) must be a primitive polynomial , φ_B (t ) being the characteristic polynomial of

The twist transformation improves the classical GFSR with the following key properties:

Period reaches the theoretical upper limit 2^nw ^− r − 1 (except if initialized with 0)
Equidistribution in n dimensions (e.g. linear congruential generators can at best manage reasonable distribution in 5 dimensions)

As like TGFSR(R), the Mersenne Twister is cascaded with a tempering transform to compensate for the reduced dimensionality of equidistribution (because of the choice of A being in the rational normal form), which is equivalent to the transformation A = R → A = T ⁻¹ RT , T invertible. The tempering is defined in the case of Mersenne Twister as

y := x ⊕ (x >> u )

y := :y ⊕ ((y << s ) & b )

y := :y ⊕ ((y << t ) & c )

z := y ⊕ (y >> l )

with <<, >> as the bitwise left and right shifts, and & as the bitwise and . The first and last transforms are added in order to improve lower bit equidistribution. From the property of TGFSR, is required to reach the upper bound of equidistribution for the upper bits.

The coefficients for MT19937 are:

(w , n , m , r ) = (32, 624, 397, 31)
a = 9908B0DF₁₆
u = 11
(s , b ) = (7, 9D2C5680₁₆ )
(t , c ) = (15, EFC60000₁₆ )
l = 18

Pseudocode

The following generates uniformly 32-bit integers in the range [0, 2³² − 1] with the MT19937 algorithm:

// Create a length 624 array to store the state of the generator

int [0..623] MT

int index = 0

// Initialize the generator from a seed

function initializeGenerator(int seed) {

MT[0] := seed

for i from 1 to 623 { // loop over each other element

MT[i] := last 32 bits of (1812433253 * (MT[i-1] xor (right shift by 30 bits (MT[i-1]))) + i) // 0x6c078965

}

// Extract a tempered pseudorandom number based on the index-th value,

// calling generateNumbers() every 624 numbers

function extractNumber() {

if index == 0 {

generateNumbers()

}

int y := MT[index]

y := y xor (right shift by 11 bits (y))

y := y xor (left shift by 7 bits (y) and (2636928640)) // 0x9d2c5680

y := y xor (left shift by 15 bits (y) and (4022730752)) // 0xefc60000

y := y xor (right shift by 18 bits (y))

index := (index + 1) mod 624

return y

}

// Generate an array of 624 untempered numbers

function generateNumbers() {

for i from 0 to 623 {

int y := 32nd bit of (MT[i]) + last 31 bits of (MT[(i+1) mod 624])

MT[i] := MT[(i + 397) mod 624] xor (right shift by 1 bit (y))

if (y mod 2) == 1 { // y is odd

MT[i] := MT[i] xor (2567483615) // 0x9908b0df

}

SFMT

SFMT, the SIMD -oriented Fast Mersenne Twister, is a variant of Mersenne Twister, introduced in 2006^[5] , designed to be fast when it runs on 128-bit SIMD.

It is roughly twice as fast as Mersenne Twister.^[6]
It has a better equidistribution property of v-bit accuracy than MT but worse than WELL ("Well Equidistributed Long-period Linear").
It has quicker recovery from zero-excess initial state than MT, but slower than WELL.
It supports various periods from 2⁶⁰⁷ -1 to 2²¹⁶⁰⁹¹ -1.

Intel SSE2 and PowerPC AltiVec are supported by SFMT. It is also used for games with the Cell BE in the Playstation 3 .^[7]

References

External links

Implementations

这里有这种算法三个版本的C/C++实现代码：http://www.cppblog.com/Files/Chipset/random.zip

a better random generator // random.hpp // Copyright (C) 2008 Chipset // // This program is free software: you can redistribute it and/or modify // it under the terms of the GNU Affero General Public License as // published by the Free Software Foundation, either version 3 of the // License, or (at your option) any later version. // // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU Affero General Public License for more details. // // You should have received a copy of the GNU Affero General Public License // along with this program. If not, see < http://www.gnu.org/licenses/ >. // #ifndef RANDOM_HPP_ #define RANDOM_HPP_ #include < ctime > class random { public : explicit random(unsigned long s = 0 ) : seed(s) { if ( 0 == seed) seed = std::time( 0 ); randomize(); } void reset(unsigned long s = 0 ) { seed = s; if ( 0 == seed) seed = std::time( 0 ); randomize(); } unsigned long rand( ) { // returns a random integer in the range [0, -1UL) randomize(); return seed ; } double real() { // returns a random real number in the range [0.0, 1.0) randomize(); return double (seed) / -1UL; } private : unsigned long seed; void randomize() { seed = 1103515245UL * seed + 123456789UL ; } }; class rand_help { static random r; public : rand_help() {} void operator ()(unsigned long s) { r.reset(s); } unsigned long operator ()() const { return r.rand() ; } double operator ()( double ) { return r.real(); } }; random rand_help:: r; extern void srandom(unsigned long ns = 0 ) { rand_help()(ns); } // reset seed extern unsigned long irand() { return rand_help()(); } // negative numbers disallowed extern double drand() { return rand_help()( 1.0 ); } // for real numbers #endif // RANDOM_HPP_ // 以上随机数产生器产生的随机数比rand()产生的随机数更加随机(可以用数学方法检验)， //范围更大(一目了然)，速度更快(测试一下便知，稍加修改我还可以让它再快一些，如果有必要)。 /* 假设随机数均匀分布为理想分布, 粗略估计随机性 */ #include < iostream > #include < vector > #include < iomanip > #include " random.hpp " int main() { srand(time( 0 )); // SZ分别取值2^3, 2^4, , 2^15进行测试 const size_t SZ = 1 << 15 ; std::vector < unsigned > v1(SZ), v3(SZ); std::vector < unsigned > v2(SZ), v4(SZ); for (size_t i = 0 ; i < SZ; ++ i) { ++ v1[rand() % SZ]; // 对应元素计数 ,理论上v1[0] ~ v1[SZ - 1]之间每个都应该是1 ++ v2[irand() % SZ]; } for (size_t i = 0 ; i < SZ; ++ i) { ++ v3[v1[i]]; // 统计频度 ++ v4[v2[i]]; } // 假设［0, SSZ)之间不存在间断点, (即使有间断点也无所谓，我们只做粗略模糊统计，因为没有必要那么精确) // 最理想的显示结果应该是0: 0, 1: SZ, 2: 0, 3: 0, 4: 0, other: 0 // 0:表示间断，1:表示均匀，2:, 3:, 4:, other: 都表示不同程度的重复 const size_t SSZ = 5 ; std::cout.fill( ' ' ); for (size_t i = 0 ; i < SSZ; ++ i) std::cout << i << " : " << std::setw(SSZ) << v3[i] << " , " ; std::cout << " other: " << std::setw(SSZ) << v3.size() - v3[ 0 ] - v3[ 1 ] - v3[ 2 ] - v3[ 3 ] - v3[ 4 ] << ' /n ' ; for (size_t i = 0 ; i < SSZ; ++ i) std::cout << i << " : " << std::setw(SSZ) << v4[i] << " , " ; std::cout << " other: " << std::setw(SSZ) << v4.size() - v4[ 0 ] - v4[ 1 ] - v4[ 2 ] - v4[ 3 ] - v4[ 4 ] << ' /n ' ; system( " pause " ); } // 做速度测试 #include < iostream > #include < ctime > #include < cstdlib > #include < windows.h > #include " random.hpp " int main() { const size_t SZ = 1 << 27 ; std::cout << " generating random numbers in progress /n " ; std::cout << SZ << " random numbers generated./n " ; std::cout << " random used: " ; Sleep( 1000 ); std::clock_t time = clock(); for (size_t i = 0 ; i < SZ; ++ i) irand(); std::cout << clock() - time << " ms, " ; std::cout << " rand() used: " ; Sleep( 1000 ); std::srand(std::time( 0 )); std::clock_t t = clock(); for (size_t i = 0 ; i < SZ; ++ i) std::rand(); std::cout << clock() - t << " ms/n " ; std::system( " PAUSE " ); return 0 ; }

产生伪随机数常用的两种算法我们讲的随机数其实暗指伪随机数。提及随机数，不少朋友可能想到C语言的库函数rand()，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的一种实现版本如下： http://www.cppblog.com/Chipset/archive/2008/12/20/69918.html 如果需要高质量的伪随机数，内存充足(约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 一种实现版本如下： 1 // ********************************************************************** 2 // This is a slightly modified version of Equamen mersenne twister. 3 // 4 // Copyright (C) 2009 Chipset 5 // 6 // This program is free software: you can redistribute it and/or modify 7 // it under the terms of the GNU Affero General Public License as 8 // published by the Free Software Foundation, either version 3 of the 9 // License, or (at your option) any later version. 10 // 11 // but WITHOUT ANY WARRANTY; without even the implied warranty of 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 // GNU Affero General Public License for more details. 14 // 15 // You should have received a copy of the GNU Affero General Public License 16 // along with this program. If not, see < http://www.gnu.org/licenses/ >. 17 // ********************************************************************** 18 19 // Original Coyright (c) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura 20 // 21 // Functions for MT19937, with initialization improved 2002/2/10. 22 // Coded by Takuji Nishimura and Makoto Matsumoto. 23 // This is a faster version by taking Shawn Cokus's optimization, 24 // Matthe Bellew's simplification, Isaku Wada's real version. 25 // C++ version by Lyell Haynes (Equamen) 26 // 27 // Redistribution and use in source and binary forms, with or without 28 // modification, are permitted provided that the following conditions 29 // are met: 30 // 31 // 1. Redistributions of source code must retain the above copyright 32 // notice, this list of conditions and the following disclaimer. 33 // 34 // 2. Redistributions in binary form must reproduce the above copyright 35 // notice, this list of conditions and the following disclaimer in the 36 // documentation and/or other materials provided with the distribution. 37 // 38 // 3. The names of its contributors may not be used to endorse or promote 39 // products derived from this software without specific prior written 40 // permission. 41 // 42 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 43 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 44 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 45 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 46 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 47 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 48 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 49 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 50 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 51 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 52 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 53 // 54 55 #ifndef mtrandom_HPP_ 56 #define mtrandom_HPP_ 57 58 #include < stddef.h > 59 60 class mtrandom 61 { 62 public : 63 mtrandom() : left( 1 ) { init(); } 64 65 explicit mtrandom(size_t seed) : left( 1 ) { init(seed); } 66 67 mtrandom(size_t * init_key, int key_length) : left( 1 ) 68 { 69 int i = 1 , j = 0 ; 70 int k = N > key_length ? N : key_length; 71 init(); 72 for (; k; -- k) 73 { 74 state[i] = (state[i] ^ ((state[i - 1 ] ^ (state[i - 1 ] >> 30 )) * 1664525UL )) + init_key[j] + j; // non linear 75 state[i] &= 4294967295UL ; // for WORDSIZE > 32 machines 76 ++ i; 77 ++ j; 78 if (i >= N) 79 { 80 state[ 0 ] = state[N - 1 ]; 81 i = 1 ; 82 } 83 if (j >= key_length) 84 j = 0 ; 85 } 86 87 for (k = N - 1 ; k; -- k) 88 { 89 state[i] = (state[i] ^ ((state[i - 1 ] ^ (state[i - 1 ] >> 30 )) * 1566083941UL )) - i; // non linear 90 state[i] &= 4294967295UL ; // for WORDSIZE > 32 machines 91 ++ i; 92 if (i >= N) 93 { 94 state[ 0 ] = state[N - 1 ]; 95 i = 1 ; 96 } 97 } 98 99 state[ 0 ] = 2147483648UL ; // MSB is 1; assuring non-zero initial array 100 } 101 102 void reset(size_t rs) 103 { 104 init(rs); 105 next_state(); 106 } 107 108 size_t rand() 109 { 110 size_t y; 111 if ( 0 == -- left) 112 next_state(); 113 y = * next ++ ; 114 // Tempering 115 y ^= (y >> 11 ); 116 y ^= (y << 7 ) & 0x9d2c5680UL ; 117 y ^= (y << 15 ) & 0xefc60000UL ; 118 y ^= (y >> 18 ); 119 return y; 120 } 121 122 double real() { return ( double )rand() / - 1UL ; } 123 124 // generates a random number on [0,1) with 53-bit resolution 125 double res53() 126 { 127 size_t a = rand() >> 5 , b = rand() >> 6 ; 128 return (a * 67108864.0 + b) / 9007199254740992.0 ; 129 } 130 131 private : 132 void init(size_t seed = 19650218UL ) 133 { 134 state[ 0 ] = seed & 4294967295UL ; 135 for ( int j = 1 ; j < N; ++ j) 136 { 137 state[j] = ( 1812433253UL * (state[j - 1 ] ^ (state[j - 1 ] >> 30 )) + j); 138 // See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. 139 // In the previous versions, MSBs of the seed affect 140 // only MSBs of the array state[]. 141 // 2002/01/09 modified by Makoto Matsumoto 142 state[j] &= 4294967295UL ; // for >32 bit machines 143 } 144 } 145 146 void next_state() 147 { 148 size_t * p = state; 149 int i; 150 151 for (i = N - M + 1 ; -- i; ++ p) 152 * p = (p[M] ^ twist(p[ 0 ], p[ 1 ])); 153 154 for (i = M; -- i; ++ p) 155 * p = (p[M - N] ^ twist(p[ 0 ], p[ 1 ])); 156 * p = p[M - N] ^ twist(p[ 0 ], state[ 0 ]); 157 left = N; 158 next = state; 159 } 160 161 size_t mixbits(size_t u, size_t v) const 162 { 163 return (u & 2147483648UL ) \| (v & 2147483647UL ); 164 } 165 166 size_t twist(size_t u, size_t v) const 167 { 168 return ((mixbits(u, v) >> 1 ) ^ (v & 1UL ? 2567483615UL : 0UL )); 169 } 170 171 static const int N = 624 , M = 397 ; 172 size_t state[N]; 173 size_t left; 174 size_t * next; 175 }; 176 177 class mtrand_help 178 { 179 static mtrandom r; 180 public : 181 mtrand_help() {} 182 void operator ()(size_t s) { r.reset(s); } 183 size_t operator ()() const { return r.rand(); } 184 double operator ()( double ) { return r.real(); } 185 }; 186 mtrandom mtrand_help:: r; 187 188 extern void mtsrand(size_t s) { mtrand_help()(s); } 189 extern size_t mtirand() { return mtrand_help()(); } 190 extern double mtdrand() { return mtrand_help()( 1.0 ); } 191 192 #endif // mtrandom_HPP_ 193