談談梅森旋轉：演算法及其爆破

01-29

原文鏈接：談談梅森旋轉：演算法及其爆破
知乎的編輯器真是太爛了，公式沒法看……所以……建議移步原文。

現代編程語言，大都在標準庫中包含了隨機庫。例如，C++ 在 C++11 標準中添加了 random 頭文件，提供了現代的隨機庫；Python 則有 random。C++11 的隨機庫將生成隨機數的過程在邏輯上切分成了兩個步驟：隨機數生成引擎和分布。在學習 C++11 的 random 庫時，std::mt19937 這一隨機數生成引擎的名字看起來十分奇怪，成功吸引了我的注意力。

查詢後得知，std::mt19937 中的 MT 是 Mersenne Twister 的縮寫，這是偽隨機數生成演算法的名字（梅森旋轉演算法）；而 19937 則取自演算法中用到的梅森素數 $2^{19937} - 1$ 這裡，梅森素數是演算法生成偽隨機數的循環長度（period），而旋轉則說的是演算法內部對定長二進位串循環位移的過程。

此篇講解梅森旋轉演算法的一些原理，並介紹對其的一個「爆破」方法。

偽隨機數發生器質量的度量—— $k$ -維 $v$ -比特準確度

梅森旋轉演算法（Mersenne Twister Algorithm，簡稱 MT）是為了解決過去偽隨機數發生器（Pseudo-Random Number Generator，簡稱 PRNG）產生的偽隨機數質量不高而提出的新演算法。該演算法由松本眞（Makoto Matsumoto）和西村拓士（Takuji Nishimura）在 1997 年提出，期間還得到了「演算法之神」高德納（Donald Ervin Knuth）的幫助。

既然 MT 是為了解決過去 PRNG 質量低下的問題，那麼首先我們就必須要有一個能夠度量 PRNG 質量的方法。否則，「公說公有理婆說婆有理」，我們就無法對 PRNG 作出統一的評價了。

這裡介紹評價 PRNG 最嚴格指標： $k$ -維 $v$ -比特準確度（ $k$ -distributed to $v$ -bit accuracy）。

假設某一 PRNG 能夠產生周期為 PP 的 ww-比特的隨機數序列 {x? i}{x→i}；同時，將 ww-比特的隨機數 x? x→ 的最高 vv 位組成的數記作 truncv(x? )truncv(x→)。現構造如下長度為 kvkv-比特的二進位數

PRNGk,v(i)=def(truncv(x? i),truncv(x? i+1),…,truncv(x? i+k?1)).PRNGk,v(i)=def(truncv(x→i),truncv(x→i+1),…,truncv(x→i+k?1)).

由於 PRNGk,v(i)PRNGk,v(i) 是長度為 kvkv-比特的二進位數，所以它可以有 2kv2kv 中不同的取值。若當 ii 遍歷 [0,P)[0,P)，PRNGk,v(i)PRNGk,v(i) 的取值在這 2kv2kv 中均勻分布。具體來說，PRNGk,v(i)PRNGk,v(i) 落在這 2kv2kv 個取值上的數量完全相等（P+12kvP+12kv），除了全零的取值會少不多不少正好一次（P+12kv?1P+12kv?1）。

顯而易見，對任意固定的 vv，若 PRNG 是 kk-維 vv-比特準確的，那麼必然也是 (k?1)(k?1)-維 vv-比特準確的，但不一定是 (k+1)(k+1)-維 vv-比特準確的。考慮 kk 是有上限的，因此，對於任意固定的 vv，必然存在最大的 k=k(v)k=k(v) 使得 PRNG 是 k(v)k(v)-維 vv-比特準確的。那麼，根據定義，顯然有 2k(v)v?1?P2k(v)v?1?P。

kk-維 vv-比特準確度也可以有密碼學角度的描述：若 PRNG 是 kk-維 vv-比特準確的，那麼即使已知 PRNG 生成的 l<kl<k 個偽隨機數的最高 vv 位，也無法推出 PRNG 生成的第 l+1l+1 個偽隨機數的最高 vv 位。

根據這樣的定義，MT 演算法具有非常優良的性能。首先 MT 演算法（MT19937）的周期非常長。其周期 P=219937?1P=219937?1，要比預估的宇宙可觀測的粒子總數（10871087）還要高出數千個數量級。其次，作為一個 3232 比特隨機數生成器，MT19937 是 623623-維 3232-比特準確的。考慮到 ?1993732?=623?1993732?=623，MT19937 在 kk-維 vv-比特準確度上的性能已達到理論上的最大值。因此，MT19937 具有非常優良的性能。

