线性同余随机发生器

来源：互联网发布：长江大学有网络教育吗编辑：程序博客网时间：2024/06/07 11:24

一. LCG(linear congruential generator)线性同余发生器

伪随机数生成器

LCG 算法数学上基于公式：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有很好的应用, 例如内存很紧张的嵌入式中, 电子游戏控制台用的小整数, 使用高位可以胜任.

二. MCG

当LCG的c = 0，它被称为multiplier congruence generator(MCG). 这种情况下不可能达到满周期，即周期一定小于M。且种子不能为0（也不能是M、2M、3M...），否则得到的序列将是全zero。

在MCG中，若a是M的Primitive Root，则此MCG的周期等于T-1。一个具体的例子，当M=2^31 − 1，则M具有以下Primitive Root: 16807, 397204094, 764261123, 630360016. 就我目前看过的代码中, 选16807的居多.

因为 M = 2^31 − 1 = 2147483647 = 127773 * 16807 + 2836

我们可以用127773把seed分成高位和低位两部分.

seed = hi* 127773 + lo ;

所以 seed * 16807 = hi* 127773 * 16807 + lo * 16807 = hi * (M - 2836) + lo * 16807

上式最右端对M取余后得到 lo * 16807 - hi * 2836

但是这个计算结果有可能是负数, 此时让它加上M. 就回到[0,M)的范围内了.

但是这个技巧其实并不高明. 因为把seed拆成hi和lo两部分的时候, 需要计算除法. 而在现代CPU上, 除法的代价是乘法的10-30倍, 很可怕.

注意，multiplier congruence generator（c=0）生成的是[1，M - 1]上的均匀分布。所以一般来说会把它减去1后再返回，这样得到的就是[0, M-2]上的均匀分布. 别看这是个小东西, leveldb的代码不就犯了这个BUG吗? 假如你拿它做随机抽样, 结果数组中第一个元素永远也抽不中.

转自: http://www.sunchangming.com/blog/post/4159.html

三. 应用

3.1 C语言中伪随机数生成方法: rand(); srand(time(0))的解析

C语言中伪随机数生成算法实际上就是使用了LCG: X(n + 1) = (A * X(n) + C) % M

其中A, C,M都是常数(一般会取质数). 当C = 0时, 叫做乘同余法. 引出一个概念叫seed, 它会被作为X0被代入上式中, 然后每次调用rand()函数都会用上一次产生的随机值来生成新的随机值. 可以看出实际上用rand()函数生成的是一个递推的序列, 一切值都来源于最初的 seed. 所以当初始的seed取一样的时候, 得到的序列都相同.

C语言里面有RAND_MAX这样一个宏, 定义了rand()所能得到的随机值的范围. 在C里可以看到RAND_MAX被定义成0x7fff, 也就是32767. rand()函数里递推式中M的值就是32767.

线性同余法生成的是伪随机数, 粗略符合均匀分布. 根据中心极限定理, 任何分布的噪声, 通过反复相加, 就可以成为高斯噪声.

3.1.1 rand()函数

rand()函数返回0 --- RAND_MAX之间均匀分布的伪随机整数. RAND_MAX必须至少为32767. rand()函数不接受参数, 默认以1为种子(就是X0). 随机数生成器总是以相同的种子开始, 所以形成的伪随机数列也相同, 失去了随机意义. (但这样便于程序调试)

3.1.2 srand()函数

srand()设置不同的种子(设置X0), 因为time()返回的时间值每时每刻都不相同, 所以可以作为随机种子使用, 也就有了srand(time());

转自: http://www.cnblogs.com/xkfz007/archive/2012/03/27/2420154.html

http://www.sunchangming.com/blog/post/4159.html

3.2. C语言实现16807随机数产生器

int z = 123;                     // z 为种子，需提前设置(GetTickCount())  int schrage_next(){    const int a = 16807;// 16807 法    const int c = 0;    const int m = 2147483647;   // MAX_INT    const int q = m / a;        // q = m / a;    const int r = m % a;        // r = m % a;    int _z = a * (z % q) - r * (int)(z / q) + c; // 计算 mod    if(_z < 0)        _z += m;             // 将结果调整到 0 ~ m    return z = _z;}

X(n+1) = (A * X(n) + C) % M

其中:

A = 16807 (就是7^5)

M = 2147483647 (就是2^31 - 1)

C = 0

巨大的疑问: X(n+1) = (A * X(n) + C) % M是怎样变成int _z = a * (z % q) - r * (int)(z / q) + c;的, 还是有一点的疑问.

