伪随机数生成算法及性能检验

本文永久链接为http://johnhany.net/2013/11/random-algorithm-and-performance/ 转载请注明出处

什么叫伪随机数

在一些问题中，比如计算机仿真和模拟、密码学等应用中，需要产生一个随机数数列来解决问题。

随机数数列分为真随机数数列和伪随机数数列。真随机数数列是完全不可预测的，可以通过放射性衰变、电子设备的热噪音、宇宙射线的触发时间等物理过程得到，但无法通过单纯的软件方式获得；伪随机数数列是可预测的，严格意义上不具有随机性质，通常用数学公式的方法获得。

由于很多应用当中对“随机”数列的要求并不十分严格，而且伪随机数的产生较真随机数更为廉价、高效，所以在大多数情况下只要产生统计分布较好的伪随机数数列就能够满足应用要求。早期的伪随机数生成算法以平方取中法为主，现在则以线性同余法为主要方式。下面我会就两种方式分别给出实例，并以数据统计和图形化的方式对伪随机数生成算法的性能进行检验。

种子

正如数列需要有首项，产生伪随机数需要一个初值用来计算整个序列，这个初值被称为“种子”。种子可以是一个固定的值，也可以是根据当前系统状态确定的值。C语言用来产生伪随机数的库函数rand()的种子是固定的值，因此每次调用该函数产生的随机数数列都是相同的。所以为了获得随机性更好的数列，种子应为一个变量，该变量可以与当前系统时间、用户对键盘或鼠标操作情况有关。这里我将根据系统时间获得种子。
代码如下：

　　　#include <time.h>
　　　unsigned long ssecond,nowtime;
unsigned long seed;
　　　long gettime()          //获得当前时间
　　　{
　　　    time_t t;
　　　    time(&t);
　　　    struct tm *local;
　　　    local=localtime(&t);
　　　    local->tm_mon++;
　　　    ssecond=(long)local->tm_sec*100+local->tm_sec+36923;
　　　    nowtime=local->tm_sec + local->tm_min*100 + local->tm_hour*10000 + local->tm_mday*1000000 + local->tm_mon*100000000;
　　　    return nowtime;
　　　}

在调用伪随机数生成函数之前通过seed=gettime();语句就完成了种子的初始化。

平方取中法

平方取中法是由冯·诺依曼在1946年提出的，其基本思想为：将数列中的第a(i)项（假设其有m位）平方，取得到的2m位数（若不足2m位，在最高位前补0）中间部分的m位数字，作为a(i)的下一项a(i+1)，由此产生一个伪随机数数列。即：

x(i+1)=(10^(-m/2)*x(i)*x(i))mod(10^m)

平方取中法计算较快，但在实际应用时会发现该方法容易产生周期性明显的数列，而且在某些情况下计算到一定步骤后始终产生相同的数甚至是零，或者产生的数字位数越来越小直至始终产生零。所以用平方取中法产生伪随机数数列时不能单纯使用公式，应该在计算过程中认为加入更多变化因素，比如根据前一个数的奇偶性进行不同的运算，如果产生的数字位数减少时通过另一种运算使其恢复成m位。

代码如下：

　//编译环境：Code::Blocks 10.05
　　　//生成的数字范围为[0,LENGTH-1]
　　　long intlen(long in)        //整数in的长度
　　　{
　　　    long count=0;
　　　    while(in!=0)
　　　    {
　　　        in/=10;
　　　        count++;
　　　    }
　　　    return count;
　　　}
　　　long power_10(long in)      //10的in次幂
　　　{
　　　    long i,out=1;
　　　    for(i=0;i<in;i++)
　　　        out*=10;
　　　    return out;
　　　}
　　　long rand_pfqz(void)      //平方取中
　　　{
　　　    long len;
　　　    while(seed<10000)			//保持数位一致
　　　        seed=seed*13+seed/10+second/3;
　　　    len=intlen(seed);
　　　    long temp=seed;
　　　    temp=((seed*seed/power_10(len/2))%power_10(len));
　　　    if(temp%2==0)   temp+=second/3+7854;		//增加改变因素，
　　　    else    temp+=second*second/2317;			//以延长计算周期
　　　    seed=temp;
　　　    return (unsigned long)(seed%10000*LENGTH)/10000;
　　　}

以下为在5月18日11时12分43秒时刻以LENGTH=100产生的随机数数列：

93 12 7 4 70 32 95 97 25 8
49 40 85 63 11 53 65 4 42 76
33 1 34 88 31 37 79 1 62 97
59 21 6 47 28 1 98 60 56 53
61 48 47 31 8 54 53 15 54 67
18 6 46 45 87 44 55 54 46 24
74 12 68 41 97 18 27 86 13 81
99 74 49 52 69 11 52 89 33 7
73 22 1 95 19 89 57 21 77 90
70 87 47 59 19 26 89 32 44 33

可见，该算法在连续计算100次时取到了0~99之间的65个不同的数。

LENGTH=500，以连续产生的两个数字作为平面上点的横坐标与纵坐标，计算2000次，所得图形如下：

