乱数を巻き戻す

はじめに

(ほぼ) 最大周期を実現している擬似乱数生成器は、状態を巻き戻す（逆方向に進める）ことができます。

例えば xoshiro256 が、以下のように一見ランダムな内部状態を持っていたとしましょう。

var state = new ulong[] { 
    0x010F4C454914CD78, 
    0x83A5678480A2B416, 
    0x2652B51299006A0A, 
    0x900FEBAD58D7C533 };

この状態は、5 回の遷移 (next() 呼び出し) を経て、綺麗に並ぶようになります。

var state = new ulong[] { 
    0x0123456789abcdef, 
    0xfedcba9876543210, 
    0xdeadbeefcafebabe, 
    0x1685819840150026 };

もちろんこれは偶然に任せて求めたわけではなく、最初に下の状態があって、そこから巻き戻して上の状態を求めています。 この記事では、どうやったら巻き戻す関数を設計できるか、を紹介します。

xoshiro256 を巻き戻す

まず、xoshiro256 の内部状態 state を次に進めるメソッド Next() のコードを示します。 遷移関数はもちろん 公式の実装 と同じです。出力はここでは不要なので返り値は void にしています。 *1

// ulong[] state = new ulong[4] {...};
public void Next()
{
    ulong t = state[1] << 17;

    state[2] ^= state[0];
    state[3] ^= state[1];
    state[1] ^= state[2];
    state[0] ^= state[3];

    state[2] ^= t;
    state[3] = rotl(state[3], 45);  // rotate left.
}

public static ulong rotl(ulong x, int k) => x << k | x >> -k;

先に答えを示しておきましょう。内部状態 state を巻き戻すメソッド Prev() のコードを示します。

public void Prev()
{
    state[3] = rotl(state[3], 19);
    state[0] ^= state[3];

    ulong t = state[1] ^ state[2];
    state[1] = t ^ t << 17 ^ t << 34 ^ t << 51;
    state[3] ^= state[1];
    state[2] ^= state[0] ^ state[1] << 17;
}

どのようにして求めたか、順を追って説明していきます。

遷移を分かりやすくするため、遷移前の state を prev として保存しておき、それが Next() を通して次の state にどう反映されるかを見ます。 state の各要素について整理すると以下のようになりました。

public void Next()
{
    // copy of initial `state`
    ulong[] prev = state.ToArray();

    state[0] = prev[0] ^ prev[1] ^ prev[3];
    state[1] = prev[0] ^ prev[1] ^ prev[2];
    state[2] = prev[0] ^ prev[1] << 17 ^ prev[2];
    state[3] = rotl(prev[1] ^ prev[3], 45);
}

まず、簡単そうな state[3] の行を見ていきます。 一番下の rotl (n ビット左回転) は、合わせて 64 ビット回転させることで簡単に元に戻せます。両辺を 64 - 45 = 19 ビット左回転させることで、以下の関係式を得ます：

rotl(state[3], 19) == prev[1] ^ prev[3]    // 64 - 45 = 19

関係式の右辺を観察すると、state[0] = ... の式に同じ要素があることに気づけます。任意の x に対して x ^ x == 0 *2 が成り立つので、両式を XOR で足し合わせると：

rotl(state[3], 19) ^ state[0] == prev[0]

prev[0] を求めることができました！

さて、ここでおもむろに state[1] = ... と state[2] = ... を見ると、やはり右辺が似たような要素で構成されています。両式を XOR すると：

state[1] ^ state[2] == prev[1] ^ prev[1] << 17

この右辺が問題になります。長いので t = prev[1]; とおいて、 t ^ t << 17 から t を得る方法を模索します。一見難しそうですが、意外と簡単に戻せます。

まず、 t ^ t << 17 を 17 bit 左シフトすると ((t ^ t << 17) << 17) == (t << 17 ^ t << 34) となります。これを元の式と XOR すると、t << 17 の項が打ち消しあって 0 になるので t ^ t << 34 になります。この操作によって、シフトを 17 → 34 と大きくできました。この調子でシフトを 64 以上にできれば、64 bit の範囲から全ビットが飛び出して 0 とみなせます。

ビット列を図にするとこんな感じです。すべてを XOR することで同じ色の部分が打ち消しあい、 t すなわち prev[1] を得ることができます。

// 左辺にも同じ操作を適用します
ulong t = state[1] ^ state[2];
prev[1] == t ^ t << 17 ^ t << 34 ^ t << 51

