Randomの内部seedの推移の逆算

背景

Javaにはjava.util.Randomという乱数クラスがある。
これは次のように使うことができる。

int s = 32198541;
Random r = new Random(s); // シード指定可能
int a = r.nextInt(24); // 0から23までの整数が何かしら帰る
int b = r.nextInt(24);
int c = r.nextInt(24);

さて、ここでaはsからどう求められるだろうか。
Random#nextIntは同じ値を喰わせても異なる値が帰ってくるため、Randomには内部的な変数が入っていることは明らかである。ここではその変数をseedと名付ける（内部的にも実際そういう名前である）。

コンストラクタの挙動

Randomのコンストラクタが行っていることをほぼJavaな疑似コードで示す。

long seed;

Random(long seed)
{
    this.seed = (seed ^ 0x5DEECE66DL) & 0xFFFFFFFFFFFFL;
}

実は実質たったのこれだけである。
与えられたseedになんか変な値をビットXORしたあと下位48bitだけを抜き出している。
0x5DEECE66DLは10進数では25214903917で素因数分解は7 * 443 * 739 * 11003である。

nextIntの挙動

Random#nextIntの挙動を大胆に簡略化すると大体こんな感じである。

int nextInt(int bound)
{
    seed = (seed * 0x5DEECE66DL + 0xBL) & 0xFFFFFFFFFFFFL;
    return ((int) (seed >> 17)) % bound;
}

シードの更新部分と現在のシードから求められる乱数値の出力部分に分かれている。

まとめ

これらをまとめると、最初に与えたシードとnextIntの関係はおおよそ次のようになる。

long seed = 32198541; // 与えるシード値

// 最初のシード攪乱
seed = (seed ^ 0x5DEECE66DL) & 0xFFFFFFFFFFFFL;

// nextInt
long seed2 = seed;
seed = (seed * 0x5DEECE66DL + 0xBL) & 0xFFFFFFFFFFFFL;
long seed3 = seed;
int a = ((int) (seed >> 17)) % 24;

// nextInt
seed = (seed * 0x5DEECE66DL + 0xBL) & 0xFFFFFFFFFFFFL;
int b = ((int) (seed >> 17)) % 24;

// nextInt
seed = (seed * 0x5DEECE66DL + 0xBL) & 0xFFFFFFFFFFFFL;
int c = ((int) (seed >> 17)) % 24;

問題

ここで問題である。もし上のコードのseed3が分かっているとき、seed2を簡単に求めることは可能だろうか？
多分乱数生成の理論を学べば一撃で終わる課題だがそんなもの無視して無駄に頑張る。

アプローチ

問題の整理

次の式があった時にx = f(z)を作る。

z = (x * 0x5DEECE66DL + 0xBL) & 0xFFFFFFFFFFFFL

ここで+ 0xBLは値を少し弄ればどうにかなるので取り除く。すると次の問題に簡単化できる。

z = x * 0x5DEECE66DL % 0x1000000000000L
x = f(z)

このfが分かれば1ステップ前のシードを逆算できる。

仮定1

　　0 <= x <= 0xFFFFFFFFFFFFL、
　　z = x * 0x5DEECE66DL % 0x1000000000000Lのとき、
　　xとzは全単射である。

全単射とは番号付きカードをシャッフルしたときの番号と順番の関係みたいなやつである。
この辺は数学とかに強ければ一撃で終わる問題だが、よくわからないので適当に考える。

まずz = x * 7 % 10の場合はどうだろうか。0 >= x >= 9でのx * 7 % 10は0 7 4 1 8 5 2 9 6 3となり全単射である。z = x * 5 % 10の場合では0 5 0 5 0 5 0 5 0 5で全単射でない。z = x * y % 10としたとき、xとzが全単射になるyは1 3 7 9の4つであった。

仮定2

ここで、次のことを仮定する。

　　0 <= x < m、z = x * y % mのとき、
　　yとmの最大公約数が1ならば、xとzは全単射である。

一つのxから複数のzが出ることはありえないため、複数のxが一つのzに対応するときに全単射でなくなる。z = x * 10 % 100が全単射でないのはxが10の倍数のときにz = 0でx = 0の場合とかぶるためである。
これはxが0から1ずつ上げたときに10ごとにループすると考えることができる。すなわち、z = x * 3 % 100の場合はxが値域を超えて100になったときにループが発生するが、z = x * 10 % 100のときはループが10できてしまうため全単射ではない。

x * yはどういったときにループがmより短くなる、言い換えるとx * y % m == 0に戻るのだろうか。x * yがmで割り切れるということは、x * yがmの倍数であることを意味する。そして、mのすべての因数がx * yにも含まれている場合にmの倍数となる。例えば12(= 2 * 2 * 3) * 15(= 3 * 5)は60(= 2 * 2 * 3 * 5)の倍数である。
すると、yがmの約数を1個も含んでいないとき、xはmごとにループするということになる。もしyがmの約数を1個含んでいたとき、xはその分の約数が足りなくてもyの約数と合わせて足りてしまう。

まとめると、yとmの最大公約数が1ならばxとzは全単射になる。これでとりあえず仮定2は正しいということにしておく。
すると、0x5DEECE66DL = 25214903917 = 7 * 443 * 739 * 11003なので、0x1000000000000Lとの最大公約数は1であるため仮定1も正しいということになる。

仮定3

