64bit LCGの検索

Pokémon RNG Advent Calendar 2017の19日の記事です。

第5世代のポケモンでは次の64bitのLCG(線形合同法)が使われています。

X_{n+1} = (X_n * 0x5d588b656c078965 + 0x269ec3) % 2^64

64bitとなると全seedをしらみつぶすということが現実的な時間では不可能です。

16bitの乱数値がわかる状況を考えます。この場合でも調べるseedの個数は2^48個あってやはり全探索は時間がかかりそうです。(2^32個の探索が1秒で済むと仮定しても18時間かかる)

しかし、実は16bitの乱数情報があればそれなりに速い時間で(たとえば私のパソコンではシングルスレッドでも19秒で)条件を満たすseedをすべて求められます。
正確には乱数値16bitとその10個先の乱数値16bitです。
すなわち

f(x) = (x * A + B) % 2^64
g(x, n) = ((x >> 32) * n) >> 32
where A = 0x5d588b656c078965, B = 0x269ec3

とおいたとき与えられたa, bに対して

g(x, 65536) = a
g(f^10(x), 65536) = b

を満たすxをすべて求められます。

仕組みはこうです。

xstep := 0x05114191

という定数を定義します。

この定数により次の式は

(xstep * A^10) % 2^64 = 0xbf4a64c9 =: ystep

という小さい値になります。 (そうなるようなxstepと10という定数を探したのです)

するとxをxstepずつ大きくしたときにy = f^10(x)はystepずつ大きくなることより、yの値が目標の値に到達するところまでスキップできることになります。

実際xをxstepずつ大きくすると以下のようになります。

i = 00, x = 0000000000000000, y = 67795501267f125a
i = 01, x = 0000000005114191, y = 67795501e5c97723
i = 02, x = 000000000a228322, y = 67795502a513dbec
i = 03, x = 000000000f33c4b3, y = 67795503645e40b5
i = 04, x = 0000000014450644, y = 6779550423a8a57e
i = 05, x = 00000000195647d5, y = 67795504e2f30a47
i = 06, x = 000000001e678966, y = 67795505a23d6f10
i = 07, x = 000000002378caf7, y = 677955066187d3d9
i = 08, x = 00000000288a0c88, y = 6779550720d238a2
i = 09, x = 000000002d9b4e19, y = 67795507e01c9d6b
i = 10, x = 0000000032ac8faa, y = 677955089f670234
i = 11, x = 0000000037bdd13b, y = 677955095eb166fd
i = 12, x = 000000003ccf12cc, y = 6779550a1dfbcbc6
i = 13, x = 0000000041e0545d, y = 6779550add46308f
i = 14, x = 0000000046f195ee, y = 6779550b9c909558
i = 15, x = 000000004c02d77f, y = 6779550c5bdafa21

yがちょっとずつしか進まないので目的の値の範囲にくるまでがばっとスキップできることがわかると思います。

このアイディアはCracking Pseudorandom Sequences Generators in Java Applicationsからいただきました。

さて、これを実装します。
0個先と10個先の乱数値だけだと全部出力するとかなりの個数(2^32くらい)になるので1個先と2個先の乱数値も条件に加えることにします。

#include <cstdio>
#include <cstdint>
#include <cmath>
#include <algorithm>
#include <chrono>
#include <random>

typedef int64_t s64;
typedef uint64_t u64;
typedef __int128 s128;

const u64 A = 0x5d588b656c078965ll;
const u64 B = 0x269ec3ll;
const s128 M = (s128)1 << 64;
const int BITS = 64;
const int P = 16;

class LCGOperator {
public:
    u64 a, b;
    LCGOperator(u64 aa, u64 bb) : a(aa), b(bb) {}

    LCGOperator o(LCGOperator y) {
        LCGOperator x = *this;
        LCGOperator res(x.a * y.a, x.a * y.b + x.b);
        return res;
    }

