梅森旋转随机数生成实例

来源：互联网发布：js div显示隐藏编辑：程序博客网时间：2024/04/30 11:02

梅森旋转演算法（Mersenne twister）是一个伪随机数发生算法。由松本真和西村拓士^[1]在1997年开发，基于有限二进制字段上的矩阵线性递归 $F_{2}$ 。可以快速产生高质量的伪随机数，修正了古典随机数发生算法的很多缺陷。

Mersenne Twister这个名字来自周期长度取自梅森素数的这样一个事实。这个算法通常使用两个相近的变体，不同之处在于使用了不同的梅森素数。一个更新的和更常用的是MT19937, 32位字长。还有一个变种是64位版的MT19937-64。对于一个k位的长度，Mersenne Twister会在 $[0,2^k-1]$ 的区间之间生成离散型均匀分布的随机数。

优点：
最为广泛使用Mersenne Twister的一种变体是MT19937，可以产生32位整数序列。具有以下的优点:
周期非常长，达到2¹⁹⁹³⁷−1。尽管如此长的周期并不必然意味着高质量的伪随机数，但短周期（比如许多旧版本软件包提供的2³²）确实会带来许多问题。
在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[N];     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 linearstate[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 linearstate[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++;// Temperingy ^= (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 resolutiondouble 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 Matsumotostate[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[N];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.0cout << 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);  }

0 0