注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Kandy

长风破浪会有时,直挂云帆济沧海!

 
 
 

日志

 
 

【转】产生伪随机数常用的两种算法  

2010-09-08 17:55:50|  分类: 默认分类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
我们讲的随机数其实暗指伪随机数。提及随机数,不少朋友可能想到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[M] ^ twist(p[0], p[1]));
153 
154         for(i = M; --i; ++p)
155             *= (p[M - N] ^ twist(p[0], p[1]));
156         *= 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 
  评论这张
 
阅读(878)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017