一、解密隨機數產生器（2）——從java源碼看線性同餘演算法

標籤：oca   聲明   getc   匯入   mod   物理   from   數字   代碼   

上篇部落格中，我們瞭解了基於物理現象的真隨機數產生器，然而，真隨機數產生速度較慢，為了實際計算需要，電腦中的隨機數都是由程式演算法，也就是某些公式函數產生的，只不過對於同一隨機種子與函數，得到的隨機數列是一定的，因此得到的隨機數可預測且有周期，不能算是真正的隨機數，因此稱為偽隨機數（Pseudo Random Number）。

 

不過，別看到偽字就瞧不起，這裡面也是有學問的，看似幾個簡簡單單的公式可能是前輩們努力了幾代的成果，相關的研究可以寫好幾本書了！順便提一下，亞裔唯一圖靈獎得主姚期智，研究的就是偽隨機數產生論（The pseudo random number generating theory）。在這裡，我重點介紹兩個常用的演算法：同餘法（Congruential method）和梅森旋轉演算法（Mersenne twister）

 

1、同餘法

 

同餘法（Congruential method）是很常用的一種隨機數產生方法，在很多程式設計語言中有應用，最明顯的就是java了，java.util.Random類中用的就是同餘法中的一種——線性同餘法（Linear congruential method），除此之外還有乘同餘法（Multiplicative congruential method）和混合約餘法（Mixed congruential method）。好了，現在我們就開啟java的原始碼，看一看線性同餘法的真面目！

 

在Eclipse中輸入java.util.Random，按F3轉到Random類的原始碼：

 

首先，我們看到這樣一段說明：

 

翻譯過來是：

 

  這個類的一個實現是用來產生一串偽隨機數。這個類用了一個48位的種子，被線性同餘公式修改用來產生隨機數。（見Donald Kunth《電腦編程的藝術》第二卷，章節3.2.1）

 

顯然，java的Random類使用的是線性同餘法來得到隨機數的。

 

接著往下看，我們找到了它的建構函式與幾個方法，裡麵包含了獲得48位種子的過程：

 

Java代碼  

1. /** 
2.      * Creates a new random number generator. This constructor sets 
3.      * the seed of the random number generator to a value very likely 
4.      * to be distinct from any other invocation of this constructor. 
5.      */  
6.     public Random() {  
7.         this(seedUniquifier() ^ System.nanoTime());  
8.     }  
9.    
10.     private static long seedUniquifier() {  
11.         // L‘Ecuyer, "Tables of Linear Congruential Generators of  
12.         // Different Sizes and Good Lattice Structure", 1999  
13.         for (;;) {  
14.             long current = seedUniquifier.get();  
15.             long next = current * 181783497276652981L;  
16.             if (seedUniquifier.compareAndSet(current, next))  
17.                 return next;  
18.         }  
19.     }  
20.    
21.     private static final AtomicLong seedUniquifier  
22.         = new AtomicLong(8682522807148012L);  
23.     public Random(long seed) {  
24.         if (getClass() == Random.class)  
25.             this.seed = new AtomicLong(initialScramble(seed));  
26.         else {  
27.             // subclass might have overriden setSeed  
28.             this.seed = new AtomicLong();  
29.             setSeed(seed);  
30.         }  
31.     }  
32.     private static long initialScramble(long seed) {  
33.         return (seed ^ multiplier) & mask;  
34.     }  
35.     。。。  

 

這裡使用了System.nanoTime()方法來得到一個納秒級的時間量，參與48位種子的構成，然後還進行了一個很變態的運算——不斷乘以181783497276652981L，直到某一次相乘前後結果相同——來進一步增大隨機性，這裡的nanotime可以算是一個真隨機數，不過有必要提的是，nanoTime和我們常用的currenttime方法不同，返回的不是從1970年1月1日到現在的時間，而是一個隨機的數——只用來前後比較計算一個時間段，比如一行代碼的已耗用時間，資料庫匯入的時間等，而不能用來計算今天是哪一天。

 

好了，現在我不得不佩服這位工程師的變態了：到目前為止，這個程式已經至少進行了三次隨機：

 

1、獲得一個長整形數作為“初始種子”（系統預設的是8682522807148012L）

 

2、不斷與一個變態的數——181783497276652981L相乘（天知道這些數是不是工程師隨便滾鍵盤滾出來的-.-）直到某一次相乘前後數值相等

 

3、與系統隨機出來的nanotime值作異或運算，得到最終的種子

再往下看，就是我們常用的得到隨機數的方法了，我首先找到了最常用的nextInt（）函數，代碼如下：

 

Java代碼  

1. public int nextInt() {  
2.     return next(32);  
3. }  

 

代碼很簡潔，直接跳到了next函數：

 

Java代碼  

1. protected int next(int bits) {  
2.     long oldseed, nextseed;  
3.     AtomicLong seed = this.seed;  
4.     do {  
5.         oldseed = seed.get();  
6.         nextseed = (oldseed * multiplier + addend) & mask;  
7.     } while (!seed.compareAndSet(oldseed, nextseed));  
8.     return (int)(nextseed >>> (48 - bits));  
9. }  

 

 

OK,祝賀一下怎麼樣，因為我們已經深入到的線性同餘法的核心了——沒錯，就是這幾行代碼！

 

在分析這段代碼前，先來簡要介紹一下線性同餘法。

 

在程式中為了使運算式的結果小於某個值，我們常常採用取餘的操作，結果是同一個除數的餘數，這種方法叫同餘法（Congruential method）。

    線性同餘法是一個很古老的隨機數產生演算法，它的數學形式如下：

Xn+1 = (a*Xn+c)(mod m) 

其中，

m>0,0<a<m,0<c<m

 

這裡Xn這個序列產生一系列的隨機數，X0是種子。隨機數產生的品質與m，a，c三個參數的選取有很大關係。這些隨機數並不是真正的隨機，而是滿足在某一周期內隨機分布，這個周期的最長為m。根據Hull-Dobell Theorem，若且唯若：

1. c和m互素;

 

2. a-1可被所有m的質因數整除;

 

3. 當m是4的整數倍，a-1也是4的整數倍

 

時，周期為m。所以m一般都設定的很大，以延長周期。

現在我們回過頭來看剛才的程式，注意這行代碼：

Java代碼  

1. nextseed = (oldseed * multiplier + addend) & mask;  

和Xn+1=(a*Xn+c)(mod m)的形式很像有木有！

沒錯，就是這一行代碼應用到了線性同餘法公式！不過還有一個問題：怎麼沒見取餘符號？嘿嘿，先讓我們看看三個變數的數值聲明：

 

Java代碼  

1. private static final long multiplier = 0x5DEECE66DL;  
2. private static final long addend = 0xBL;  
3. private static final long mask = (1L << 48) - 1;  

 

其中multiplier和addend分別代表公式中的a和c，很好理解，但mask代表什麼呢？其實，x & [(1L << 48)–1]與 x（mod 2^48）等價。解釋如下：

 

x對於2的N次冪取餘，由於除數是2的N次冪，如：

0001，0010，0100，1000。。。。

相當於把x的二進位形式向右移N位，此時移到小數點右側的就是餘數，如：

13 = 1101    8 = 1000

13 / 8 = 1.101，所以小數點右側的101就是餘數，化成十進位就是5

然而，無論是C語言還是java,位元運算移走的數顯然都一去不複返了。（什麼，你說在CF寄存器中？好吧，太高端了點，其實還有更給力的方法）有什麼好辦法保護這些即將逝去的資料呢？

學著上面的mask，我們不妨試著把2的N次冪減一：

0000，0001，0011，0111，01111，011111。。。

怎麼樣，有啟發了嗎？

我們知道，某個數（限0和1）與1作與（&）操作，結果還是它本身；而與0作與操作結果總是0，即：

a & 1 = a,  a & 0 = 0

