一、简介
随机数是编程中经常要用到的东西,但很遗憾的是在windows下vc和MinGW中的 RAND_MAX 都是32767,
也就是说调用系统的rand()函数只能产生范围在[0, 32767]之间的随机数。
那如何得到范围达到[0, 2147483647]的随机数呢?
二、最直观的方法
最直接的想法显然是调用rand()生成更大区间的随机数,比如两个rand()相加或者相乘,但是,
两个rand()相加时,越大的数出现的概率越高,因为2=2+0=1+1=0+2,而5=0+5=1+4=2+3=3+2=4+1=5+0;
两个rand()相乘时,是永远不可能得到更大区间范围内的质数的,比如两个[0, 4]范围内的随机数相乘的所有组合都不可能得到7。
综上,似乎要寻找一种更靠谱的方法。
三、基于按位随机的随机数产生方法
首先我们知道,由一个较大的区间产生较小区间的随机数是十分方便的,只要把多余的数剔除掉,然后等分区间就可以了。
比如,由[0, 17]之间的随机数发生器A产生[0, 2]之间的随机数发生器B,只需要将 A 做一个映射:
[0, 4] --> 0
[5, 9] --> 1
[10, 14] --> 2
[15, 17] --> 丢弃
这样0、1、2出现的概率都是相等的,对于A而言都是5/18,对于B而言都是1/3。
特别的,给定任意一个随机数发生器A(产生[0, m]之间的随机数,m>0),一定能用这个随机数发生器产生[0, 1]的随机数发生器B:
m=1时A即为[0, 1]随机数发生器,m>1时,
若m为奇数,则将A产生的数直接%2;若m为偶数,若A产生的数等于m则丢弃,其它情况直接%2。
这样就得到了[0, 1]随机数发生器B。
那么,无论我们需要的随机数范围有多大,只要按位随机出0或者1,那么这个数一定随机的。
比如需要[0, 2147483647]之间的随机数,那么只要按位随机出31个二进制位,那么就得到了[0, 2147483647]之间的随机数。
实现代码如下:
1 View Code 2 /*生成随机数*/ 3 int randrange(int min, int max) 4 { 5 int range = max - min + 1; 6 7 if (range < 1) 8 { 9 throw std::range_error("The upper limit is less than the lower limit!"); 10 return 0; 11 } 12 //系统rand()函数产生的随机数范围是0到RAND_MAX(32767) 13 //这里产生0到2147483647之间的随机数 14 std::bitset<31> randomset; 15 for (int i = 0; i != 31; ++i) 16 { 17 randomset[i] = (rand() % 2 == 1); 18 } 19 int randomnumber = static_cast<int>(randomset.to_ulong()); 20 21 return (randomnumber % range + min); 22 }
看上去这个思路应该是没有问题的,因为每一位都是随机的,那整个数一定是随机的,但是测试结果不尽如人意。
测试代码:
1 std::bitset<2147483648> randomSet; 2 3 int main(void) 4 { 5 unsigned int seed = (unsigned int)time(0) + rand(); 6 long long maxCount = 21474836480L; 7 srand(seed); 8 randomSet.reset(); 9 10 std::cout << randomSet.count() << std::endl; 11 std::cout << "seed = " << seed << std::endl; 12 13 for (long long i = 0; i != maxCount; ++i) 14 { 15 int randNum = randRange(0, 2147483647); 16 randomSet[randNum] = true; 17 } 18 19 std::cout << randomSet.count() << std::endl << std::endl; 20 21 system("pause"); 22 return 0; 23 }
测试结果是:
randomSet.count()等于131066
测试代码中for循环执行了21474836480次,如果随机算法 “足够好”,那么randomSet中每一位都会有10次机会被置成true。
然而有测试结果来看,randomSet中却只有131066位被置为ture,这与2147483648相差甚远,看来这个随机算法真是 “足够不好”。
因此我们可以“武断”的认为,基于rand()来按位随机是完全不可行的。
那为什么这种做法结果不对呢。。。。????
四、系统的rand()函数是怎么实现的
其实这种按位随机的做法本没问题,但这基于一个假设:用于产生每位0或1的rand()必须是绝对随机的。
我对绝对随机(假设在[0, 99]范围内)的理解是:
1、要保证在随机次数足够大的情况下,随机出来的数在[0, 99]区间内服从均匀分布。这意味着如果[0, 99]的数不是等概率取到就不是绝对随机的。
2、每次产生的随机数必须是相互独立的。这意味着{0, 1, 2, ……, 99}这种有规律的数列虽然满足第一个条件,但它也不是绝对随机的。
那我们系统随机函数rand()是否满足绝对随机呢?
经过google,并查看vs2010和GLibC2.16.0中对rand()函数的定义可以知道,rand()能满足第一个条件却不满足第二个条件,下面看源码:
vs2010的 rand.c 默认在 “C:\Program Files\Microsoft Visual Studio 10.0\VC\crt\src”下,其源码如下:
rand.c
1 void __cdecl srand ( 2 unsigned int seed 3 ) 4 { 5 _getptd()->_holdrand = (unsigned long)seed; 6 } 7 8 int __cdecl rand ( 9 void 10 ) 11 { 12 _ptiddata ptd = _getptd(); 13 /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/ 14 return( ((ptd->_holdrand = ptd->_holdrand * 214013L 15 + 2531011L) >> 16) & 0x7fff ); 16 }
GLibC2.16.0在这里下载,其rand()函数的一套源码如下:
/stdlib/rand.c
1 /* /stdlib/rand.c */ 2 int 3 rand () 4 { 5 return (int) __random (); 6 }
/stdlib/random.c
1 /* /stdlib/random.c */ 2 long int 3 __random () 4 { 5 int32_t retval; 6 7 __libc_lock_lock (lock); 8 9 (void) __random_r (&unsafe_state, &retval); 10 11 __libc_lock_unlock (lock); 12 13 return retval; 14 } 15 16 weak_alias (__random, random)
/stdlib/random_r.c
1 /* /stdlib/random_r.c */ 2 int 3 __random_r (buf, result) 4 struct random_data *buf; 5 int32_t *result; 6 { 7 int32_t *state; 8 9 if (buf == NULL || result == NULL) 10 goto fail; 11 12 state = buf->state; 13 14 if (buf->rand_type == TYPE_0) 15 { 16 int32_t val = state[0]; 17 /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/ 18 val = ((state[0] * 1103515245) + 12345) & 0x7fffffff; 19 state[0] = val; 20 *result = val; 21 } 22 else 23 { 24 int32_t *fptr = buf->fptr; 25 int32_t *rptr = buf->rptr; 26 int32_t *end_ptr = buf->end_ptr; 27 int32_t val; 28 29 val = *fptr += *rptr; 30 /* Chucking least random bit. */ 31 *result = (val >> 1) & 0x7fffffff; 32 ++fptr; 33 if (fptr >= end_ptr) 34 { 35 fptr = state; 36 ++rptr; 37 } 38 else 39 { 40 ++rptr; 41 if (rptr >= end_ptr) 42 rptr = state; 43 } 44 buf->fptr = fptr; 45 buf->rptr = rptr; 46 } 47 return 0; 48 49 fail: 50 __set_errno (EINVAL); 51 return -1; 52 } 53 54 weak_alias (__random_r, random_r)
注意代码中加了/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/标记的语句,可以看到库函数产生随机数用的是线性同余法。
所谓线性同余法是基于 X(n+1) = (a * X(n) + c) % m 这样的公式递推产生随机数,其中a、c、m是常数,随机数种子也就是用以确定X(0)的数。
可以看到,在vs2010中,a=214013, c=2531011, m=32768(&0x7fff),在GLibC中a=1103515245, c=12345, m=2147483648(&0x7fffffff)。
所以vs2010中产生的随机数范围是[0, 32767], 而GLibC中则可达到[0, 2147483647]。后面所说的rand()默认指vs2010中的rand()。
因为线性同余法的随机数是由递推公式产生的,虽然能够保证大样本下的近似均匀分布(后面会有测试),但显然不能满足连续随机数之间的相互独立。
一般而言,采用这种方法会用time(0)作为随机数种子,time(0)的返回值是从1970年1月1日00:00到当前时刻所经历的时间(单位秒)。
对于程序的运行而言,如果不要求程序开始运行的时刻(不受控制的),我们可以近似认为随机数种子是不可预知的,因此,即使我们知道a、c、m的值也
无法预测会产生什么样的随机数。实际中我们可以认为rand()产生的数是“随机数”正是基于以上的原因(不可预测性)。
然而如果我们知道a、c、m的值和X(n),那么之后递推产生的所有的数都是可预知的。
所以,由于前面的方法需要连续调用rand()函数来随机出每一个二进制位,这本质上和调用一次rand()函数没有区别,除了第一次rand()之外后面都是有规律并直接取决于前一次结果。
那每次rand()都重置随机数种子呢?同样也是不行的,在编程实现上,除了程序开始运行的时刻是不可控(近似认为)的之外,后面的都是完全确定的,所以每次重置随机数
种子是在编程实现上要么一直srand(),而time(0)的值在1s内是不变的,要么每次Sleep(rand())一段时间,而这个rand()的结果已经由上次的随机数种子确定了。
综上所述,若要利用rand()产生我们自己的随机数,那么产生每个数的时候rand()最多只能用一次。
遗憾的是,从小空间(比如[0, 32767])到大空间(比如[0, 2147483647])的映射不可能是 一 一 映射,所以用rand()产生[0, 2147483647]范围内的随机数至少
要连续调用两次甚至更多的rand()。基于前面的结论,调用系统函数rand()产生[0, 2147483647]范围内均匀分布的随机数是不可能的。
再次遗憾的是,在MinGW(可以近似认为是windows下的GCC)下RAND_MAX仍然是32767,好吧,只要把GLibC中的rand()代码手动拿来用了。
五、GLibC中rand()函数的简单实现
下面是一个简单的版本,忽略了临界区访问控制等代码(如果不需要并行执行的话)。
实现代码:
1 unsigned int latestRandNum = 1; 2 3 /*设置随机数种子*/ 4 void srandGLibC(unsigned int seed) 5 { 6 latestRandNum = (seed == 0 ? 1 : seed); 7 } 8 9 /*生成[0, 2147483647]之间的随机数*/ 10 int randGLibC(void) 11 { 12 int randNum = ((latestRandNum * 1103515245) + 12345) & 0x7fffffff; 13 latestRandNum = randNum; 14 15 return (latestRandNum); 16 } 17 18 /*在指定范围内生成随机数*/ 19 int randRange(int min, int max) 20 { 21 unsigned int range = max - min + 1; 22 23 if (range == 1 || range > 2147483648) 24 { 25 return min; 26 } 27 else if (range < 1) 28 { 29 throw std::range_error("The upper limit is less than the lower limit!"); 30 return -1; 31 } 32 33 return (randGLibC() % range + min); 34 }
测试代码:
1 std::bitset<2147483648> randomSet; 2 3 int main(void) 4 { 5 unsigned int maxCount = 2147483648L; 6 7 for (int count = 0; count != 10; ++count) 8 { 9 unsigned int seed = (unsigned int)time(0) + randGLibC(); 10 srandGLibC(seed); 11 randomSet.reset(); 12 std::cout << randomSet.count() << std::endl; 13 std::cout << "seed = " << seed << std::endl; 14 15 for (long long i = 0; i != maxCount; ++i) 16 { 17 int randNum = randRange(0, 2147483647); 18 randomSet[randNum] = true; 19 } 20 21 std::cout << randomSet.count() << std::endl << std::endl; 22 } 23 24 system("pause"); 25 return 0; 26 }
测试结果:
1 0 2 seed = 2456847069 3 2147483648 4 5 0 6 seed = 3328132247 7 2147483648 8 9 0 10 seed = 2984931043 11 2147483648 12 13 0 14 seed = 1708563292 15 2147483648 16 17 0 18 seed = 2508745454 19 2147483648 20 21 0 22 seed = 2926070301 23 2147483648 24 25 0 26 seed = 2299831591 27 2147483648 28 29 0 30 seed = 2616227439 31 2147483648 32 33 0 34 seed = 3271328796 35 2147483648 36 37 0 38 seed = 2417430222 39 2147483648
一共进行了10次测试,每次用不同的随机数种子连续在[0, 2147483647]区间内rand 2147483648次。
而结果是,每次测试中这2147483648个数都正好各被随机到一次。
而把实现代码中
int randNum = ((latestRandNum * 1103515245) + 12345) & 0x7fffffff;
改成
int randNum = ((latestRandNum * 1103515245) + 12344) & 0x7fffffff;
同样用不同的随机数种子进行10次测试,测试结果是:
1 0 2 seed = 2456847067 3 536870912 4 5 0 6 seed = 3268585404 7 268435456 8 9 0 10 seed = 1792994210 11 134217728 12 13 0 14 seed = 2662066861 15 536870912 16 17 0 18 seed = 2955665000 19 268435456 20 21 0 22 seed = 2643173644 23 268435456 24 25 0 26 seed = 2339863270 27 33554432 28 29 0 30 seed = 2380245690 31 134217728 32 33 0 34 seed = 2764481028 35 268435456 36 37 0 38 seed = 1370651278 39 67108864
由上面两次测试结果可知,对于线性同余法的递推公式 X(n+1) = (a * X(n) + c) % m
当m=2147483648时,a=1103515245, c=12345是一组能在[0, 2147483647]保证均匀分布的参数,至于为什么,我不知道。
------------------------------- End -----------------------------