線形合同法であるseedが0からいくつ進めたものかを得る

sid-searchを作り直してみた - oupoの日記でも使ったアルゴリズムを解説する。

問題

線形合同法の漸化式に使われる関数(x * A + B) mod 2^Nをf(x)とおき、この線形合同法の数列は最大周期であるとする。
を初項を0としたこの線形合同法の数列とする。すなわちa_0 = 0, a_{n+1} = f(a_n)
与えられたseed、xに対してf^i(0) = xとなるi ∈ 0...2^Nを求めたい。

思考

a_nは偶数、奇数を繰り返すことからxの偶奇からiの偶奇がわかる。iの第0ビットがわかったので、同じことを繰り返して第1ビット、第2ビット...と次々に求めていけないだろうかと考える。

  • iが偶数の場合、0にf^2を何回適用すればxになるかを求めて、それを2倍すれば0にfを何回適用すればxになるかが得られる。
  • iが奇数の場合、0にf^2を何回適用すればf(x)になるかを求めて、それを2倍して1引けば0にfを何回適用すればxになるかが得られる。

この考えを使って繰り返せばよい。

ここで、f^i(0) = xとなるiの偶奇はxの偶奇から得られたが、fの代わりにf^2を当てはめた場合のiの偶奇は得られるか？

xの第1ビットから分かるのではないかと予想できる。

実際:

n a_n a_n mod 4
0 0x00000000 0b00
2 0x6c078966 0b10
4 0xdbffe6dc 0b00
6 0x992415e2 0b10
8 0x895277f8 0b00
10 0x014d329e 0b10
12 0x75ab8f54 0b00
14 0xe0d45b9a 0b10
16 0xe69948f0 0b00
18 0x964b4cd6 0b10

一般に、g=f^{2^k}とおいたとき、g^i(0) = xとなるiを2で割った余りはxの第kビットに一致する。

これはの下位kビットだけ注目したとき周期が2^kであることから証明できる。

k=1を例にとって解説してみる。

の周期は4で、この表を考える。

n a_n mod 2 a_n mod 4
0 0b0 0b00
1 0b1
2 0b0
3 0b1

a_n mod 2の値からa_n mod 4の第0ビットは分かる。

n a_n mod 2 a_n mod 4
0 0b0 0b00
1 0b1 0b?1
2 0b0 0b?0
3 0b1 0b?1

a_2 mod 4が0b00だとすると周期が4なことに反してしまうのでa_2 mod 4 = 0b10である。

n a_n mod 2 a_n mod 4
0 0b0 0b00
1 0b1 0b?1
2 0b0 0b10
3 0b1 0b?1

よって0を初項、f^2を漸化式とする数列はmod 4で見たとき0b00, 0b10を繰り返す。

プログラム

# f^m(0) = sとなるm ∈ 0...2^Nを求める
def seed_to_index(s)
  seed_to_index0(s, 0)
end

# g = f^{2^n}としてg^m(0) = sとなるm ∈ 0...2^{N-n}を求める
def seed_to_index0(s, n)
  if n == N
    0
  elsif ((s >> n) & 1) != 0
    (seed_to_index0(step_seed(s, 2**n), n + 1) * 2 - 1) % 2**(N-n)
  else
    seed_to_index0(s, n + 1) * 2
  end
end

前回の記事では下から5行目がseed_to_index0(step_seed(s, -2**n), n + 1) * 2 + 1となっていたのだが、-2^n進める定数をあらかじめ表にしておくより、2^n進める定数を表にしておく方が自然なので今回の方が好みだ。

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