梅森旋轉演算法的描述

旋轉

32 位的梅森旋轉演算法能夠產生周期為 PP 的 ww-比特的隨機數序列 {x? i}{x→i}；其中 w=32w=32。這也就是說，每一個 x? x→ 是一個長度為 32 的行向量，並且其中的每一個元素都是二元數域 F2=def{0,1}F2=def{0,1} 中的元素。現在，我們定義如下一些記號，來描述梅森旋轉演算法是如何進行旋轉（線性移位）的。

nn：參與梅森旋轉的隨機數個數；
rr：[0,w)[0,w) 之間的整數；
mm：(0,n](0,n] 之間的整數；
AA：w×ww×w 的常矩陣；
x? (u)x→(u)：x? x→ 的最高 w?rw?r 比特組成的數（低位補零）；
x? (l)x→(l)：x? x→ 的最低 rr 比特組成的數（高位補零）。

梅森旋轉演算法，首先需要根據隨機數種子初始化 nn 個行向量：

x? 0,x? 1,…,x? n?1.x→0,x→1,…,x→n?1.

而後根據下式，從 k=0k=0 開始依次計算 x? nx→n：

x? k+n=defx? k+m⊕(x? (u)k∣x? (l)k+1)A.(1)(1)x→k+n=defx→k+m⊕(x→k(u)∣x→k+1(l))A.

其中，x? ∣x? ′x→∣x→′ 表示兩個二進位數按位或；x? ⊕x? ′x→⊕x→′ 表示兩個二進位數按位半加（不進位，也就是按位異或）；x? Ax→A 則表示按位半加的矩陣乘法。在 MT 中，AA 被定義為

?????????aw?11aw?21aw?3??1a0?????????(11?1aw?1aw?2aw?3?a0)

因此

