63変数連立2次方程式が解ければTinyMTのstateは連続した64個の乱数値の下位2bitから得られる
Pokémon RNG Advent Calendar 2016 13日目の記事です。
疑似乱数TinyMTにおいて得られた乱数値から乱数生成器の状態を逆算する話です。
SM孵化乱数列を計算する : ただの雑記byさきによってTinyMTのstateは連続した127個の乱数値の下位1bitがあれば得られることが示されました。
これを連続した64個の乱数値の下位2bitから得られるともっと楽できそうです。
これは63変数63本のF_2係数連立2次方程式が解ければ可能です。
たとえば、連続する64個の下位2bitが01, 00, 00, …, 00であるstateは次の連立方程式を解けば得られます。
あとは、この方程式のグレブナー基底を求めるなり、SAT-based Constraint Solverにかけるなりして解けばいいですね！
どちらも試しましたが、変数の個数の指数のオーダーの時間がかかって63変数は無理そうでした！
ソースコードは以下の通り。SageMathCloudで動かすことができます。
mat1 = 0x8f7011ee mat2 = 0xfc78ff1f tmat = 0x3793fdff F2 = GF(2) # 整数を32次元F2ベクトルに変換する def int_to_f2(x): return vector([F2((x >> i) & 1) for i in range(32)]) # 整数列をF2ベクトルに変換する def ints_to_f2(xs): return vector([F2((xs[floor(i / 32)] >> (i % 32)) & 1) for i in range(32*len(xs))]) # 32次元F2ベクトルを整数に変換する def f2_to_int(vec): x = 0 for i in range(32): x |= Integer(vec[i]) << i return x # F2ベクトルを整数列に変換する def f2_to_ints(vec): return [f2_to_int(vec[i*32:i*32+32]) for i in range(floor(len(vec) / 32))] # 行列を横につなげる def join_matrix_horizontal(mats): m = mats[0].nrows() n = mats[0].ncols() return matrix([[mats[floor(j / n)][i, j % n] for j in range(n*len(mats))] for i in range(m)]) # 行列を縦につなげる def join_matrix_vertical(mats): m = mats[0].nrows() n = mats[0].ncols() return matrix([[mats[floor(i / m)][i % m, j] for j in range(n)] for i in range(m*len(mats))]) # 行列からそのi番目の行だけを得る def nth_row(mat, i): m = mat.nrows() n = mat.ncols() return matrix([[mat[i, j] for j in range(n)]]) # 128次元ベクトルを32次元ベクトルを4つつなげたものと思ってそのi番目をとってくるという変換の行列 def elem_matrix(i): mats = [Mat(GF(2), 32, 32).zero() for j in range(4)] mats[i] = Mat(GF(2), 32, 32).identity_matrix() return join_matrix_horizontal(mats) # 真なら1、偽なら0を返す def xi(x): if x: return 1 else: return 0 # 32次元ベクトルに対して右シフトを表す行列 def rshift_matrix(k): return matrix([[F2(xi(j - i == k)) for j in range(32)] for i in range(32)]) # 32次元ベクトルに対して左シフトを表す行列 def lshift_matrix(k): return rshift_matrix(-k) # 32次元ベクトルに対して第kビットを得るという変換の行列 def get_bit_matrix(k): return matrix([[F2(xi(j == k)) for j in range(32)]]) get_0bit = get_bit_matrix(0) # n次元ベクトルをn×1行列に変換する def vect_to_matrix(vec): return matrix([[vec[i]] for i in range(len(vec))]) # 32次元ベクトルに対してaとandをとるという変換を表す行列 def mask_matrix(a): return matrix([[F2(a[i] * xi(i == j)) for j in range(32)] for i in range(32)]) # 128次元ベクトルから調律した32次元ベクトルを返す変換を表す行列 def temper_matrix(): t0 = elem_matrix(3) t1 = elem_matrix(0) t2 = rshift_matrix(8) * elem_matrix(2) return t0 + (t1 + t2) + vect_to_matrix(int_to_f2(tmat)) * get_0bit * (t1 + t2) # 128次元ベクトルから調律に使われるt1の値を返す変換を表す行列 def temper_t1_matrix(): return elem_matrix(0) # 128次元ベクトルから調律に使われるt2の値を返す変換を表す行列 def temper_t2_matrix(): return rshift_matrix(8) * elem_matrix(2) temper = temper_matrix() temper_t1 = temper_t1_matrix() temper_t2 = temper_t2_matrix() # 128次元のstateから次のstateを得る変換を表す行列 def next_state_matrix(): y = elem_matrix(3) x = mask_matrix(int_to_f2(0x7fffffff)) * elem_matrix(0) + elem_matrix(1) + elem_matrix(2) x += lshift_matrix(1) * x y += rshift_matrix(1) * y + x st0 = elem_matrix(1) st1 = elem_matrix(2) + vect_to_matrix(int_to_f2(mat1)) * get_0bit * y st2 = x + lshift_matrix(10) * y + vect_to_matrix(int_to_f2(mat2)) * get_0bit * y st3 = y return join_matrix_vertical([st0, st1, st2, st3]) next_state = next_state_matrix() # 128次元ベクトルを4つの32次元ベクトルをつなげたものとみたとき第0要素の最上位bitを無視して127bitを得るという変換の行列 def from_128bit_to_127bit_matrix(): m = [[F2(0) for j in range(128)] for i in range(127)] for i in range(31): m[i][i] = F2(1) for i in range(31, 127): m[i][i+1] = F2(1) return matrix(m) # 逆に127bitから1bit水増しして128次元ベクトルを返すという変換を表す行列 def from_127bit_to_128bit_matrix(): m = [[F2(0) for j in range(127)] for i in range(128)] for i in range(31): m[i][i] = F2(1) for i in range(31, 127): m[i+1][i] = F2(1) return matrix(m) from_128bit_to_127bit = from_128bit_to_127bit_matrix() from_127bit_to_128bit = from_127bit_to_128bit_matrix() # 127次元のstateから次のstateを得る変換を表す行列 next_state_127 = from_128bit_to_127bit * next_state * from_127bit_to_128bit # 127次元のstateから連続する127個の乱数のそれぞれ最下位bitを得る変換の行列 def from_state_to_random_stream_matrix(): mats = [None for i in range(127)] step = Mat(GF(2), 127, 127).identity_matrix() A = get_0bit * temper * from_127bit_to_128bit for i in range(127): mats[i] = A * step step *= next_state_127 return join_matrix_vertical(mats) # stateから調律後の第1bitの連続した127個分を得る def from_state_to_tempered_bit1_stream_matrix(): mats = [None for i in range(127)] step = Mat(GF(2), 127, 127).identity_matrix() A = get_bit_matrix(1) * temper * from_127bit_to_128bit for i in range(127): mats[i] = A * step step *= next_state_127 return join_matrix_vertical(mats) # stateから調律に使われるt1の第0bitの連続した127個分を得る def from_state_to_t1_stream_matrix(): mats = [None for i in range(127)] step = Mat(GF(2), 127, 127).identity_matrix() A = get_0bit * temper_t1 * from_127bit_to_128bit for i in range(127): mats[i] = A * step step *= next_state_127 return join_matrix_vertical(mats) # stateから調律に使われるt2の第0bitの連続した127個分を得る def from_state_to_t2_stream_matrix(): mats = [None for i in range(127)] step = Mat(GF(2), 127, 127).identity_matrix() A = get_0bit * temper_t2 * from_127bit_to_128bit for i in range(127): mats[i] = A * step step *= next_state_127 return join_matrix_vertical(mats) from_state_to_random_stream = from_state_to_random_stream_matrix() # 127個の連続した第0bitからstateを得る行列 （逆行列を計算） from_random_stream_to_state = from_state_to_random_stream.inverse() from_state_to_tempered_bit1_stream = from_state_to_tempered_bit1_stream_matrix() from_state_to_t1_stream = from_state_to_t1_stream_matrix() from_state_to_t2_stream = from_state_to_t2_stream_matrix() # 連立一次方程式の解のパラメータから解を得る変換の行列 def from_param_to_state_matrix(solution0, basis): m = [[0 for j in range(len(basis)+1)] for i in range(127)] for i in range(127): for j in range(len(basis)+1): if j == 0: m[i][j] = solution0[i] else: m[i][j] = basis[j-1][i] return matrix(m) # 係数がベクトルaの成分の一次式を文字列で返す def poly(a): poly = [] for i in range(len(a)): if a[i] != 0: if i == 0: poly.append("1") else: poly.append("x"+str(i-1)) if poly != []: return "+".join(poly) else: return "0" # a, b, cから決まる方程式を文字列で返す def equation(a, b, c): return poly(a) + "+(" + poly(b) + ")*(" + poly(c) + ")=0" # 乱数列の最下位bitが1,0,…,0,…(最初に1、続いて0が63個)であるようなstateを求めるという一次方程式を解き、解空間(solution1, basis)を得る I = matrix([[xi(i == j) for j in range(127)] for i in range(64)]) vec = vector([xi(i == 0) for i in range(64)]) solution0 = (I * from_random_stream_to_state).solve_right(vec) basis = kernel(((I * from_random_stream_to_state).transpose())).basis() from_param_to_state = from_param_to_state_matrix(solution1, basis) for i in range(len(basis)): a = list(nth_row(from_state_to_tempered_bit1_stream, i) * from_param_to_state)[0] b = list(nth_row(from_state_to_t1_stream, i) * from_param_to_state)[0] c = list(nth_row(from_state_to_t2_stream, i) * from_param_to_state)[0] print(equation(a, b, c))
得られた連立方程式は以下にアップロードしてあります。
