1年以上昔にこんな記事を書いた。

久しぶりに読み返してみると、疑似乱数が決定的・周期的なことまではなんとなく示せているものの具体的な予測方法は全然書いてないじゃん！！！と内なる声に突っ込まれた。「予測できるか試したい」なんてタイトル掲げておいて、なんて中途半端な奴なんだ。

多分このときは周期性を実験するところで気力が尽きたのかもしれない、許してあげてね。

そういうわけで今回こそは実際に予測してみようと思う。

xorshift 32bit

今回は python を使いたかったので前の記事で使った疑似乱数生成器を移植した。

また、今回出力は self.s をそのまま使わずに下位 16bit だけを取り出す形にした。これは「出力 = 内部状態」だと予測する必要がなくなってしまうためである。¹

ちなみに内部状態と出力を異なる bit 数にする効果については以下の記事で昔考えている。

疑似乱数（xorshift）の内部状態と生成される数値の関係について考えてみる - Buri Memo:

class XorShift():
    def __init__(self, seed):
        self.s = seed

    # python の int 型は無制限なので 32bit に切り詰める
    def mask(self, x):
        return x & 0xffffffff

    def rand(self):
        # https://ja.wikipedia.org/wiki/Xorshift の xorshift32
        self.s = self.mask(self.s ^ (self.s << 13))
        self.s = self.mask(self.s ^ (self.s >> 17))
        self.s = self.mask(self.s ^ (self.s << 5))
        
        # 出力は内部状態変数 s の下位 16bit
        return self.s & 0xffff

これを実行すると、以下のように 0~65535 の乱数を得られる。今回はこの 10個の値から次の乱数出力（7649）を予想できるか試してみたいと思う。

sequence = [rng.rand() for i in range(10)]

# [17433, 63802, 48521, 8888, 60923, 59364, 53581, 56036, 61202, 34977]
print("Randoms", sequence)

# 7649 （これを予想する）
print(f"Next random is {rng.rand()}")

総当りで予測する

一番簡単に思いつく方法はこれ。適当なシード値から [17433, 63802, 48521, 8888, 60923, 59364, 53581, 56036, 61202, 34977] の並びが出現するまで乱数を作りまくるというシンプルな方法。

この乱数生成器の周期は 2³²-1 であるので平均的には 21億回ほどで見つかるはず。実際にやってみると 50分くらいかかったが次の出力 7649 を予測することができた。

st = time.time()
rng = XorShift(1)
samples = [rng.rand() for i in range(len(sequence))]
while samples != sequence:
    samples = [*samples[1:], rng.rand()]
print(f"Predict: {rng.rand()} ({time.time() - st}s)")
# Predict: 7649 (2990.971052646637s)

流石に時間かかりすぎなので C言語でもやってみると 5秒くらいで終わった。やっぱ python は総当りに向いてないなぁ。
ところで Math.random() で使われる xorshift 128+ は周期が 2¹²⁸ - 1 なので、平均でその半分の 8 * 10²⁸ 倍もかかることになる。現実的な時間で終わるのは厳しそう。

blog-assets/xorshift32-prediction/main.c at main · buri83/blog-assets · GitHub

Z3 を使用して予測する

Z3 は Microsoft Research の SMT solver らしい。正直 SMT について全く理解できてないが、とりあえずは与えられた論理式（制約）を満たす解を1個高速に探してくれるソフトウェアというイメージでいいかなと思う。

github.com

例えば z3 を使って以下のような x, y を求めてみる。

x, y は自然数
x + 2y < 5 を満たす

以下のスクリプトを実行すると [x = 1, y = 1] [x = 2, y = 1] と正しい答えが得られた。スゲ〜〜

import z3

x = z3.Int("x")
y = z3.Int("y")

solver = z3.Solver()
solver.add(x > 0, y > 0)
solver.add(x + 2 * y < 5)

# [x = 1, y = 1] [x = 2, y = 1]
while solver.check() == z3.sat:
    answer = solver.model()
    print(answer)
    # 解を1個ずつしか表示できないので、見つけた解を除外してもう一度解いてもらう
    solver.add(z3.Not(z3.And(x == answer[x], y == answer[y])))

乱数生成器もただの数式なので z3 を使って同じように計算できる。（そういえば逆方向の変換できないことを目的にしているハッシュ関数はどうなるんだろう？）

実行すると一瞬で seed を求めることができた！（たったの 63 ms ！！）

import z3

st = time.time()
solver = z3.Solver()

s = z3.BitVec("s", 32)
for i in sequence:
    s = s ^ (s << 13)
    s = s ^ z3.LShR(s, 17)
    s = s ^ (s << 5)
    solver.add(s & 0xffff == i)

if solver.check() != z3.sat:
    raise Exception("unsat !!!")

model = solver.model()
states = {str(s): model[s] for s in model}

seed = states["s"].as_long()

rng = XorShift(seed)
for _ in range(len(sequence)):
    rng.rand()

print(f"Predict: {rng.rand()} ({time.time() - st}s)")
print("Predicted initial seed", seed)
# Predicted initial seed 777
# Predict: 7649 (0.06324553489685059s)

ちなみにこの方のスクリプトや解説動画をかなり参考にさせてもらった。こちらでは Node.js の Math.random() を実際に予想しており周期が 128bit にもかかわらず一瞬で求められている。また、Math.random() の出力は少数の値なので float を vector に変換するような処理もありより実用的な実装もされていたりもする。

気になる方は是非見てほしい。

github.com