x? A={x? >>1(x? >>1)⊕a? if x0=0if x0=1.x→A={x→>>1if x0=0(x→>>1)⊕a→if x0=1.

此處，若 r=0r=0，則 MT 退化為 TGFSR (Matsumoto and Kurita 1992, 1994)；若再加上 A=IwA=Iw，則又從 TGFSR 退化為 GFSR (Lewis and Payne 1973)。

因此，梅森旋轉 11 式完全由位運算組成（移位、按位與、按位或、按位異或）。

線性反饋移位寄存器、旋轉之名、周期

上一小節我們介紹了 MT 演算法當中的「旋轉」。但只憑抽象的數學公式（尤其是二進位的邏輯數學），很難看出它為什麼是「旋轉」。這一節我們首先介紹線性反饋移位寄存器（Linear Feedback Shifting Register，簡稱 LFSR），看到它是如何「旋轉」的；最後再將 LFSR 和 MT 演算法當中的旋轉步驟統一起來。

反饋移位寄存器是對二進位序列進行等長變換的一種特殊函數。它包括兩個部分：

級。等長変換的長度即是反饋移位寄存器的級。
反饋函數。若反饋函數是線性的，則稱線性反饋移位寄存器；否則是非線性反饋移位寄存器。

一般來說，LFSR 是基於異或運算的。一個 LFSR 的工作步驟是這樣的：

將原寄存器狀態的最低位作為輸出。
執行線性反饋函數；也就是選取其中若干位，從高位到低位迭代異或。
將元寄存器狀態向低位移位 1 位，並以上述迭代異或的結果作為填充值填入最高位。

在不斷的迭代異或、填充高位的過程中，二進位位在寄存器中循環旋轉。這就是「旋轉」的來由。

對於一個 4 級的 LFSR 來說，假設其反饋函數是 f(x)=x4+x2+x+1f(x)=x4+x2+x+1。則 LFSR 每次從最低位取出結果，將最高位（x4x4）和倒數第二低位（x2x2）取異或後，再與最低位（xx）取異或後，填入移位後的最高位。

因此，若初始狀態為 (1000)2(1000)2，則有

1000n1100n1110n0011n0001n1000n

如此，我們就構建了一個循環長度為 5 的 LFSR。

考慮一個 ww 級的 LFSR，其最多共有 2w2w 種狀態。又考慮對於異或來說，全零的狀態是不可用的（因為不論如何運算都是全零），因此一個 ww 級的 LFSR，其有效狀態共有 2w?12w?1 個。因此，理論上，一個 LFSR 的循環長度最大為 2w?12w?1。可以證明，當且僅當反饋函數是 F2F2 上的本原多項式（素多項式）時，LFSR 的循環長度達到最大值。

回過頭來看 11 式，不難發現，這其實相當於一個 nw?rnw?r 級的線性反饋移位寄存器（取 x? (u)kx→k(u) 的最高 w?rw?r 位與 x? (l)k+1x→k+1(l) 的最低 rr 位進行迭代異或，再經過一個不影響周期的線性變換 AA）。只不過，11 式每一次運算，相當於 LFSR 進行了 ww 輪計算。若 ww 與 nw?rnw?r 互素，那麼這一微小的改變是不會影響 LFSR 的周期的。考慮到 LFSR 的計算過程像是在「旋轉」，這即是「梅森『旋轉』」名字的來由。

對這個等效的 nw?rnw?r 級 LFSR 來說，當且僅當反饋函數是 F2F2 上的本原多項式（素多項式）時，MT 的循環周期長度 PP 達到最大值（2nw?r?12nw?r?1）。

提取（tempering）輸出

MT19937 有兩個主要特性。一是周期很長，達到 219937?1219937?1，二是滿足 623623-維 3232-比特準確性。上述「旋轉」的過程，幫助我們達成了較長的周期。接下來，我們需要將每次旋轉的結果提取（tempering）輸出，以保證 MT 是 623623-維 3232-比特準確的。

提取的方法很簡單，只需要將每次旋轉得到的輸出右乘一個可逆矩陣 TT 即可。將 x? ?x? Tx→?x→T 表述成位運算，則有

y? y? y? z? ←x? ⊕(x? >>u)←y? ⊕((y? <<s)ANDb? )←y? ⊕((y? <<t)ANDc? )←y? ⊕(y? >>l)(2)(3)(4)(5)(2)y→←x→⊕(x→>>u)(3)y→←y→⊕((y→<<s)AND?b→)(4)y→←y→⊕((y→<<t)AND?c→)(5)z→←y→⊕(y→>>l)

此處，uu, ss, tt, ll 是整數參數；b? b→ 和 c? c→ 是 ww-比特的整數，用作比特遮罩（bit mask）；最終能得到的 z? z→ 即是當前輪次的隨機數輸出。

演算法描述

這樣一來，MT 演算法的主要有兩個部分：旋轉、提取。

旋轉部分有參數

ww：生成的隨機數的二進位長度；
nn：參與旋轉的隨機數個數（旋轉的深度）；
mm：參與旋轉的中間項；
rr：x? (u)kx→k(u) 和 x? (l)k+1x→k+1(l) 的切分位置；
a? a→：矩陣 AA 的最後一行。

提取部分有參數

uu, ss, tt, ll：整數參數，移位運算的移動距離；
b? b→, c? c→：比特遮罩。

於是，我們得到 MT 演算法的完整描述如下。

常量初始化。
lower→lower→upper→a? i←0←(lower→<<1)AND1for r times.←COMPLlower→←aw?1aw?2?a1a0←0lower→←0lower→←(lower→<<1)AND?1for r times.upper→←COMPL?lower→a→←aw?1aw?2?a1a0i←0
工作區非零初始化。
x? [0],x? [1],…,x? [n?1]x→[0],x→[1],…,x→[n?1]
旋轉

t? x? [i]←(x? [i]ANDupper→)OR(x? [(i+1)modn]ANDlower→)←x? [(i+m)modn]XOR(t? >>1)XOR{0? a? if t0=0otherwiset→←(x→[i]AND?upper→)OR?(x→[(i+1)modn]AND?lower→)x→[i]←x→[(i+m)modn]XOR?(t→>>1)XOR?{0→if t0=0a→otherwise
提取輸出。
y? y? y? y? y? ←x? [i]←y? XOR(y? >>u)←y? XOR((y? <<s)ANDb? )←y? XOR((y? <<t)ANDc? )←y? XOR(y? >>l)output y? y→←x→[i]y→←y→XOR?(y→>>u)y→←y→XOR?((y→<<s)AND?b→)y→←y→XOR?((y→<<t)AND?c→)y→←y→XOR?(y→>>l)output y→
更新循環變數 i←(i+1)modni←(i+1)modn，而後跳轉至步驟 3 繼續進行。

再探梅森旋轉

至此，我們已經探索了梅森旋轉演算法表面上的全部內容：我們已經知道梅森旋轉演算法是什麼，也知道梅森旋轉演算法為什麼這樣起名，也有了完整的演算法描述。但是，關於梅森旋轉演算法還有很多深層的問題我們未曾探索。比如說，對於 nn, ww 和 rr 的組合，我們是否有必要追求最長周期 PP 使得 P=wnw?r?1P=wnw?r?1？又比如說，我們提到 LFSR 取得最長周期的充要條件是反饋函數是 F2F2 上的素多項式，那麼怎樣驗證反饋函數是否是素的？

這一節，我們來討論這些問題。

關於周期

前文提到，梅森旋轉的過程，實際上是對長度為 nw?rnw?r 的二進位串做了一個 LFSR 的變體。這裡，我們將它記作 BB。

我們已經知道，這個 LFSR 的變體，其周期的上限是 2nw?r?12nw?r?1。這樣一來，整個序列的周期達到這一上限就意味著除去全零的狀態，整個序列每一種可能狀態都被遍歷了一遍；而全零的狀態則被遍歷了 1 遍。考慮在這 nw?rnw?r 比特的序列中，x? nx→n 有 n?1n?1個完整的 ww-比特向量；因此，x? nx→n 顯然是 (n?1)(n?1)-維的。這也就是說，選擇不同的隨機數種子，至多只能改變 x? nx→n 序列的起始相位。

這樣一來，我們有：當梅森旋轉達到最大周期時，若 nn 確定，n?1n?1 就確定了，進而整個序列同分布的維數 n?1n?1 也就確定了。因此，對於梅森旋轉而言，提升維數是很容易的事情。

這即是努力使梅森旋轉達到最大周期的意義。

多項式素檢測與參數調優

在梅森旋轉演算法中，反饋函數（BB 的特徵多項式）的素檢測是很容易的。這是因為，對於 p=nw?rp=nw?r（其中 pp 是梅森素數的冪）級的 LFSR 來說，其反饋函數在 F2F2 上的素檢測的複雜度是 O(p2)O(p2)。這一方面得益於梅森素數的性質，另一方面得益於 MT 是工作在 F2F2 上的演算法。(Matsumoto and Nishimura, 1997) 這一特性的證明，牽扯到很多抽象代數和數論方面的知識；此處我們按下不表，留待後續用專門的文章來證明。

梅森旋轉演算法中，要實現 PRNG 的最佳性能，需要對旋轉和提取兩部分參數做細緻的調整。調整這部分參數，尋得最優參數組合，是有特定演算法可尋的。這部分內容十分繁瑣，此處也不表。有興趣的用戶可閱讀梅森旋轉演算法原始論文第四節、第五節。

梅森旋轉演算法的 Python 實現

此處給出一個 Python 實現的梅森旋轉演算法（mt19937），為後續對演算法的「爆破」提供素材。

#! coding: utf-8nnclass MersenneTwister:n __n = 624n __m = 397n __a = 0x9908b0dfn __b = 0x9d2c5680n __c = 0xefc60000n __kInitOperand = 0x6c078965n __kMaxBits = 0xffffffffn __kUpperBits = 0x80000000n __kLowerBits = 0x7fffffffnn def __init__(self, seed = 0):n self.__register = [0] * self.__nn self.__state = 0nn self.__register[0] = seedn for i in range(1, self.__n):n prev = self.__register[i - 1]n temp = self.__kInitOperand * (prev ^ (prev >> 30)) + in self.__register[i] = temp & self.__kMaxBitsnn def __twister(self):n for i in range(self.__n):n y = (self.__register[i] & self.__kUpperBits) + n (self.__register[(i + 1) % self.__n] & self.__kLowerBits)n self.__register[i] = self.__register[(i + self.__m) % self.__n] ^ (y >> 1)n if y % 2:n self.__register[i] ^= self.__an return Nonenn def __temper(self):n if self.__state == 0:n self.__twister()nn y = self.__register[self.__state]n y = y ^ (y >> 11)n y = y ^ (y << 7) & self.__bn y = y ^ (y << 15) & self.__cn y = y ^ (y >> 18)nn self.__state = (self.__state + 1) % self.__nnn return ynn def __call__(self):n return self.__temper()nnif __name__ == "__main__":n mt = MersenneTwister(0)n tank = set()n kLen = 100n for i in range(kLen):n t = mt()n tank.add(t)n print(t)n print(len(tank) == kLen)n

爆破梅森旋轉演算法

梅森旋轉演算法的設計目的是優秀的偽隨機數發生演算法，而不是產生密碼學上安全的隨機數。從梅森旋轉演算法的結構上說，其提取演算法 __temper 完全基於二進位的按位異或；而二進位按位異或是可逆的，故而 __temper 是可逆的。這就意味著，攻擊者可以從梅森旋轉演算法的輸出，逆推出產生該輸出的內部寄存器狀態 __register[__state]。若攻擊者能夠獲得連續的至少 __n個寄存器狀態，那麼攻擊者就能預測出接下來的隨機數序列。

現在我們遵循這個思路，爆破梅森旋轉演算法。

逆向 `__temper`

我們以向右移位後異或為例，首先觀察原函數。

def right_shift_xor(value, shift):n result = valuen result ^= (result >> shift)n return resultn

簡單起見，我們觀察一個 8 位二進位數，右移 3 位後異或的過程。

value: 1101 0010nshifted: 0001 1010 # 010 (>> 3)nresult: 1100 1000n

首先，觀察到 result 的最高 shift 位與 value 的最高 shift 位是一樣的。因此，在 result 的基礎上，我們可以將其與一個二進位遮罩取與，得到 value 的最高 shift 位。這個遮罩應該是：1111 1111 << (8 - 3) = 1110 0000。於是我們得到 1100 0000。

其次，注意到對於異或運算有如下事實：a ^ b ^ b = a。依靠二進位遮罩，我們已經獲得了 value 的最高 shift 位。因此，我們也就能得到 shifted 的最高 2 * shift 位。它應該是 1100 0000 >> 3 = 0001 1000。將其與 result 取異或，則能得到 value 的最高 2 * shift 位。於是我們得到 1101 0000。

如此往複，即可復原 value。據此有代碼

def inverse_right_shift_xor(value, shift):n i, result = 0, 0n while i * shift < 32:n part_mask = ((0xffffffff << (32 - shift)) & 0xffffffff) >> (i * shift)n part = value & part_maskn value ^= part >> shiftn result |= partn i += 1n return resultn

對左移後取異或，也有類似分析。於是，得到對 __temper 的完整求逆代碼。

class TemperInverser:n __b = 0x9d2c5680n __c = 0xefc60000n __kMaxBits = 0xffffffffnn def __inverse_right_shift_xor(self, value, shift):n i, result = 0, 0n while i * shift < 32:n part_mask = ((self.__kMaxBits << (32 - shift)) & self.__kMaxBits) >> (i * shift)n part = value & part_maskn value ^= part >> shiftn result |= partn i += 1n return resultnn def __inverse_left_shift_xor(self, value, shift, mask):n i, result = 0, 0n while i * shift < 32:n part_mask = (self.__kMaxBits >> (32 - shift)) << (i * shift)n part = value & part_maskn value ^= (part << shift) & maskn result |= partn i += 1n return resultnn def __inverse_temper(self, tempered):n value = temperedn value = self.__inverse_right_shift_xor(value, 18)n value = self.__inverse_left_shift_xor(value, 15, self.__c)n value = self.__inverse_left_shift_xor(value, 7, self.__b)n value = self.__inverse_right_shift_xor(value, 11)n return valuenn def __call__(self, tempered):n return self.__inverse_temper(tempered)n

爆破

逆向 __temper() 之後，只要獲得足夠的狀態，即可構建出梅森旋轉內部的寄存器狀態。因此有如下驗證代碼。

class MersenneTwisterCracker:n __n = 624nn def __init__(self, mt_obj):n inverser = TemperInverser()n register = [inverser(mt_obj()) for i in range(self.__n)]n self.__mt = MersenneTwister(0)n self.__mt.load_register(register)nn def __call__(self):n return self.__mt()nnif __name__ == "__main__":n mt = MersenneTwister(0)n for i in range(100):n mt()n mtc = MersenneTwisterCracker(mt)n for i in range(100):n assert(mt() == mtc())n

運行後，Python 沒有拋出異常，順利推出。這說明 mtc 已能夠成功預測 mt 之後的任意順序輸出。

總結

梅森旋轉演算法，是一個優秀的偽隨機數發生演算法。在偽隨機數的評價體系中，它是一個相當優秀的演算法：周期長、均勻性好、速度快（基本都是位運算）。在條件允可的情形下，若有使用隨機數的需求，應首先考慮梅森旋轉演算法。

同時也應該注意到，梅森旋轉演算法不是為了密碼學隨機而設計的——在獲得足夠連續輸出的情況下，梅森旋轉演算法接下來的輸出值是可以準確預測的。梅森旋轉演算法容易被爆破的根源在於，其提取輸出函數是可逆的，因此暴露了其內部狀態。若要產生密碼學上的隨機數，可考慮在梅森旋轉演算法之後，拼接一值得信賴的單向雜湊函數（如 sha256）。否則，若直接用梅森旋轉演算法的輸出值作密碼學用途，則有信息泄露的風險，應引起注意。

錯誤應用梅森旋轉演算法，導致高危漏洞的一個典型是 Discuz! 的密碼重置漏洞。

擴展閱讀：梅森旋轉演算法的原始論文