如何产生 [0, 2147483647] 之间的随机数

一、简介

随机数是编程中经常要用到的东西,但很遗憾的是在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]之间的随机数。

实现代码如下:

View Code
 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  }

看上去这个思路应该是没有问题的,因为每一位都是随机的,那整个数一定是随机的,但是测试结果不尽如人意。

测试代码:

View Code
 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

View Code
 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

View Code
1 /*   /stdlib/rand.c   */
2 int
3 rand ()
4 {
5   return (int) __random ();
6 }

/stdlib/random.c

View Code
 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

View Code
 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()函数的简单实现

下面是一个简单的版本,忽略了临界区访问控制等代码(如果不需要并行执行的话)。

实现代码:

View Code
 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 }

测试代码:

View Code
 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 }

测试结果:

View Code
 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次测试,测试结果是:

View Code
 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 -----------------------------

posted on 2012-11-20 20:12  Snser  阅读(2180)  评论(3编辑  收藏  举报