很容易想到一个简单的做法
x = rand () % RANGE; /* Poor! return a random number in [0,RANGE) */
这种做法有三个问题:
1 rand()等概率返回[0, RAND_MAX]间这RAND_MAX+1个数中的任意一个,如果RAND_MAX+1不能被RANGE整除,那么这个[0, RANGE)分布就blatantly不均匀了。举个例子RAND_MAX=32767,目标区间为[0,10),32767%10=7,则(7,10)间数的取值概率小于[0,7]间数。
2 rand()的伪随机性。大多数随机数发生器产生随机数的低位bits都distressingly不随机(有一个short cyclical pattern或者bits间有依赖关系)。
3 RAND_MAX在<stdlib.h>中定义,其值最小为32767(2^15-1),最大为2147483647(2^31-1),假如RAND_MAX+1 < RANGE呢?
先讨论RANGE<=RAND_MAX+1的情况
方法1:
这种方法相当于取高位bits
(int)(rand() / (RAND_MAX + 1.0) * RANGE)
但是这种方法解决了问题2,却并没有解决问题1,因为如果RAND_MAX+1不能被RANGE整除,不管按照什么规则,RAND_MAX+1个数始终是无法均分放到RANGE个格子里的,输入数据量较大时这个bias就会比较明显。
方法2:
去除x=rand()所产生最大的(RAND_MAX+1)%RANGE个数,再(int)x/(int)((RANG_MAX+1.0)/RANGE)
int myrandom_range(int start,int end)//[start,end)
{
int N = end - start;
int bound = (RAND_MAX + 1) / N * N;
int r;
do
{
r = rand();
} while(r >= bound);
return start + int(r / (RAND_MAX + 1.0) * N);
}
这种方法解决了问题1&2,但是这个做法在RAND_MAX%RANGE (< RANGE)比较大的时候,可能会很浪费时间(极端情况是接近1/2概率需要重试),一般情况下还是够用的。
Let p = (RAND_MAX % RANGE) / (RAND_MAX + 1.0). This is the probability that any given call to rand() will require a retry. Note that this value is maximized when RANGE*2 = RAND_MAX+3, and which will yield a value of p roughly equal to 1/2.
- The average number of times that rand() will be called is 1/(1-p) (with a worst case of about 2).
- The probability that at least n calls to rand() will be required beyond the first one is pn.
So for most values of p which will be much less than 1/32 say, performance should not be a concern at all.
再讨论RANGE>RAND_MAX+1的情况
现在我们来考虑一个相对简单的问题,RANGE能被RAND_MAX+1整除的情况,例如数据范围[0,2^32)。
方法1加强版:
前面提到,RAND_MAX其值最小为(2^15-1),最大为(2^31-1),按最坏的情况算至少需要调用三次rand(),才能获得32位随机数,我们采取方法1中取高位的做法,三次rand()分别取低15位中的前10,11,11位。
inline unsigned __int32 rand32()
{
return ((rand()&0x00007FE0)>>5) + ((rand()&0x00007FF0)<<6) + ((rand()&0x00007FF0)<<17);
}
方法3:
可以考虑分段抽样,分成[n/(RNAD_MAX+1)]段,先等概率得到段再得到每段内的某个元素,这样分段也类似地有一个尾数问题,不是每次都刚好分到整数段,一定或多或少有一个余数段,这部分的值如何选取?
选到余数段的数据拿出来选取,先进行一次选到余数段概率的事件发生,然后进行单独选取:r = N % (RAND_MAX+1); //余数if ( happened( (double)r/N ) )//选到余数段的概率result = N-r+myrandom(r); // myrandom可以用情况1中的代码实现elseresult = rand()+myrandom(N/(RAND_MAX+1))*(RAND_MAX+1); // 如果选不到余数段再进行分段选取const double MinProb=1.0/(RAND_MAX+1);
bool happened(double probability)//probability 0~1
{
if(probability<=0)
{
return false;
}
if(probability
{
return rand()==0&&happened(probability*(RAND_MAX+1));
}
if(rand()<=probability*(RAND_MAX+1))
{
return true;
}
return false;
}
long myrandom(long n)//产生0~n-1之间的等概率随机数
{
t=0;
if(n<=RAND_MAX)
{
long R=RAND_MAX-(RAND_MAX+1)%n;//尾数
t = rand();
while ( t > r )
{
t = rand();
}
return t % n;
}
else
{
long r = n%(RAND_MAX+1);//余数
if( happened( (double)r/n ) )//取到余数的概率
{
return n-r+myrandom(r);
}
else
{
return rand()+myrandom(n/(RAND_MAX+1))*(RAND_MAX+1);
}
}
}
最后,使用之前,别忘了设随机种子srand((unsigned int)time(0))哦。
如果你懒得自己写,干脆用,
produces integers in the range [0, 232-1].
参考: