Redis源碼中看偽隨機數產生演算法

標籤：drand48   lrand48   偽隨機數產生演算法   隨機數種子   redis源碼   

導讀--------------------------------------------------------------------------------------------------------------------------------------------------------------

        Redis源碼中有一個rand.c的源檔案，很明顯這是一個和（偽）隨機數有關的檔案。細看該檔案代碼只有寥寥50行，不過涉及到的演算法原理卻不簡單，讀起來雖然有些晦澀，但對於深入理解48位空間中的偽隨機數演算法是不可多得的範本。作者在該檔案的注釋中寫道：這個偽隨機數產生函數是從pysam源碼中的drand48()派生過來的。關於pysam是什麼項目，並不是重點，其實很多Unix系統中都存在drand48這個函數（SVr4，POSIX.1-2001），我們可在終端中man一下drand48。

        可以看到在標頭檔stdlib.h中定義了一系列和隨機數有關的函數，除drand48()外，還有erand48()、lrand48()、nrand48()等等。所謂的48指的是隨機數種子的二進位位元。它們的主要操作使用的都是同一個隨機數產生演算法，但是後續的處理不同，比如drand48()、erand48()返回的是 [0.0, 1.0)之間的浮點型，而lrand48()和nrand48()返回的是[0, 2^31)之間的整數。rand.c檔案中有一個重要的函數——redisLrand48()。雖然作者是受了drand48()源碼的啟發編寫的該函數，但實際上redisLrand48()和lrand48()更像，擁有幾乎一樣的編程介面，並且設定相同的隨機數種子，產生的偽隨機數序列相同。

        繼續讀rand.c開頭的注釋，可以發現作者編寫這個檔案的目的是用於替換Lua中實現的math.random()函數。作者認為Lua隨機數產生演算法的預設實現中使用了libc的rand()函數，而libc的rand()並不能保證在設定了相同種子的情況下，不同平台能產生相同的隨機數序列。因此作者重寫這一檔案，保證跨平台性。

--------------------------------------------------------------------------------------------------------------------------------------------------------------

演算法原理線性同餘方程        該偽隨機數產生演算法使用了線性同餘方程，顧名思義：所謂 “線性”，指的是自變數（方程中的x）和因變數（方程中的y）是線性關係；所謂 “同餘”，表示它們與同一個數（比如m）相除，餘數相同。可以表示為：y = f(x)%m 或 y = f(x) mod m 【f(x)是線性函數，即x的次數為1】        在Linux中，通過lrand48的手冊頁，可以看到lrand48()所使用的線性同餘方程的具體內容（Redis使用的是同一個方程）：

Xn+1 = (aXn + c) mod m, where n >= 0//預設情況下a = 0x5DEECE66Dc = 0xBm = 2^48 

X0、X1、X2……Xn就是偽隨機數序列。它們每一個都要由上一個數字通過運算來產生。
種子        隨機數產生演算法中都有種子的概念，無論是普通的rand()還是lrand48()。我們前面多次提到 偽隨機數這一術語，之所以說“偽”，是因為當種子確定之後，此後產生到隨機數序列也就確定了。換句話說只要是相同的種子，產生的隨機數序列總是固定的（相同的平台下），並非完全隨機。        通過上一小節中的線性同餘方程，我們觀察到隨機數序列的產生受到X0、a、c的影響（m總是2^48，表示48位空間中的運算，所以這些函數尾碼才會有48）。因而X0、a、c其本質就是種子。POSIX的lrand48()的實現中和Redis的實現中給這三個參數提供了相同預設值。但不同之處是POSIX標準提供了lcong48()函數來修改這三個參數的預設值，而Redis中僅僅能修改X0（表示隨機數序列的第一個數）的值，所以在Redis中真正的只有一個——X0.源碼實現變數        rand.c中定義了三個靜態全域變數：x、a、c 與線性同餘方程中的參數相對應

static uint32_t x[3] = { X0, X1, X2 }, a[3] = { A0, A1, A2 }, c = C;

這裡涉及到幾個宏定義：

