出力関数が線形な場合のTinyMTにおいて連続する4つの乱数値からstateを得る

疑似乱数TinyMTの逆算について考えます。https://odanpoyo.github.io/2016/12/10/Tiny-Mersenne-twister%E3%81%AE%E9%80%86%E7%AE%97/の続きともいえるような内容です。
上の記事はstateから一つ前のstateを計算するという話でしたが、この記事では連続する4つの乱数値からstateを得ることを考えます。

TinyMTは何も設定しなければ出力関数が非線形ですが、LINEARITY_CHECKをONにすると線形になります。その場合、連続する4つの乱数値からstateを得るのはただの（F_2係数の）連立一次方程式を解くだけなので簡単にできます。

SageMath Cloudで動かせるプログラムを作りました。
動かす場合、メニューのKernel→Change Kernel→SageMath 7.4と選ぶ必要があるかもしれません。

なお、TinyMTの実装はTinyMT/tinymt32.h at master · MersenneTwister-Lab/TinyMT · GitHubのtinymt32_next_stateとtinymt32_temperを参照。

LINEARITY_CHECKがOFFの場合、つまり出力関数が非線形な場合は乱数値からstateを得るのは現実的な実行時間じゃ無理だろうなと思っていました。しかしそれをさき (@water_blow)さんがくつがえしてくれました。

そうなんです。LINEARITY_CHECKがOFFのとき以下のように算術加算が使われていて非線形なので逆算は不可能に見えます。しかし下位1bitを見ればxorしていることと変わらないので上のツイートのように逆算が可能になります。

	inline static uint32_t tinymt32_temper(tinymt32_t * random) {
		uint32_t t0, t1;
		t0 = random->status[3];
#if defined(LINEARITY_CHECK)
		t1 = random->status[0]
			^ (random->status[2] >> TINYMT32_SH8);
#else
		t1 = random->status[0]
			+ (random->status[2] >> TINYMT32_SH8);
#endif
		t0 ^= t1;
		t0 ^= -((int32_t)(t1 & 1)) & random->tmat;
		return t0;
	}

全く気が付きませんでした。すごいですね。さきさんがこの件について記事を書いてくれることを期待してみます。

追記(2016/12/11)
さきさんが書いてくれました。

追記(2016/12/11)
こちらでも連続した127個の乱数の下位1bitからstateを算出する行列を計算しました。
さきさんのツールと結果が一致します。

32bit LCGのseedを手計算で求めよう (半分ネタ)


X_{n+1} = (a X_n + b) \bmod 2^{32}
　　　　 a = \mathrm{0x41c64e6d},\;b = \mathrm{0x6073}

で与えられる線形合同法疑似乱数において、その上位16bit

Y_n = \lfloor X_n / 2^{16} \rfloor

が観測できるという状況を考えます。
そのとき $X_0$ の値を求めるにはどうすればよいでしょうか。

$Y_0$ を $\alpha$ とおくとき $X_0$ は


X_0 = 2^{16} \alpha + x,\;0 \le x < 2^{16}

とかけます。よって $2^{16}$ 通りの $x$ すべてをしらみつぶして、 $Y_1$ の値が一致するものを調べればよさそうです。$Y_1$ の値が一致するものは複数個ある可能性がある場合があるのでその場合は $Y_2$ の値を調べます。これで疑似乱数のseedを特定できるでしょう。

しかし、$2^{16}$ 通り調べるのは手計算ではできないですね。
そこで手計算できる方法を考えます。

まず、


X_{n+2^{16}} = (\mathrm{0xdfa40001} X_n + \mathrm{0x47310000}) \bmod 2^{32}

であることに注意します。
特に係数がそれぞれ $2^{16}$ の倍数に1足した値と$2^{16}$ の倍数になっているので、それらを$2^{16} A + 1$, $2^{16} B$とおいてみます。

すると上の $X_0$ の候補の値について


\begin{align*} X_{2^{16}} &= (2^{16} A + 1) (2^{16} \alpha + x) + 2^{16} B \\ &= 2^{16}(Ax+\alpha+B) + x\end{align}

が成り立ちます。ここで $Y_{2^{16}}$ を $\beta$ とおけば

Ax+\alpha+B \equiv \beta \pmod{2^{16}}

なのでこれを解くことにより $x$ を得ます！
ただし $A$ は4の倍数になっているので、得られるのは $x \bmod 2^{14}$ です。したがってこの方法で得られるseedは4通りあるので、その後 $Y_1$ などを使って一通りに絞らないといけません。