而我們將x對2^N取餘操作希望達到的目的可以理解為：

1、所有比2^N位（包括2^N那一位）高的位全都為0

2、所有比2^N低的位保持原樣

因此， x & （2^N-1）與x（mod 2^N）運算等價，還是13與8的例子：

1101 % 1000 = 0101    1101 & 0111 = 0101

二者結果一致。

嘿嘿，講明白了這個與運算的含義，我想上面那行代碼的含義應該很明了了，就是線性同餘公式的直接套用，其中a = 0x5DEECE66DL, c = 0xBL, m = 2^48，就可以得到一個48位的隨機數，而且這個謹慎的工程師進行了迭代，增加結果的隨機性。再把結果移位，就可以得到指定位元的隨機數。

接下來我們研究一下更常用的一個函數——帶參數n的nextInt：

 

Java代碼  

1. public int nextInt(int n) {  
2.     if (n <= 0)  
3.         throw new IllegalArgumentException("n must be positive");  
4.   
5.     if ((n & -n) == n)  // i.e., n is a power of 2  
6.         return (int)((n * (long)next(31)) >> 31);  
7.   
8.     int bits, val;  
9.     do {  
10.         bits = next(31);  
11.         val = bits % n;  
12.     } while (bits - val + (n-1) < 0);  
13.     return val;  
14. }  

 

顯然，這裡基本的思路還是一樣的，先調用next函數產生一個31位的隨機數（int類型的範圍），再對參數n進行判斷，如果n恰好為2的方冪，那麼直接移位就可以得到想要的結果；如果不是2的方冪，那麼就關於n取餘，最終使結果在[0,n)範圍內。另外，do-while語句的目的應該是防止結果為負數。

你也許會好奇為什麼(n & -n) == n可以判斷一個數是不是2的次方冪，其實我也是研究了一番才弄明白的，其實，這主要與補碼的特性有關：

眾所周知，電腦中負數使用補碼儲存的（不懂什麼是補碼的自己百度惡補），舉幾組例子：

2 ：0000 0010      -2 ：1111 1110

8 ：0000 1000      -8 ：1111 1000

18 ：0001 0010     -18 ：1110 1110

20 ：0001 0100     -20 ：1110 1100

不知道大家有沒有注意到，補碼有一個特性，就是可以對於兩個相反數n與-n，有且只有最低一個為1的位元字相同且都為1，而更低的位全為0，更高的位各不相同。因此兩數作按位與操作後只有一位為1，而能滿足這個結果仍為n的只能是原本就只有一位是1的數，也就是恰好是2的次方冪的數了。

不過個人覺得還有一種更好的判斷2的次方冪的方法：

n & (n-1) == 0

感興趣的也可以自己研究一下^o^。

好了，線性同餘法就介紹到這了，下面簡要介紹一下另一種同餘法——乘同餘法（Multiplicative congruential method）。

上文中的線性同餘法，主要用來產生整數，而某些情景下，比如科研中，常常只需要（0，1）之間的小數，這時，乘同餘法是更好的選擇，它的基本公式和線性同餘法很像：

Xn+1=（a*Xn ）(mod m ）

其實只是令線性公式中的c=0而已。只不過，為了得到小數，我們多做一步：

Yn = Xn/m  

由於Xn是m的餘數，所以Yn的值介於0與1之間，由此到（0，1）區間上的隨機數列。

除此之外，還有混合約餘法，二次同餘法，三次同餘法等類似的方法，公式類似，也各有優劣，在此不詳細介紹了。

同餘法優勢在計算速度快，記憶體消耗少。但是，因為相鄰的隨機數並不獨立，序列關聯性較大。所以，對於隨機數品質要求高的應用，特別是很多科研領域，並不適合用這種方法。

不要走開，下篇部落格介紹一個更給力的演算法——梅森旋轉演算法（Mersenne Twister），持續關注啊！

一、解密隨機數產生器（2）——從java源碼看線性同餘演算法（轉）