あとは XOR のパズルで消していくだけです。 一番最初の式の両辺に、既知である prev[1] を XOR すれば prev[3] が得られます。

rotl(state[3], 19) ^ prev[1] == prev[3]

そして state[2] の式の両辺に、既知である prev[0] ^ prev[1] << 17 を XOR すれば prev[2] が得られます。

state[2] ^ prev[0] ^ prev[1] << 17 == prev[2]

以上によって、 prev を全て求めることができました。

"xor-shift" 逆算について

上の解説や図では t ^ t << k の逆算で t ^ (t << k) ^ (t << 2k) ^ (t << 3k) ^ ... とシフトを k ずつ増やしていましたが、以下のように代入を挟みながらシフトの k を 2 倍にしていくことで、小さな k に対して効率的に求められます。 *3

t ^= t << k;
t ^= t << (k << 1);
t ^= t << (k << 2);

また、逆算手法は右論理シフト t ^ t >> k に対しても、シフト方向を変えるだけで同様に適用可能です。

では右算術シフトの場合はどうでしょうか。 次式について考えてみましょう。

long t = some_value; 
long u = t ^ t >> k;

この場合、どのような t, k を選択しても u の最上位ビットは 0 になります。 t が正 (最上位ビットが 0) なら t >> k の最上位ビットも 0 、t が負 (最上位ビットが 1) なら t >> k の最上位ビットも 1 になります。これを XOR すれば、正負ともに 0 になってしまいます。

そうなると、異なる t から同一の u が得られるので、 u から t の値を確実に求めることができず、逆算はできなくなってしまいます。

ただし算術シフトが絡むとダメなのかといえばそうではなく、例えば t << k ^ t >> a は特定の k, a のペアで逆算可能になります。 shioi128 で使用されている t << 2 ^ t >> 19 などが該当します。 どのようにして逆算するかは、自分の目で確かめてみてください。 *4

PCG / 線形合同法を巻き戻す

さて次は PCG32 を巻き戻してみましょう。 ここでは出力関数を無視するので、結局のところ線形合同法と同じになります。

// ulong state = ..., ulong incr = ... | 1;
public void Next()
{
    state = state * 6364136223846793005 + incr;
}

例によって結果を先に書くと、 Prev() はこのようになります。

public void Prev()
{
    state = (state - incr) * 13877824140714322085;
}

加算 (incr) のほうは移項するだけなので分かりやすいとして、乗算のほうの定数はどこから来たのでしょうか。

実は、ulong では *5 6364136223846793005 * 13877824140714322085 == 1 となります。したがって (state * 6364136223846793005) * 13877824140714322085 == state となり、逆算できます。

では肝心の定数の求め方はというと、以下の MultiplicativeInverse() メソッドによって求めることができます。任意の奇数に対して上記の性質を満たす値を計算できます。 *6

public static ulong MultiplicativeInverse(ulong x)
{
    if ((x & 1) == 0)
        throw new ArgumentException("x must be odd");

    ulong y = x;
    y = y * (2 - y * x);
    y = y * (2 - y * x);
    y = y * (2 - y * x);
    y = y * (2 - y * x);
    y = y * (2 - y * x);

    return y;
}

メソッド内部の不思議な計算は ニュートン法 に基づいています。 時間と体力が厳しかったので解説はリンク先を参照してください。

おわりに

ここでは乱数を巻き戻す手法について紹介しました。 あまり使いどころはないかもしれませんが、実際ひとつ前の記事ではこれを調査に使っていたりします。

また、計算手法自体は出力関数の逆算にも用いることができます。例えば、Mersenne Twister の tempering (出力関数) を逆に進めて内部状態を復元するときに使えたりもします。

あえて xoshiro をチョイスしたのは、 Xorshift など有名なもの・実用されているものは先行研究が充実していて盛大に被ったためです。 それらの逆算に興味がある方は以下の記事が大変参考になります。

• TinyMT の更新関数 tinymt32_next_state(tinymt32_t*) の逆関数に関する考察 - Plus Le Blog
• Xorshift の逆関数に関する考察／正逆双方向に計算できる Xorshift - Plus Le Blog
• xoroshiro128+の状態遷移の逆算の考察 - HackMD

*1:念のため説明しておくと、C# の ulong は符号なし 64 bit 整数を表します。C の uint64_t に相当するものです。

*2:XOR のほうです、念のため

*3: を にする感じです

*4:いつか気力があったら書きます……

*5:数学っぽく言えば、 を法として

*6:偶数は何を掛けても偶数になるため、 1 になりません。したがって計算不能です。