上の方法をまとめると


1) 乱数を1つ生成し、その値(16bit)を得る ($\alpha$ とする)
2) 乱数を2^{16}-1 = 65535個進める
3) 乱数を1つ生成し、その値(16bit)を得る ($\beta$ とする)
4) $\alpha$ と $\beta$ からseedが4通りに絞れる
となります。

さて、よく見ると 0xdfa40001 は 2^{16} の倍数 + 1 であるだけでなく、 2^{18} の倍数 + 1 です。ここに気付くと上の方法は改良できます。
2^{16}個先を考えるのではなく2^{14}個先を考えるのです。


X_{n+2^{14}} = (\mathrm{0xf7e90001} X_n + \mathrm{0xf1cc4000}) \bmod 2^{32}

ここに出てくる係数をあらためて $2^{16} A + 1 = 0xf7e90001, 2^{14} B = 0xf1cc4000$ とおきます。

すると $X_0$ の候補の値について


\begin{align*} X_{2^{14}} &= (2^{16} A + 1) (2^{16} \alpha + x) + 2^{14} B \\ &= 2^{14}(4Ax+4\alpha+B) + x\end{align}

が成り立ちます。ここで $X_{2^{14}}$ を $2^{16} \beta + y$ とおけば

\begin{align*} 2^{14}(4Ax+4\alpha+B)+x = 2^{16} \beta + y.\end{align}

$x = 2^{14} x' + x'', y = 2^{14} y' + y'', 0 \ge x', y' < 4, 0 \ge x'', y'' < 2^{14}$ とおけば

\begin{align*} 2^{14}(4Ax+4\alpha+B+x')+x'' = 2^{14} (4 \beta + y') + y''.\end{align}

上位18bitをとることにより

4Ax+4\alpha+B + x' \equiv 4 \beta + y' \pmod{2^{18}}

$y' - x' = \delta$とけば

4Ax+4\alpha+B \equiv 4 \beta + \delta \pmod{2^{18}} (1)
　　　　 -4 < \delta < 4

この式に適する$\delta$を選び、$x$について解くと解の候補が得られます。$x$と$\delta$の整合性をチェックし、解になっているか確かめられます。

上の方法をまとめると


1) 乱数を1つ生成し、その値(16bit)を得る ($\alpha$ とする)
2) 乱数を2^{14}-1 = 16383個進める
3) 乱数を1つ生成し、その値(16bit)を得る ($\beta$ とする)
4) $\alpha$ と $\beta$ からseedが最大2通り(最小0通り)に決定される
となります。

試しにやってみましょう。$\alpha = 11111, \beta = 22222$としてみます。
式(1)は

253860x+44444+247601 \equiv 88888 + \delta \pmod{2^{18}}

です。整理すると
253860x \equiv -203157 + \delta \pmod{2^{18}}.

$-203157 + \delta$は4の倍数にならざるを得ないことから$\delta = 1, -3$です。

$\delta = 1$のとき：

253860x \equiv -203156 \pmod{2^{18}}.

両辺を4で割り、
63465x \equiv -50789 \pmod{2^{16}}.

63465の2^{16}を法とした逆元はユークリッドの互除法により20569と求まるので
x \equiv 20569 \cdot (-50789) \pmod{2^{16}}.

計算すると
x \equiv 30435 \pmod{2^{16}}

と求まります。
このとき$x' = 1$なので$\delta = 1$と整合します （$\delta = y' - x'$をみたす$0 \ge y' < 4$があります）。

$\delta = -3$のとき：
同様の計算で

x \equiv 51004 \pmod{2^{16}}

と求まります。
このとき$x' = 3$なので$\delta = 1$と整合しません （$\delta = y' - x'$をみたす$0 \ge y' < 4$がありません）。

よって解は30435の一つです。
ゆえにseedは 2^{16} \cdot 11111 + 30435 = 728200931 = \mathrm{0x2b6776e3}です。

2^{14}回乱数を進めるなんて2^{16}通り機械でしらみつぶしするよりずっと大変じゃないかということで、半分ネタ記事でした。

追記 (2016/12/8)

2^{14}の方で必ず解が一通りになるという勘違いをしていてその間違いを残したまま公開していました。現在は修正済みです。