#define X00x330E#define X10xABCD#define X20x1234#define A00xE66D#define A10xDEEC#define A20x5#define C0xB

        可以看出數組a中的三個元素儲存的就是線性同餘方程中的常數 0x5DEECE66D的三個部分。注意的是a中每個元素是uint32_t類型，就是說4個位元組，實際上我們把0x5DEECE66D拆成三個部分（0x5、0xDEEC、0xE66D），每個部分最多隻佔用了2個位元組。另外要注意的是宏定義中的X0和前文（上一節）中的X0語義是不同的，這裡的X0表示的是數組X預設初始化時的第一個元素的內容，而上一節中提到的X0表示的是第一個隨機數，相當於本節中的X2*2^32 + X1*2^16 + X0。
對外函數        rand.c對外只提供了兩個函數：

• redisLrand48()：返回一個隨機數
• redisSrand48()：設定隨機數種子

另外還有一個靜態函數next()，它是為redisLrand48()服務的，對外不可見。redisSrand48()函數

void redisSrand48(int32_t seedval) {    SEED(X0, LOW(seedval), HIGH(seedval));}

        設定隨機數的種子，用到了宏函數SEED。看一下定義：

#define SET3(x, x0, x1, x2)((x)[0] = (x0), (x)[1] = (x1), (x)[2] = (x2))#define SEED(x0, x1, x2) (SET3(x, x0, x1, x2), SET3(a, A0, A1, A2), c = C)

        可以看成宏函數SEED進行了三個賦值操作，分別是給數組x、數組a和變數c。不過a和c賦的都是預設值，所以SEED(X0,LOW(seedval),HIGH(seedval))實際進行的操作就是用X0（0x330E）、seedval的低16位，seedval的高16位分別賦值給x[0]、x[1]、x[2]。LOW和HIGH也是兩個宏函數

#define N16 #define MASK((1 << (N - 1)) + (1 << (N - 1)) - 1) //MASK=65535(16個二進位位1)#define LOW(x)((unsigned)(x) & MASK)                //LOW(x)獲得3x的低16位（2個位元組）#define HIGH(x)LOW((x) >> N)                         //HIGH(x)獲得x的高16位（2個位元組）

注意上面所指的高16位和低16位指的是32位整型資料的高低。最終每次要保留的結果是48位的（儲存在數組x[ ]中），它分為高16位(x[2])、中16位(x[1])和低16位(x[0])。不要混淆。

redisLrand48()函數

int32_t redisLrand48() {    next();    return (((int32_t)x[2] << (N - 1)) + (x[1] >> 1));}

        next函數容後再稟，先看一下返回值，把N=16帶進去，就是：

((int32_t)x[2] << 15) + (x[1] >> 1)

就是x[2]*2^15 + x[1]/2。這裡沒有什麼道理可講，隨機數嘛，這隻是一種方案而已（lrand48()採用的方案）。next()函數公式推導與化簡        next()是為redisLrand48()服務的static函數，也就是說對外不可見的，但是next()函數才是偽隨機數產生演算法的精髓所在。再回顧一遍線性同餘公式：

Xn+1 = (aXn + c) mod m, where n >= 0,m = 2^48

        注意：n和n+1都是X的下標。從數學角度來看很簡單啊，先乘再加，最後模數。但是電腦做起來卻不盡然，因為會溢出，所以我們是用數組儲存的x和a。所以要做的是用數組類比兩個數的乘法運算。先把方程中的X和a表示出來：


對線性同餘方程進行一下換算：

Xn+1 = (aXn + c) mod mXn+1 = ((aXn) mod m + c mod m) mod mXn+1 = ((aXn) mod m + c) mod m

        先計算一下a*Xn：
        這個多項式共有9個部分組成。再計算一下(a*Xn) mod m。因為m = 2^48，所以上述多項式中，2的冪次大於等於48的項一定是2^48的倍數，模數之後就是0，故可去掉這些項。所以(a*Xn)mod m的結果為：
        只剩下6項了。接著合并同類項，並加上c，很容易寫出 a*Xn + c的結果：
        上述算式的每個括弧中的內容就是在一次隨機數產生之後x[2]、x[1]、x[0] 裡面儲存的新的內容（即用來表示Xn+1）。當然了編程去計算的時候還要考慮進位的問題。先看一個宏定義：用到的宏

