Hatena::ブログ(Diary)

oupoの日記

2016-12-13

63変数連立2次方程式が解ければTinyMTのstateは連続した64個の乱数値の下位2bitから得られる

00:01

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は次の連立方程式を解けば得られます。

f:id:oupo:20161212215141p:image

あとは、この方程式グレブナー基底を求めるなり、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))

得られた連立方程式は以下にアップロードしてあります。

トラックバック - http://d.hatena.ne.jp/oupo/20161213/1481554861
Connection: close