May | 2010 | i, Claud

上篇结束的时候，我们提出了这个问题：

对于大于RAND_MAX的正整数n，如何通过rand()得到一个0到n之间的随机数？

由前面的讨论可知，只要构造出范围比n大的随机数生成器，这个问题就解决了。

如何用小范围生成器构造大范围生成器？第一直觉往往是：

result = rand()*rand();

请注意：这种算法下，大于RAND_MAX的素数及其倍数产生的概率为0。这是由取值的离散性导致的。然而，即便是连续的随机变量，两个均匀分布之积通常也不是均匀分布。

产生这种错误直觉的原因在于：我们会优先考虑形象直观的“范围”性质，而把不直观的“随机”性抛弃或者简化。这个原因也适用于上一次的讨论。

相乘的错误很容易被发现。而第二直觉往往是：

result = rand() + rand();

这次，所有的结果都可以取到了。但最简单的反例：结果为1的概率是结果为0的概率的两倍。事实上，两个均匀分布之和，也不是均匀分布。一个很漂亮的结果是：很多个均匀分布之和（正规化以后），会趋于一个正态分布，这是中心极限定律的最简情形。

下面这个算法的错误就不容易被发现了，也因此成为大部分程序员使用的方法：

result = int ( rand()/float(RAND_MAX) * n ) ;

它的想法是：先生成[0,1]区间上的随机浮点数，然后乘以欲得的范围n。

请思考这样一个事实：rand()一共可能生成（RAND_MAX+1）个整数，将它们除以常数、乘以常数并取整，得到的依然是（RAND_MAX+1）个整数，其分布范围比原来要广。因此，一定有一些整数得不到。

还有一种算法，错误更难被发现（我甚至见到了基于它的技术文献）。在实践中，我们往往是一次生成许多个随机数（只生成一个，不具有统计意义）。这种算法是：

result[0] = rand();
for( i = 1 ; i < m ; ++i )
{
    result[i] = ( rand() + result[i-1] )%n ;
}

即循环相加得到一组“随机”数。它的问题何在？

我曾试图揣摩设计这种算法的人是如何思考的：首先，确实可以取到0到n的所有数；其次，结果大抵上是在各个区间段均匀分布的；最后，通过相加把随机进一步“搞混乱”，得到的应该是很随机的结果。他们的思维被局限在了“区间段”上面。

依旧以一例说明它是错的：第一个数必然落在0到RAND_MAX之间，因此，m个随机结果都落在RAND_MAX+1到n之间的情况出现的概率为0。而若需要的是m个独立随机变量（这也是用于统计的基本前提），它们都落在RAND_MAX+1到n之间的概率是大于0的。

当然了，这个例子已经超出了我们讨论的范围——生成“一个”随机数。

那么，到底该如何做呢？一个容易理解的方法是：按数位逐位生成。

不妨设n = 2^m -1，那么0到n均匀分布的随机变量，其二进制表示中的任何一位是0或是1的概率都为1/2，即满足0到1的均匀分布。因此，可以通过rand()用前一篇的方法得到每一位，然后拼起来即可。

这种算法的缺点在于，时间复杂度达到了log(n)。一种改善的思路是：一块一块的生成，然后拼出比n大的随机数。

最后，让我们回到现实。在标准C中，RAND_MAX比我们想象的要小很多，它的值是0×7fff = 32767，也就是signed short的最大取值。因此，如果使用它，当我们需要的随机范围大于3万，就会遇到本文所述扩大随机范围的问题。当然，你也可以使用其他的伪随机数生成器。

稍微思考一下就会发现一个问题：C语言原生地支持32位整型数，为什么这里的RAND_MAX是2^15-1而不是2^31-1呢？

另一个常识之误是，我们总认为标准C里的rand()是使用的线性同余算法。事实上，我们拿K&R的书，翻到第46页，可以看到以下代码：

unsigned long int next = 1;
int rand( void )
{
    next = next * 1103515245 + 12345;
    return ( unsigned int )( next / 65536 ) % 32768;
}

请注意next/65536，即取高16位。这是对线性同余法的一种改进。

我个人的认为是，这种改进能消除rand()%n用法的低位短周期性（但还没给出有效的证明）。至于其他的考虑，我没有去作考证。

正因为这个改进，rand()的范围被限制在了16位而非32位。

也有人提到使用next = ( next * 1103515245 + 12345 ) & 0×7fffffff来生成32位的随机数，我以为是不可取的。至于其原因，就是另一篇文章了。

给定一个伪随机数生成器rand()，假定它已达到一个理想状态，即输出的值满足0到RAND_MAX上的均匀分布。

对于小于RAND_MAX的正整数n，如何通过rand()得到一个0到n之间的随机数？

常规方法是这样的：

result = rand()%(n+1);

事实上，得到的结果并非均匀分布。例如：当n=RAND_MAX*2/3时，result为0 — n/2中某个数的概率是其为n/2+1 — n中某个数的概率的两倍。（这是差距最大的一种情况）

应该使用这种方法：

do{

    result = rand();

}while( result &gt; n );

但这种方法有个缺点：rand()的运算次数满足几何分布，故期望为：(RAND_MAX+1)/(n+1)。当n较小时，效率太低。

一种折中方案是：

do{

    result = rand();

}while( result &gt;= (RAND_MAX+1)/(n+1)*(n+1) );

result = result%(n+1);

这种算法具有较好的效率，唯一的缺点是：在实际使用时，需考虑RAND_MAX+1是否会造成溢出。

我们的下一个问题是：

对于大于RAND_MAX的正整数n，如何通过rand()得到一个0到n之间的随机数？

i, Claud