线性同余法

线性同余方法是目前应用广泛的伪随机数生成算法，其基本思想是通过对前一个数进行线性运算并取模从而得到下一个数。即：
a(i+1)=(a(i)*b+c)mod(m)

其中b称为乘数，c称为增量，m称为模数，它们均为常数。

乘数、增量和模数的选取可以多种多样，只要保证产生的随机数有较好的均匀性和随机性即可。

线性同余法的最大周期是m，但一般情况下会小于m。要使周期达到最大，应该满足以下条件：

(1) c和m互质；

(2) m的所有质因子的积能整除b-1；

(3) 若m是4的倍数，则b-1也是；

(4) b，c，a(0)（初值，一般即种子）都比m小；

(5) b，c是正整数。

在C和VC中都定义有用于产生伪随机数的库函数rand()，而且都是利用线性同余法产生伪随机数。
在C中的rand()函数定义如下：

　　　#define RANDOM_MAX 0x7FFFFFFF
static long do_rand(unsigned long *value)
{
long quotient, remainder, t;

quotient = *value / 127773L;
remainder = *value % 127773L;
t = 16807L * remainder - 2836L * quotient;
if (t <= 0)
t += 0x7FFFFFFFL;
return ((*value = t) % ((unsigned long)RANDOM_MAX + 1));
}
　　static unsigned long next = 1;
　　int rand(void)
　　{
　　   return do_rand(&next);
　　}
　　void srand(unsigned int seed)			//赋初值为种子
　　{
　　   next = seed;
　　}

在VC中rand()函数定义如下：

　　int __cdecl rand (void) {       return(((holdrand = holdrand * 214013L + 2531011L) >> 16) & 0x7fff);   }

C程序设计语言（第二版）(Brian W. Kernighan, Dennis M. Ritchie.)中给出了一种适于C的rand()函数：

　　unsigned long next=1;
　　int rand(void)
　　{
　　    next=next*1103515245+12345;
　　    return (unsigned int)(next/65536)%32768;
　　}
　　void srand(unsigned int seed)
　　{
　　　　　next=seed;
　　}

为了提高库函数rand()的性能，可以通过以下函数进行再次运算产生数列：

　　int rrand(int n)
　　{
　　    return 1+(int)(n*rand()/(RAND_MAX+1));
　　}

以下为LENGTH=100，利用C的rand()库函数产生的随机数序列：

41 67 34 0 69 24 78 58 62 64
5 45 81 27 61 91 95 42 27 36
91 4 2 53 92 82 21 16 18 95
47 26 71 38 69 12 67 99 35 94
3 11 22 33 73 64 41 11 53 68
47 44 62 57 37 59 23 41 29 78
16 35 90 42 88 6 40 42 64 48
46 5 90 29 70 50 6 1 93 48
29 23 84 54 56 40 66 76 31 8
44 39 26 23 37 38 18 82 29 41
可见，该算法在连续计算100次时取到了0~99之间的64个不同的数。
为了提高算法的均匀性，我找到了一种均匀度非常好的线性同余算法。C代码如下：

　　　//编译环境：Code::Blocks 10.05
　　　//与平方取中法类似，种子根据当前系统时间获取
　　　unsigned long nowtime;
unsigned long seed;
　　　long gettime()          //获得当前时间
　　　{
　　　    time_t t;
　　　    time(&t);
　　　    struct tm *local;
　　　    local=localtime(&t);
　　　    local->tm_mon++;
　　　    nowtime=local->tm_sec + local->tm_min*100 + local->tm_hour*10000 + local->tm_mday*1000000 + local->tm_mon*100000000;
　　　    return nowtime;
　　　}
　　　#define LENGTH 100
　　　//生成的数字范围为[0,LENGTH-1]
　　int rand_xxty(void)     //线性同余法
　　{
　　    seed=(unsigned long)(seed*101+81)%LENGTH;
　　    return (int)seed;
　　}

以下为在5月18日12时51分27秒时刻以LENGTH=100产生的伪随机数：

56 37 18 99 80 61 42 23 4 85
66 47 28 9 90 71 52 33 14 95
76 57 38 19 0 81 62 43 24 5
86 67 48 29 10 91 72 53 34 15
96 77 58 39 20 1 82 63 44 25
6 87 68 49 30 11 92 73 54 35
16 97 78 59 40 21 2 83 64 45
26 7 88 69 50 31 12 93 74 55
36 17 98 79 60 41 22 3 84 65
46 27 8 89 70 51 32 13 94 75
可见，该算法在连续计算100次时取到了0~99之间的所有的数。

LENGTH=500，以连续产生的两个数字作为平面上点的横坐标与纵坐标，计算2000次，做出图形，可以看出该方法的特点，以及均匀性好的原因了：



从数据统计和所绘图形来看，该算法有极好的均匀性，但规律性较强，随机性较差，适合应用在对均匀度要求较高，而对随机性要求不高的问题中。
以上给出的几种伪随机数生成算法各有优势，而且性能都能达到一般应用的要求，具有一定的实用价值。