Knuthの補題Pと(Z/mZ)*の構造定理の補題

“The Art of Computer Programming”において線形合同法に関する定理を証明するのに使われる補題Pの証明を記す．また整数論の教科書で(Z/mZ)*の構造を決定するのに用いられる補題 (補題2と補題4)が補題Pから導かれることを見る．




線形合同法であるseedが0からいくつ進めたものかを得る その3 (法が2のべきとは限らない一般の場合)

法が2のべきとは限らない一般の場合で求める方法ができた。

法が素数pのべきの場合、{x_n}が法p^eで最大周期のLCG数列なら{x_{pn} / p}が法p^{e-1}で最大周期のLCGの数列になることを使う。
一般の法mに対しては、mを素因数分解したのち、各素数について計算し、中国剰余定理で結果を得る。

require "prime"

# プログラムの結果が正しいかテストする
def test()
  m = 2 * (3 ** 4) # テストする法mの値

  factors = Prime.prime_division(m)
  product = factors.map{|p,e| p }.inject(1, :*)
  if m % 4 == 0 then product *= 2 end

  # 最大周期になる全てのaとbを試す
  (0..m/product-1).each do |k|
    a = product * k + 1
    m.times do |b|
      next if gcd(b, m) != 1
      list = gen_list(a, b, m)
      list.each_with_index do |s, i|
        index = calc_index(s, a, b, m)
        if i != index
          raise "wrong result (#{s},#{a},#{b},#{m}) expected: #{i}, got: #{index}"
        end
      end
      puts "ok a=#{a}, b=#{b}"
    end
  end
end

# x_0 = 0, x_{n+1} = (ax_n + b) mod m で定義される数列の最初のm項を得る
def gen_list(a, b, m)
  list = []
  s = 0
  m.times do
    list << s
    s = (a * s + b) % m
  end
  list
end

# x_0 = 0, x_{n+1} = (ax_n + b) mod m で定義される数列でx_n = sとなるnを求める
# ただし数列の周期はmでなければならない
def calc_index(s, a, b, m)
  factors = Prime.prime_division(m)
  res = 0
  product = 1
  factors.each do |p, e|
    i = calc_index_pe(s, a, b, p, e)
    res = chinese(res, product, i, p**e)
    product *= p**e
  end
  res
end

# calc_indexでm = p^e (pは素数)の場合
def calc_index_pe(s, a, b, p, e)
  if e == 0 then
    0
  else
    l = (inverse(b, p) * s) % p
    pe = p ** e
    ss = step(a, b, s, p - l, pe)
    (aa, bb) = step_param(a, b, p, pe)
    (calc_index_pe(ss / p, aa, bb / p, p, e - 1) * p - p + l) % pe
  end
end

# f(x) = (ax + b) mod mに対してf^n(s)を求める
def step(a, b, s, n, m)
  if n == 0 then
    s
  elsif n % 2 == 0 then
    step((a * a) % m, ((a + 1) * b) % m, s, n / 2, m)
  else
    (step((a * a) % m, ((a + 1) * b) % m, s, n / 2, m) * a + b) % m
  end
end

# 関数 f(x) = (ax + b) mod m のn個の合成f^nは
# f^n(x) = (a'x + b') mod mと書ける。そのa', b'を得る
def step_param(a, b, n, m)
  if n == 0
    [1, 0]
  elsif n % 2 == 0
    step_param((a * a) % m, (a + 1) * b % m, n / 2, m)
  else
    (aa, bb) = step_param((a * a) % m, (a + 1) * b % m, n / 2, m)
    [(aa * a) % m, (aa * b + bb) % m]
  end
end

# 中国剰余定理
# x ≡ a (mod p), x ≡ b (mod q) となるxを求める (p, qは互いに素)
def chinese(a, p, b, q)
  (a * q * inverse(q, p) + b * p * inverse(p, q)) % (p*q)
end

# ax + by = k の解を得る (ただしkはaとbの最大公約数)
def extgcd(a, b)
  if b == 0 then
    return [1, 0]
  end
  (q, r) = a.divmod(b)
  (x, y) = extgcd(b, r)

  # 今 bx + ry = k が成り立っている。r = a - bqを使えば
  # ay + b(x - qy) = kを得る
  [y, x - q * y]
end

def gcd(a, b)
  if b == 0 then
    a
  else
    gcd(b, a % b)
  end
end

