32bit LCGのseedを手計算で求めよう (半分ネタ)
で与えられる線形合同法疑似乱数において、その上位16bit
が観測できるという状況を考えます。
そのとき $X_0$ の値を求めるにはどうすればよいでしょうか。
$Y_0$ を $\alpha$ とおくとき $X_0$ は
とかけます。よって $2^{16}$ 通りの $x$ すべてをしらみつぶして、 $Y_1$ の値が一致するものを調べればよさそうです。$Y_1$ の値が一致するものは複数個ある可能性がある場合があるのでその場合は $Y_2$ の値を調べます。これで疑似乱数のseedを特定できるでしょう。
しかし、$2^{16}$ 通り調べるのは手計算ではできないですね。
そこで手計算できる方法を考えます。
まず、
であることに注意します。
特に係数がそれぞれ $2^{16}$ の倍数に1足した値と$2^{16}$ の倍数になっているので、それらを$2^{16} A + 1$, $2^{16} B$とおいてみます。
すると上の $X_0$ の候補の値について
が成り立ちます。ここで $Y_{2^{16}}$ を $\beta$ とおけば
なのでこれを解くことにより $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}個先を考えるのです。
ここに出てくる係数をあらためて $2^{16} A + 1 = 0xf7e90001, 2^{14} B = 0xf1cc4000$ とおきます。
すると $X_0$ の候補の値について
が成り立ちます。ここで $X_{2^{14}}$ を $2^{16} \beta + y$ とおけば
$x = 2^{14} x' + x'', y = 2^{14} y' + y'', 0 \ge x', y' < 4, 0 \ge x'', y'' < 2^{14}$ とおけば
上位18bitをとることにより
$y' - x' = \delta$とけば
(1)
この式に適する$\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)は
です。整理すると
$-203157 + \delta$は4の倍数にならざるを得ないことから$\delta = 1, -3$です。
$\delta = 1$のとき:
両辺を4で割り、
63465の2^{16}を法とした逆元はユークリッドの互除法により20569と求まるので
計算すると
と求まります。
このとき$x' = 1$なので$\delta = 1$と整合します ($\delta = y' - x'$をみたす$0 \ge y' < 4$があります)。
$\delta = -3$のとき:
同様の計算で
と求まります。
このとき$x' = 3$なので$\delta = 1$と整合しません ($\delta = y' - x'$をみたす$0 \ge y' < 4$がありません)。
よって解は30435の一つです。
ゆえにseedはです。
2^{14}回乱数を進めるなんて2^{16}通り機械でしらみつぶしするよりずっと大変じゃないかということで、半分ネタ記事でした。
追記 (2016/12/8)
2^{14}の方で必ず解が一通りになるという勘違いをしていてその間違いを残したまま公開していました。現在は修正済みです。
線形合同法の周期を一般に決定する方法
この記事ではKnuthの補題Pと(Z/mZ)*の構造定理の補題 - oupoの日記で与えた補題を利用する.
線形合同法である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数列の性質
証明すべき命題
(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が奇数であることは最大周期であるために必要。
「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}が最大周期であることが証明された。
線形合同法であるseedが0からいくつ進めたものかを得る その2
- 前回の記事: 線形合同法であるseedが0からいくつ進めたものかを得る - oupoの日記
- 続きの記事: 線形合同法であるseedが0からいくつ進めたものかを得る その3 (法が2のべきとは限らない一般の場合) - oupoの日記
周期2^k, 初項0のLCG数列{x_n}に対して、{x_{2n} / 2}って周期2^(k-1)のLCG数列だよね、というアイディアで再実装してみました。
#include <stdio.h> #include <stdint.h> typedef uint32_t u32; /* x[0] = 0, x[n+1] = (a * x[n] + b) mod 2^k で与えられる数列{x[n]}において、x[n] = sを満たすnを求める ただし数列の周期は2^kになっていなければならない */ u32 calc_index(u32 a, u32 b, u32 s, u32 k) { if (k == 0) { return 0; } else if ((s & 1) == 0) { return calc_index(a * a, (a + 1) * b / 2, s / 2, k - 1) * 2; } else { return calc_index(a * a, (a + 1) * b / 2, (a * s + b) / 2, k - 1) * 2 - 1; } } int main(void) { u32 A = 0x41c64e6d; u32 B = 0x6073; printf("%.8x\n", calc_index(A, B, 0, 32)); printf("%.8x\n", calc_index(A, B, B, 32)); printf("%.8x\n", calc_index(A, B, A * B + B, 32)); printf("%.8x\n", calc_index(A, B, A * A * B + A * B + B, 32)); return 0; }