#define MUL(x, y, z){ int32_t l = (long)(x) * (long)(y); (z)[0] = LOW(l); (z)[1] = HIGH(l); }#define CARRY(x, y)((int32_t)(x) + (long)(y) > MASK)#define ADDEQU(x, y, z)(z = CARRY(x, (y)), x = LOW(x + (y)))

        很明顯宏函數MUL(x,y,z)表示的是乘法運算，x和y相乘，結果用z來儲存。z[0]儲存低16位，z[1]儲存高16位。        CARRY(x,y)就是判斷兩個數相加是否會進位”，這裡的“進位”指的超出2個位元組（16位）的大小。因為我們無論是數組a還是數組x，其中每個元素儲存的都是2個位元組的有效資料。也就是說，如果x[0]裡面的數字超過2個位元組，就要把多出的部分“進位”到x[1]中。每個元素其實都是4個位元組的大小，之所以沒有用到全部的4位元組來儲存，是因為當4位元組全部用來儲存資料時，相加以後可能就會溢出，編譯器會彈出警告，其結果也會變成負數。
        ADDEQU(x,y,z)執行的操作是：x和y相加，其結果儲存到x中。z中儲存是否“進位”。next()源碼

static void next(void) {    uint32_t p[2], q[2], r[2], carry0, carry1;    MUL(a[0], x[0], p);    ADDEQU(p[0], c, carry0);    ADDEQU(p[1], carry0, carry1);    MUL(a[0], x[1], q);    ADDEQU(p[1], q[0], carry0);    MUL(a[1], x[0], r);    x[2] = LOW(carry0 + carry1 + CARRY(p[1], r[0]) + q[1] + r[1] +            a[0] * x[2] + a[1] * x[1] + a[2] * x[0]);    x[1] = LOW(p[1] + r[0]);    x[0] = LOW(p[0]);}

        我們來一行一行的分析一下。聲明語句不看，從第4行開始讀。

    MUL(a[0], x[0], p);    ADDEQU(p[0], c, carry0);    ADDEQU(p[1], carry0, carry1);

• MUL運算：a[0]和x[0]相乘，結果儲存到數組p中。
• 第一個ADDEQU：然後p[0]再和c相加，carry0標識是否有“進位”。這裡完成的就是運算。
• 第二個ADDEQU之後就是把低16位加上c之後是否進位，加到p[1]，並且判斷這次相加之後有沒有產生新的進位（指2^16次多項式向2^32次多項式的進位）並儲存到carry1。

至此一次多項式相關的運算計數完畢，接下來是2^16次多項式的運算。

    MUL(a[0], x[1], q);    ADDEQU(p[1], q[0], carry0);    MUL(a[1], x[0], r);

• 第一個MUL運算：a[0]和x[1]相乘，結果儲存到數組q中。
• ADDEQU運算：將p[1]（儲存了“低16位”運算之後向“中16位”進位的值）和q[0]相加，結果儲存到p[1]，carry0標識“中16位”（2^16次多項式）對高16位（2^32次多項式）是否有新的進位。
• 第二個MUL運算：a[1]和x[0]相乘，結果儲存到數組r中。

至此完成了係數相關運算。

    x[2] = LOW(carry0 + carry1 + CARRY(p[1], r[0]) + q[1] + r[1] +            a[0] * x[2] + a[1] * x[1] + a[2] * x[0]);    x[1] = LOW(p[1] + r[0]);    x[0] = LOW(p[0]);

        倒著讀吧，先看x[0]，它儲存的應該是一次多項式的結果。        x[1]儲存2^16次多項式的係數。p[1]裡面儲存著a[0]x[1]的結果以及低16位的進位。r[0]儲存著a[1]x[0]的結果（不包括進位）。
        x[2]最複雜。回顧一下2^32次多項式的係數：再看一眼，那行代碼，這幾個兩兩直接相乘的部分讀懂了吧，那麼可以去掉他們再看其他部分：carry0 + carry1 +CARRY(p[1], r[0]) + q[1] + r[1]剩餘的多項式就是2^16次多項式在運算過程中向上的進位了。
        畫了一個圖給大家表示一下： 
        這個圖的上面兩行是表頭，下面彩色部分是所儲存的內容。比如第三列表示x[0]裡面儲存的是p[0]和c的和。        從第一幅圖到第二幅圖的轉變是因為執行了 ADDEQU(p[1], q[0], carry0);這條語句，將q[0]加到了p[1]上。然後再去看看x[2]的那段代碼，就清楚多了吧，重點是考慮進位。

Redis源碼中看偽隨機數產生演算法