    LCGOperator pow(u64 n) {
        LCGOperator x = *this;
        if (n == 0) {
            return LCGOperator(1, 0);
        } else if (n % 2 == 1) {
            return x.o(x).pow(n / 2).o(x);
        } else {
            return x.o(x).pow(n / 2);
        }
    }

    u64 apply(u64 s) {
        return a * s + b;
    }
};

void search(u64 s0, u64 s1, u64 s2, u64 s10, bool check) {
    LCGOperator op0(A, B);
    LCGOperator op = op0.pow(10);
    const u64 xstep = 0x05114191;
    const u64 ystep = (xstep * op.a) % M;

    const s128 xstart = (s128)s0 << (BITS - P);
    const s128 xend = (s128)(s0 + 1) << (BITS - P);
    const s128 ystart = (s128)s10 << (BITS - P);
    const s128 yend = (s128)(s10 + 1) << (BITS - P);
    u64 found_count = 0;

    for (u64 i = 0; i < xstep; i ++) {
        s128 x = xstart + i;
        s128 y = op.apply(x);
        s128 s;

        while (x < xend) {
            while (ystart <= y && y < yend && x < xend) {
                if ((y % M) >> (BITS - P) == s10 && (!check || ((op0.apply(x) >> (BITS - P)) == s1 && (op0.pow(2).apply(x) >> (BITS - P)) == s2))) {
                    if (check) printf("found %016lx\n", (u64)x);
                    found_count ++;
                }
                x += xstep;
                y += ystep;
            }

            s128 ynext = (y >= ystart ? M : 0) + ystart - y;
            s = ynext / ystep + (ynext % ystep != 0 ? 1 : 0);

            y = (y + ystep * s) % M;
            x += xstep * s;
        }
    }
    printf("found_count=%ld\n", found_count);
}

int main() {
    std::mt19937 rnd;
    for (int i = 0; i < 10; i ++) {
        u64 seed = (u64)rnd() << 32 | rnd();
        printf("seed=%016lx\n", seed);
        LCGOperator op0(A, B);
        LCGOperator op = op0.pow(10);
        u64 s0 = seed >> (BITS - P);
        u64 s1 = op0.apply(seed) >> (BITS - P);
        u64 s2 = op0.pow(2).apply(seed) >> (BITS - P);
        u64 s10 = op.apply(seed) >> (BITS - P);
        auto start = std::chrono::system_clock::now();
        search(s0, s1, s2, s10, true);
        auto end = std::chrono::system_clock::now();
        printf("time=%ld\n", std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count());
    }
    return 0;
}

このプログラムは解の一つとなるべきseedを一つ乱数 (MT使用)で決めて、探索します。ちゃんと決めた解が探索によって見つかっていればいいわけです。

実行結果は次のようになります。

seed=d091bb5c22ae9ef6
found d091bb5c22ae9ef6
found d091de63e71f9f06
found d09117bbd62d8f67
found_count=3
seed=e7e1faeed5c31f79
found e7e1faeed5c31f79
found_count=1
seed=2082352cf807b7df
found 2082352cf807b7df
found_count=1
seed=e9d300053895afe1
found e9d3c6ad4987bf80
found e9d300053895afe1
found e9d3230cfd06aff1
found_count=3
seed=a1e24bba4ee4092b
found a1e24bba4ee4092b
found_count=1
...

ちゃんと見つかっていますね!

また、0個先と10個先の乱数値だけを条件にして探索して得られた個数がLattE (線形不等式を満たす整数点の個数を出力するソフトウェア)によって数えられた個数と一致することは何個かの具体例について確認済みです。

で、実際問題、第5世代のポケモンで0個先と10個先の乱数の16bit値が観測できる状況があるのかというとないです。PWT (ポケモンワールドトーナメント)でレンタルポケモンのトレーナーIDが乱数で決まっていればチャンスがあったのになという感じです。

16bitじゃなくて6bitか7bitくらいの乱数値からseedを現実的な時間で見つけるアルゴリズムを見つけられれば、ペラップの音程からseedを調べられますね。興味のある方、ぜひチャレンジしてみてください。

筆者: oupo (連絡先: oupo.nejiki@gmail.com)