# a x ≡ 1 (mod m) となるxを得る (aとmは互いに素)
def inverse(a, m)
  (x, y) = extgcd(a, m)
  x
end

線形合同法数列が最大周期になる条件(ただし法が2のべきの場合) その2

線形合同法であるseedが0からいくつ進めたものかを得る その2 - oupoの日記の記事のアイディアを使ったら簡単な証明を得たので記す。
先の最大周期になる条件の証明はくどすぎた感があるので今回は簡潔に済ませる。

定義

　　x_0 = s
　　x_{n+1} = (a x_n + b) mod 2^k

で定義される数列{x_n}を(s,a,b,2^k)で定義されるLCG数列と呼ぶことにする。

この数列は、x_n = sとなる最小の整数n > 0がn = 2^kのとき最大周期であると呼ぶ。

また、a, b, 2^kが定まっている文脈において関数fを
　　f(x) = (ax + b) mod 2^k
の意味で使う。

最大周期のLCG数列の性質

  • 0 <= i, j < 2^kかつi ≠ jならばx_i ≠ x_jである
  • 数列の項x_0, x_1, …, x_{2^k-1}には2^k個の整数0, 1, …, 2^k-1が一回ずつ現れる

証明すべき命題

(s,a,b,2^k)で定義されるLCG数列が最大周期となる必要十分条件は以下の3つを全て満たすことである

  • aが奇数
  • bが奇数
  • k > 1ならばa ≡ 1 (mod 4)

ただしk > 0とする。

s = 0としても一般性を失わないこと

上記の命題を証明するにあたって、s = 0の場合のみ証明すれば十分である。
このことは、2つの初項s, s'に対して(s,a,b,2^k)で定義されるLCG数列が最大周期ならば(s',a,b,2^k)で定義されるLCG数列も最大周期であることを言えばよい。
実際、(s,a,b,2^k)で定義されるLCG数列{x_n}が最大周期ならx_i = s'となるiが存在する。そこでy_n = x_{i+n}と定義すれば{y_n}は(s',a,b,2^k)で定義されるLCG数列に一致するが、これは最大周期となる。

以下s=0とする。

bが奇数であることの必要性

bが偶数なら{x_n}の各項はすべて偶数になり最大周期とならない。よってbが奇数であることは最大周期であるために必要。

aが奇数であることの必要性

aが奇数とするとf(0) = f(2^{k-1})となってfは単射とならないので、{x_n}は最大周期にならない。

「k > 1ならば、a ≡ 1 (mod 4)」の必要性

(0,a,b,2^k)で定義されるLCG数列{x_n}が最大周期と仮定する。
数列の漸化式より
　　x_{n+2} = (a^2 x_n + (a + 1)b) mod 2^k
が得られる。
ここでnが偶数なら上の両辺は偶数である(∵ x_0 = 0とa + 1が偶数であることから帰納法)。
よって両辺を2で割って
　　x_{n+2} / 2 = (a^2 (x_n / 2) + (a + 1)b / 2) mod 2^{k-1}
を得るので数列{x_{2n} / 2}は(0, a^2, (a + 1)b / 2), 2^{k-1})で定義されるLCG数列だと分かる。
法を2^{k-1}とするこのLCG数列は仮定により最大周期であるので(a + 1)b / 2は奇数でなければならない。(この部分で仮定k > 1を使っている。k = 1なら2^{k-1} = 1を法とするLCGは"加える数"が奇数である必要はない)

しかしa ≡ 3 (mod 4)なら(a + 1)b / 2は偶数となるのでa ≡ 1 (mod 4)でなければならない

十分性の証明

(s,a,b,2^k)で定義されるLCG数列{x_n}が

  • aが奇数
  • bが奇数
  • k > 1ならば、a ≡ 1 (mod 4)

を満たすならば、{x_n}が最大周期なることを証明する。

kに関する帰納法による。
k = 1なら明らかである。
k > 1としてk - 1のとき成り立つとする。
すると{x_{2n} / 2}は(0, a^2, (a + 1)b / 2, 2^{k-1})で定義されるLCG数列なので帰納法の仮定によりこれは最大周期である。
よってx_{2i} / 2 = 0となる最小のi > 0は2^{k-1}なので、偶数iでx_i = 0を満たす最小のi > 0は2^kだと分かる。また奇数iに対してはx_iは奇数となるのでx_i = 0とならない。
したがって{x_n}が最大周期であることが証明された。

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