转自:

http://zhidao.baidu.com/link?url=6nXiMmBihBuFip5E7lG-whzTxSKPRiAgZtthXUZ8lCyLco8rFCC1uWX3UpQ_nvzaA_WOKhrSzsRymi8BwlJNda

3.3 std::hash_map中的hash_compare

template<class _Kty,class _Ty,class _Tr = hash_compare<_Kty, less<_Kty> >,class _Alloc = allocator<pair<const _Kty, _Ty> > >class hash_map: public _Hash<_Hmap_traits<_Kty, _Ty, _Tr, _Alloc, false> >{......};/////////////////////////////////////////////////////////////////////////////////// TEMPLATE CLASS hash_comparetemplate<class _Kty,class _Pr = less<_Kty> >class hash_compare{// traits class for hash containerspublic:......size_t operator()(const _Kty& _Keyval) const{// hash _Keyval to size_t value by pseudorandomizing transformlong _Quot = (long)(hash_value(_Keyval) & LONG_MAX);ldiv_t _Qrem = _CSTD ldiv(_Quot, 127773);                // 这里明显也是使用了MCG_Qrem.rem = 16807 * _Qrem.rem - 2836 * _Qrem.quot;if (_Qrem.rem < 0)_Qrem.rem += LONG_MAX;return ((size_t)_Qrem.rem);}bool operator()(const _Kty& _Keyval1, const _Kty& _Keyval2) const{// test if _Keyval1 ordered before _Keyval2return (comp(_Keyval1, _Keyval2));}protected:_Pr comp;// the comparator object};

四. 一些评价

线性同余发生器:

X(n+1) = (A * X(n) + C) % M

这方法有4个参数:

M: 模数 m > 0, 一般去 2^31 - 1

A: 乘数 0 <= A < M, 例如 A = 7^5 = 16807

C: 增量 0 <= C < M

X0: 种子 0 <= X0 < M

随机数序列{Xn}通过下列迭代方程得到: X(n+1) = (A * X(n) + C) % M

如果A, C, M和X0都是整数, 得到的序列都是整数, 并且 0 <= Xn < M

A, C, M, X0的选择对于算法非常关键.

例如:

M = 32, A = 7, C = 0, X0 = 1, 产生的序列为: {7, 17, 23, 1, 7, 17, 23, ......}, 周期是4;

如果A = 5, 序列为: {1, 5, 25, 29, 17, 21, 9, 13, 1, 5, ......}, 周期是8;

一般M选取为计算机所能表示的最大非负整数;

评价随机数产生器的准则

1). 这个函数应该是一个完整的周期的产生函数, 这个函数的在重复之前产生出0到M之间的所有书;

2). 产生的序列应该看起来是随机的;

3). 这个函数应该用32位算数高效实现;

在线性同余发生器中

在A的超过20亿的可选数字钟, 只有不多的几个可以通过所有3个准则, 其中一个是 A = 7^5 = 16807

优点:

如果对乘数和模数进行适当的选择, 得到的随机数序列和从集合{1, 2, ..., M}中随机抽取的数字序列式无法分辨的;

速度快, 每位只需要很少的操作;

缺点:

线性同余伪随机数具有可预测性, 随机仅在于X0的选取;

敌对方一旦知道采用的是线性同余算法, 并且参数已知(只要知道序列中3个数就可以求解出A, C, M), 那么所有的数字就都知道了.

参考:

http://www.cnblogs.com/xkfz007/archive/2012/03/27/2420154.html

http://www.sunchangming.com/blog/post/4159.html

http://wenku.baidu.com/link?url=LCHN3NqeHXKN7VbQnrslV7E4O4U2_qy9-bT-86HooAmMzBLdVvEkEXnqaawAkAtfP9n8sYOI6byk_Ytdnwis-POfUFCNWGqM4HryBe7sw9S

http://zhidao.baidu.com/link?url=1jBpjRyS5cGycUXO1H4eUnPHxtVv5v3CgiNLe0VQIxmg_CG-qRxwVRQtqSowzAXg8S0hU3AVxI_p5ubZd41yvq

参考网址:

http://www.cnblogs.com/xkfz007/archive/2012/03/27/2420154.html

http://www.sunchangming.com/blog/post/4159.html

http://wenku.baidu.com/link?url=LCHN3NqeHXKN7VbQnrslV7E4O4U2_qy9-bT-86HooAmMzBLdVvEkEXnqaawAkAtfP9n8sYOI6byk_Ytdnwis-POfUFCNWGqM4HryBe7sw9S

0 0