回 帖 发 新 帖 刷新版面

主题:产生伪随机数两种常用算法

我们讲的随机数其实暗指伪随机数。不少朋友可能想到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++实现版本如下:
//************************************************************************
//  Copyright (C) 2008 - 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/>.
//************************************************************************
#ifndef LCRANDOM_HPP_
#define LCRANDOM_HPP_
#include <ctime>

class lcrandom
{
public:
  explicit lcrandom(size_t s = 0) : seed(s)
  {
    if (0 == seed) seed = std::time(0);
    randomize();
  }

  void reset(size_t s)
  {
    seed = s;
    if (0 == seed) seed = std::time(0);
    randomize();
  }

  size_t 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:
  size_t seed;
  void randomize() { seed = 1103515245UL * seed + 12345UL; }
};

class lcrand_help
{
  static lcrandom r;
public:
  lcrand_help() {}
  void operator()(size_t s) { r.reset(s); }
  size_t operator()() const { return r.rand(); }
  double operator()(double) { return r.real(); }
};
lcrandom lcrand_help:: r;

extern void lcsrand(size_t s) { lcrand_help()(s); }
extern size_t lcirand() { return lcrand_help()(); }
extern double lcdrand() { return lcrand_help()(1.0); }

#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算法原理编写的。

回复列表 (共10个回复)

沙发

我最近也写了一个版本(修改别人的),如果您发现任何bug请告诉我,谢谢!
//************************************************************************
//  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/>.
//************************************************************************

#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_

板凳

以上两个版本优点很明显,首先随机性好,而且速度快,如果您比较懒,还不必初始化就能用。不过以上两个实现都不是线程安全的,而且都不能用于加密,不过第二个实现稍加修改后可以用于加密的。

3 楼

这东西很强大...拜一个

4 楼

记得还有一种人字型算法(不记得是不是这么说的),它的非周期轨道点的分布密度函数是人字映射与线性同余法结合,可产生统计性质优良的均匀随机数。

for(int i=1;i<=n;i++)//线性同余法生成随机数
{
  f[i]=fmod((a*f[i-1]),m);
    if(f[i]<=m/2)//与人字映射结合生成随机数
      {
         f[i]=2*f[i];
      }
    else
       {
         f[i]=2*(m-f[i])+1;
       }
}


随机数产生的算法其实蛮多的,好像还有一种写法可以生成近似正态分布的随机数,不过我没掌握。

5 楼

[quote]记得还有一种人字型算法(不记得是不是这么说的),它的非周期轨道点的分布密度函数是人字映射与线性同余法结合,可产生统计性质优良的均匀随机数。

for(int i=1;i<=n;i++)//线性同余法生成随机数
{
  f[i]=fmod((a*f[i-1]),m);
    if(f[i]<=m/2)//与人字映射结合生成随机数
      {
         f[i]=2*f[i];
      }
    else
       {
         f[i]=2*(m-f[i])+1;
       }
}


随机数产生的算法其实蛮多的,好像还有一种写法可以生成近似正态分布的随机数,不过我没掌握。

[/quote]

主要问题是随机性怎么样,以及速度怎么样。

6 楼

谢谢LZ,学习了

7 楼


感谢分享  回头细看 。。。。

--------------------------
我看帖 我骄傲。。。

8 楼

to Chipset
生成随机数的方法很多,好像最主要的是在随机种子的选定上。
不知道你的种子是基于什么选定的?

9 楼

8f:
对,种子选择是很重要的问题,对于传统古老的LCG(linear congruential generator)尤其重要。其实随机数的产生是个数学问题,本贴从数学方面进行讨论的,但也给出了两种数学方法的代码实现示例。

10 楼

原来基础的问题蕴藏着丰富的智慧。

我来回复

您尚未登录,请登录后再回复。点此登录或注册