libc随机数生成器有缺陷吗？

考虑一种算法来测试在特定次数的尝试之后从一组N个唯一数字中挑选某个数字的概率（例如，N = 2，轮盘中的概率是什么（没有0），它需要X尝试黑赢？）

对此的正确分布是pow（1-1 / N，X-1）*（1 / N）。

但是，当我使用下面的代码测试时，X = 31处始终存在深沟，独立于N，并且独立于种子。

这是一个内在的缺陷，由于PRNG的实施细节在使用中无法防止，这是一个真正的错误，还是我忽略了一些明显的东西？

// C #include  #include  #include  int array[101]; void main(){ int nsamples=10000000; double breakVal,diffVal; int i,cnt; // seed, but doesn't change anything struct tms time; srandom(times(&time)); // sample for(i=0;i<nsamples;i++){ cnt=1; do{ if((random()%36)==0) // break if 0 is chosen break; cnt++; }while(cnt<100); array[cnt]++; } // show distribution for(i=1;i<100;i++){ breakVal=array[i]/(double)nsamples; // normalize diffVal=breakVal-pow(1-1/36.,i-1)*1/36.; // difference to expected value printf("%d %.12g %.12g\n",i,breakVal,diffVal); } } 

使用libc6软件包2.15-0ubuntu20和Intel Core i5-2500 SandyBridge测试了最新的Xubuntu 12.10，但几年前我在一台较旧的Ubuntu机器上发现了这一点。

我也在Windows 7上使用Unity3D / Mono进行了测试（虽然不知道哪个Mono版本），这里使用System.Random时X = 55的沟渠，而Unity的内置Unity.Random没有可见的沟渠（至少没有）对于X <100）。

分布：

差异：

Answer 1

这是因为glibc的random()函数不够随机。 根据这个页面 ，对于random()返回的随机数，我们有：

o i = (o i-3 + o i-31 ) % 2^31

要么：

o i = (o i-3 + o i-31 + 1) % 2^31 。

现在取x i = o i % 36 ，并假设上面的第一个等式是使用的那个（这种情况发生时每个数字的概率为50％）。 现在，如果x i-31 =0且x i-3 !=0 ，则x i =0的机会小于1/36。 这是因为o i-31 + o i-3 50％时间将小于2 ^ 31，当发生这种情况时，

x i = o i % 36 = (o i-3 + o i-31 ) % 36 = o i-3 % 36 = x i-3 ，

这是非零的。 这会导致您在0样本后看到31个样本的沟渠。

Answer 2

在该实验中测量的是伯努利实验的成功试验之间的间隔，其中成功被定义为对于某些k （在OP中为36 random() mod k == 0 。 不幸的是，由于random()的实施意味着伯努利试验在统计上不是独立的，因此可能会受到损害。

我们将为`random（）’的i th输出写rnd i ，我们注意到：

rnd i = rnd i-31 + rnd i-3 ，概率为0.75

rnd i = rnd i-31 + rnd i-3 + 1 ，概率为0.25

（见下面的校样大纲。）

让我们假设rnd i-31 mod k == 0我们现在正在看rnd i 。 然后必须是rnd i-3 mod k ≠ 0 ，因为否则我们将周期计为长度k-3 。

但是（大部分时间） (mod k): rnd i = rnd i-31 + rnd i-3 = rnd i-3 ≠ 0 。

因此，目前的试验在统计学上并不依赖于之前的试验，成功后的第31次试验成功的可能性要小于无偏见的伯努利试验系列。

使用线性同余生成器的常用建议（实际上并不适用于random()算法）是使用高阶位而不是低阶位，因为高阶位是“更随机”（也就是说，与连续值的相关性较小）。 但是在这种情况下也不会起作用，因为上述标识对于函数high log k bits同样适用于函数mod k == low log k bits 。

实际上，我们可能期望线性同余生成器更好地工作，特别是如果我们使用输出的高阶位，因为尽管LCG在蒙特卡罗模拟中不是特别好，但它不会受到线性反馈的影响。 random() 。

---

random算法，默认情况下：

让state成为unsigned longs的向量。 使用种子，一些固定值和混合算法初始化state 0 ...state 30 。 为简单起见，我们可以认为状态向量是无限的，尽管只使用了最后的31个值，所以它实际上是作为环形缓冲区实现的。

要生成rnd i : (Note: ⊕ is addition mod 2 ³² .)

state i = state i-31 ⊕ state i-3

rnd i = (state i - (state i mod 2)) / 2

现在，请注意：

(i + j) mod 2 = i mod 2 + j mod 2 if i mod 2 == 0或j mod 2 == 0

(i + j) mod 2 = i mod 2 + j mod 2 - 2如果i mod 2 == 1且j mod 2 == 1

如果i和j均匀分布，则第一种情况将发生在75％的时间，第二种情况发生在25％。

因此，通过生成公式中的替换：

rnd i = (state i-31 ⊕ state i-3 - ((state i-31 + state i-3 ) mod 2)) / 2

= ((state i-31 - (state i-31 mod 2)) ⊕ (state i-3 - (state i-3 mod 2))) / 2或

= ((state i-31 - (state i-31 mod 2)) ⊕ (state i-3 - (state i-3 mod 2)) + 2) / 2

这两种情况可以进一步简化为：

rnd i = rnd i-31 ⊕ rnd i-3

rnd _i = rnd _i- 31⊕rnd _i-3 + 1

如上所述，第一种情况发生在75％的时间，假设rnd _i-31和rnd _i-3独立地从均匀分布中抽出（它们不是，但它是合理的第一近似）。

Answer 3

正如其他人所指出的， random()不是随机的。

在这种情况下，使用较高位而不是较低位无效。 根据手册（ man 3 rand ）， rand() 旧实现在低位有问题。 这就是为什么建议使用random() 。 但是， rand()的当前实现使用与random()相同的生成器。

我尝试了推荐正确使用旧rand() ：

 if ((int)(rand()/(RAND_MAX+1.0)*36)==0) 

……在X = 31处得到了同样的深沟

如果我将rand()的数字与另一个序列混合在一起，我会摆脱沟渠：

 unsigned x=0; //... x = (179*x + 79) % 997; if(((rand()+x)%36)==0) 

我使用的是旧的线性同余发生器 。 我从素数表中随机选择了79,179和997。 这应该产生长度为997的重复序列。

也就是说，这个技巧可能会引入一些非随机性，一些足迹…由此产生的混合序列肯定会失败其他统计测试。 x在连续迭代中从不采用相同的值。 实际上，重复每个值需要997次迭代。

”[..]不应使用随机选择的方法生成随机数。 应该使用一些理论。“（DEKnuth，”计算机程序设计的艺术“，第2卷）

对于模拟，如果您想确定，请使用Mersenne Twister