Monthly Archives: May 2010

5月中旬新书观察

先看计算机方面。本周排行版的两大亮点是产品、架构。

四月份出版的《人人都是产品经理》自出版后成为最畅销的书,与同期出版的另一本《结网:互联网产品经理改变世界》(销量第二)遥想辉映。这是一个可喜的现象,我们可以大胆地预计,在接下来一年里同类图书市场将会更加繁荣。自2000年左右计算机图书出版市场化以来,先后经历了语言、技术、模式、架构的变迁,终于一步步坚定地向用户需求靠近,这也是行业逐渐成熟的一种标志。

Continue reading

伪随机数使用之误(二)

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

对于大于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之间的随机数?