上篇结束的时候,我们提出了这个问题:
对于大于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位的随机数,我以为是不可取的。至于其原因,就是另一篇文章了。
“至于其原因,就是另一篇文章了”,随便点一下呗
线性同余的伪随机性不好,取低几位也是一样的。这就是为什么K&R里面是取高几位了。具体可以参考TAOCP vol2的相关小节。