Randomの内部seedの推移の逆算


# 背景

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

```java:
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な疑似コードで示す。

```java
long seed;

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

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

## `nextInt`の挙動

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

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

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

## まとめ

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

```java
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)`を作る。

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

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

```java
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`を求める関数として使うことができる。

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

27 * 63 = 1701
```

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

```java
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になる。

```java
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`である。これで乗算部分の逆関数を求めることができた。

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

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

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

# 結論

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

```java
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
```