C++随机数--——生成任意范围内等概率随机数“足够好”的做法
很容易想到一个简单的做法
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中的代码实现
else
result = 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<MinProb)
{
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))哦。
如果你懒得自己写，干脆用boost库，mt19937 produces integers in the range [0, 232-1].
参考:
