梅森旋转随机数生成实例
2015-05-06 14:35
405 查看
梅森旋转演算法(Mersenne twister)是一个伪随机数发生算法。由松本真和西村拓士[1]在1997年开发,基于有限二进制字段上的矩阵线性递归
。可以快速产生高质量的伪随机数,
修正了古典随机数发生算法的很多缺陷。
Mersenne Twister这个名字来自周期长度取自梅森素数的这样一个事实。这个算法通常使用两个相近的变体,不同之处在于使用了不同的梅森素数。一个更新的和更常用的是MT19937, 32位字长。 还有一个变种是64位版的MT19937-64。 对于一个k位的长度,Mersenne
Twister会在
的区间之间生成离散型均匀分布的随机数。
梅森旋转算法实现:
应用举例:
。可以快速产生高质量的伪随机数,
修正了古典随机数发生算法的很多缺陷。
Mersenne Twister这个名字来自周期长度取自梅森素数的这样一个事实。这个算法通常使用两个相近的变体,不同之处在于使用了不同的梅森素数。一个更新的和更常用的是MT19937, 32位字长。 还有一个变种是64位版的MT19937-64。 对于一个k位的长度,Mersenne
Twister会在
的区间之间生成离散型均匀分布的随机数。
优点:
最为广泛使用Mersenne Twister的一种变体是MT19937,可以产生32位整数序列。具有以下的优点:
周期非常长,达到219937−1。尽管如此长的周期并不必然意味着高质量的伪随机数,但短周期(比如许多旧版本软件包提供的232)确实会带来许多问题。
在1 ≤ k ≤ 623的维度之间都可以均等分布(参见定义).
除了在统计学意义上的不正确的随机数生成器以外,在所有伪随机数生成器法中是最快的(当时)
梅森旋转算法实现://************************************************************************ // This is a slightly modified version of Equamen mersenne twister. // // Copyright (C) 2009 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/>. //************************************************************************ // Original Coyright (c) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura // // Functions for MT19937, with initialization improved 2002/2/10. // Coded by Takuji Nishimura and Makoto Matsumoto. // This is a faster version by taking Shawn Cokus's optimization, // Matthe Bellew's simplification, Isaku Wada's real version. // C++ version by Lyell Haynes (Equamen) // // Redistribution and use in source and binary forms, with or without // modification, are permitted provided that the following conditions // are met: // // 1. Redistributions of source code must retain the above copyright // notice, this list of conditions and the following disclaimer. // // 2. Redistributions in binary form must reproduce the above copyright // notice, this list of conditions and the following disclaimer in the // documentation and/or other materials provided with the distribution. // // 3. The names of its contributors may not be used to endorse or promote // products derived from this software without specific prior written // permission. // // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. // #ifndef mtrandom_HPP_ #define mtrandom_HPP_ #include <stddef.h> class mtrandom { public: mtrandom() : left(1) { init(); } explicit mtrandom(size_t seed) : left(1) { init(seed); } mtrandom(size_t* init_key, int key_length) : left(1) { int i = 1, j = 0; int k = N > key_length ? N : key_length; init(); for(; k; --k) { state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525UL))+ init_key[j] + j; // non linear state[i] &= 4294967295UL; // for WORDSIZE > 32 machines ++i; ++j; if(i >= N) { state[0] = state[N - 1]; i = 1; } if(j >= key_length) j = 0; } for(k = N - 1; k; --k) { state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941UL)) - i; // non linear state[i] &= 4294967295UL; // for WORDSIZE > 32 machines ++i; if(i >= N) { state[0] = state[N - 1]; i = 1; } } state[0] = 2147483648UL; // MSB is 1; assuring non-zero initial array } void reset(size_t rs) { init(rs); next_state(); } size_t rand() { size_t y; if(0 == --left) next_state(); y = *next++; // Tempering y ^= (y >> 11); y ^= (y << 7) & 0x9d2c5680UL; y ^= (y << 15) & 0xefc60000UL; y ^= (y >> 18); return y; } double real() { return (double)rand() / -1UL; } // generates a random number on [0,1) with 53-bit resolution double res53() { size_t a = rand() >> 5, b = rand() >> 6; return (a * 67108864.0 + b) / 9007199254740992.0; } private: void init(size_t seed = 19650218UL) { state[0] = seed & 4294967295UL; for(int j = 1; j < N; ++j) { state[j] = (1812433253UL * (state[j - 1] ^ (state[j - 1] >> 30)) + j); // See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. // In the previous versions, MSBs of the seed affect // only MSBs of the array state[]. // 2002/01/09 modified by Makoto Matsumoto state[j] &= 4294967295UL; // for >32 bit machines } } void next_state() { size_t* p = state; int i; for(i = N - M + 1; --i; ++p) *p = (p[M] ^ twist(p[0], p[1])); for(i = M; --i; ++p) *p = (p[M - N] ^ twist(p[0], p[1])); *p = p[M - N] ^ twist(p[0], state[0]); left = N; next = state; } size_t mixbits(size_t u, size_t v) const { return (u & 2147483648UL) | (v & 2147483647UL); } size_t twist(size_t u, size_t v) const { return ((mixbits(u, v) >> 1) ^ (v & 1UL ? 2567483615UL : 0UL)); } static const int N = 624, M = 397; size_t state ; size_t left; size_t* next; }; class mtrand_help { static mtrandom r; public: mtrand_help() {} void operator()(size_t s) { r.reset(s); } size_t operator()() const { return r.rand(); } double operator()(double) { return r.real(); } }; mtrandom mtrand_help:: r; extern void mtsrand(size_t s) { mtrand_help()(s); } extern size_t mtirand() { return mtrand_help()(); } extern double mtdrand() { return mtrand_help()(1.0); } #endif // mtrandom_HPP_
应用举例:
#include <cstdio> #include <cstdlib> #include <iostream> #include <stddef.h> #include <ctime> using namespace std; class mtrandom { public: mtrandom() : left(1) { init(); } explicit mtrandom(size_t seed) : left(1) { init(seed); } mtrandom(size_t* init_key, int key_length) : left(1) { int i = 1, j = 0; int k = N > key_length ? N : key_length; init(); for(; k; --k) { state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525UL))+ init_key[j] + j; // non linear state[i] &= 4294967295UL; // for WORDSIZE > 32 machines ++i; ++j; if(i >= N) { state[0] = state[N - 1]; i = 1; } if(j >= key_length) j = 0; } for(k = N - 1; k; --k) { state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941UL)) - i; // non linear state[i] &= 4294967295UL; // for WORDSIZE > 32 machines ++i; if(i >= N) { state[0] = state[N - 1]; i = 1; } } state[0] = 2147483648UL; // MSB is 1; assuring non-zero initial array } void reset(size_t rs) { init(rs); next_state(); } size_t rand() { size_t y; if(0 == --left) next_state(); y = *next++; // Tempering y ^= (y >> 11); y ^= (y << 7) & 0x9d2c5680UL; y ^= (y << 15) & 0xefc60000UL; y ^= (y >> 18); return y; } double real() { return (double)rand() / -1UL; } // generates a random number on [0,1) with 53-bit resolution double res53() { size_t a = rand() >> 5, b = rand() >> 6; return (a * 67108864.0 + b) / 9007199254740992.0; } public: void init(size_t seed = 19650218UL) { state[0] = seed & 4294967295UL; for(int j = 1; j < N; ++j) { state[j] = (1812433253UL * (state[j - 1] ^ (state[j - 1] >> 30)) + j); // See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. // In the previous versions, MSBs of the seed affect // only MSBs of the array state[]. // 2002/01/09 modified by Makoto Matsumoto state[j] &= 4294967295UL; // for >32 bit machines } } void next_state() { size_t* p = state; int i; for(i = N - M + 1; --i; ++p) *p = (p[M] ^ twist(p[0], p[1])); for(i = M; --i; ++p) *p = (p[M - N] ^ twist(p[0], p[1])); *p = p[M - N] ^ twist(p[0], state[0]); left = N; next = state; } size_t mixbits(size_t u, size_t v) const { return (u & 2147483648UL) | (v & 2147483647UL); } size_t twist(size_t u, size_t v) const { return ((mixbits(u, v) >> 1) ^ (v & 1UL ? 2567483615UL : 0UL)); } static const int N = 624, M = 397; size_t state ; size_t left; size_t* next; }; class mtrand_help { public: static mtrandom r; public: mtrand_help() {} void operator()(size_t s) { r.reset(s); } size_t operator()() const { return r.rand(); } double operator()(double) { return r.real(); } }; mtrandom mtrand_help:: r; void mtsrand(size_t s) { mtrand_help()(s); } size_t mtirand() { return mtrand_help()(); } double mtdrand() { return mtrand_help()(1.0); } int Random(int low, int high); int Random2(double start, double end); int Random3(double start, double end); int main(int argc, char *argv[]) { mtsrand(time(NULL)); int t = mtdrand()*10;//0~1.0 cout << t << endl; for(int i=0;i<100;i++) { if(i%5 == 0) cout << endl; cout << Random(1,50)<<" "; } cout << endl; return 0; } //[low high] int Random(int low, int high) { return low + mtirand() % (high - low + 1); } //[start end) int Random2(double start, double end) { return start+(end-start)*rand()/(RAND_MAX + 1.0); }
相关文章推荐
- 【随机数生成算法系列】线性同余法和梅森旋转法
- jQuery - 综合实例 - 生成在两个边界间的随机数
- php生成N个不重复的随机数实例
- 数组应用实例(生成并打印随机数和统计随机数的分布)
- PHP生成随机数的方法实例分析
- js生成随机数的方法实例
- python3生成随机数实例
- C#生成随机数实例
- mysql 存储过程实例 (日期以小时递增 while loop循环嵌套 随机数生成)
- python3生成随机数实例
- 一起talk C栗子吧(第六回:C语言实例--生成随机数)
- PHP随机数生成代码与使用实例分析
- js生成随机数的方法实例总结
- 伪随机数生成——梅森旋转(Mersenne Twister/MT)算法笔记
- C#生成不重复随机数列表实例
- 生成随机数方法:java.util.Random.nextInt(int n)方法实例
- C++随机数生成实例讲解
- java生成两数之间随机数实例
- C#--实例选号器生成随机数