z = x * 7 % 10の系列0 7 4 1 8 5 2 9 6 3を見ると、3倍した0 21 12 3 24 15 6 27 18 9の1の位が順番であることに気づく。すなわちx = z * 3 % 10である。

そこで、次のことを仮定する。

　　0 <= x < m、z = x * y % mのとき、
　　x = z * w % mを満たすwが存在する。

定義の対称性により、yに対する唯一の相方wについて、wに対する相方も唯一yである。

代入するとz = (z * w % m) * y % mとなる。(z * w % m) * y部分はmずつ増えても同じであるため、z * w % mのところで% mしても意味はない。なのでz = z * w * y % mと書ける。zに対して* w * y % mしてもzであるということは多分w * y % m == 1を意味する。

確かめるためにとりあえずスクリプトでm = 100, 64のときのy:wを列挙してみた。

m=100

1:1|3:67|7:43|9:89|11:91|13:77|17:53|19:79|21:81|23:87|27:63|29:69
31:71|33:97|37:73|39:59|41:61|43:7|47:83|49:49|51:51|53:17|57:93|59:39
61:41|63:27|67:3|69:29|71:31|73:37|77:13|79:19|81:21|83:47|87:23|89:9
91:11|93:57|97:33|99:99

m=64

1:1|3:43|5:13|7:55|9:57|11:35|13:5|15:47|17:49|19:27
21:61|23:39|25:41|27:19|29:53|31:31|33:33|35:11|37:45|39:23
41:25|43:3|45:37|47:15|49:17|51:59|53:29|55:7|57:9|59:51
61:21|63:63

1:1|3:2B|5:D|7:37|9:39|B:23|D:5|F:2F
11:31|13:1B|15:3D|17:27|19:29|1B:13|1D:35|1F:1F
21:21|23:B|25:2D|27:17|29:19|2B:3|2D:25|2F:F
31:11|33:3B|35:1D|37:7|39:9|3B:33|3D:15|3F:3F

y:wのときにw:yであることが確認できる。また、y * w % mを計算するとすべて1であった。そのためm = 100のときはyとwの1の位をかけると必ず1の位が1になる。そのため1の位は必ず1:1|3:7|7:3|9:9と対応している。九九表の中で1の位が1である場所はそこしかないのだ。

そして、次のようにちゃんとxを求める関数として使うことができる。

52 = 76 * 27 % 100
76 = 52 * 63 % 100

27 * 63 = 1701

もっと大きな数でも同じ計算ができる。

3596 = 7852 * 273 % 10000
7852 = 3596 * 6337 % 10000

273 * 6337 = 1730001

仮定3はとりあえず正しいことにしておく。おそらくy * w % m == 1となるwを探すことがx = f(z)を計算する十分条件である。

桁ごとの分割

さて、ここで0x5DEECE66DLの相方をGPGPUごり押しで探してもいいのだが、先ほど「m = 100のときに1の位が固定」という面白い性質を確認した。これは16進数でも言えるはずである。
実際、m = 64のときに1の位には1:1|3:B|5:D|7:7|9:9|B:3|D:5|F:Fという対応が存在する。この組み合わせでしか0001という2進数を下位4ビットに発生させることができないのであろう。実際掛け合わせると次のように16で割った余りが1になる。

0x1 * 0x1 =   1 =  0 * 16 + 1
0x3 * 0xB =  33 =  2 * 16 + 1
0x7 * 0x7 =  49 =  3 * 16 + 1
0x5 * 0xD =  65 =  4 * 16 + 1
0x9 * 0x9 =  81 =  5 * 16 + 1
0xF * 0xF = 225 = 14 * 16 + 1

ここから言えることは、0x5DEECE66DLの相方は1の位が0x5ということである。これだけで1/16の計算量で済む。だがm = 64のときに16進数表記でできて256進数表記でできない道理は見つからない。恐らく16進数で2桁ごとに区切っても良いはずだ。

そこで、z = x * 0x6D % 0x100の相方wを探してみる。相方は0x65で0x6D * 0x65 = 0x2B01である。これで0x5DEECE66DLの相方は下位2桁が0x65であるという予想ができた。

z = x * 0xE66D % 0x10000の相方は0x1365で、かけると0x11750001になる。この計算は相方の候補を順番に生成して* 0xE66D % 0x10000して1になるかどうかを確認すればよいが、相方の下2桁が0x65であることはわかっているため、残りの2桁だけをループさせればよいのだ。

z = x * 0x5DEECE66DL % 0x1000000000000Lの相方は0xDFE05BCB1365Lであり、0x5DEECE66DL * 0xDFE05BCB1365L = 0x86E9000000000001である。これで乗算部分の逆関数を求めることができた。

z = x * 0x5DEECE66DL & 0xFFFFFFFFFFFFL
x = z * 0xDFE05BCB1365L & 0xFFFFFFFFFFFFL

あとは+ 0xBL部分を少し追加してやれば完了である。

z = (x * 0x5DEECE66DL + 0xBL) & 0xFFFFFFFFFFFFL
x = (z - 0xBL) * 0xDFE05BCB1365L & 0xFFFFFFFFFFFFL

結論

Randomのシードは次のように逆算できる。

long seed = 0x32522523;
seed = (seed * 0x5DEECE66DL + 0xBL) & 0xFFFFFFFFFFFFL;
System.out.println(String.format("%x", seed));
seed = (seed - 0xBL) * 0xDFE05BCB1365L & 0xFFFFFFFFFFFFL;
System.out.println(String.format("%x", seed));

結果

86e8d09b41f2
32522523