偽隨機數生成算法及性能檢驗

本文永久鏈接為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次，做出圖形，可以看出該方法的特點，以及均勻性好的原因了：

從數據統計和所繪圖形來看，該算法有極好的均勻性，但規律性較強，隨機性較差，適合應用在對均勻度要求較高，而對隨機性要求不高的問題中。

以上給出的幾種偽隨機數生成算法各有優勢，而且性能都能達到一般應用的要求，具有一定的實用價值。