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

2016-09-14

線形合同法の周期が必ず法の値から定まる数の約数になること

2015-02-01

線形合同法である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

2015-01-28

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

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

線形合同法である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}が最大周期であることが証明された。

2015-01-28

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

周期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;
}

2014-11-05

バトルファクトリー必勝(?)乱数調整

ポケモンHGSSにおけるバトルファクトリーの必勝(?)乱数調整の動画を撮影しました。

作った(or 使った)ツールは以下

必勝seedを探すのとNPCのシミュレータを作るのどっちが手間かかっただろう…

オープンレベル5周目、ボーナスポケモンなし、seedは0x2263499bです。

この乱数調整はさきさんによるバトルステージ乱数調整とバトルトレイン乱数調整に影響を受けてやってみました。

必勝seedの探索

必勝seedは2^32通りしらみつぶしです。
まず次の2匹が最初の6匹に出るseedだけに絞ります。

ムクホーク技: ブレイブバード, ギガインパクト, インファイト, とんぼがえりアイテム:こだわりハチマキ
ムウマージ技: シャドーボール, サイコキネシス, １０まんボルト, みちづれアイテム: きあいのタスキ

大変な火力を持つ鉢巻ムクホークとタスキ道連れムウマージ。この2匹の力でなんとか必勝できるものを探そうということですね。ちなみに鉢巻ムクホークはさきさんのバトルトレイン乱数でも活躍しています。

この2匹を選び、あとの1匹を4匹の中からどれを取るか、全て試します。
必勝判定関数を用意しておき、7回の戦闘に全て勝てるか判定します。なお、二つの戦闘の間の交換は一切しないものとします。

この必勝判定関数は必勝かどうかを完全に判定するものではありません。すなわち判定関数が必勝でないといっても本当は必勝な可能性はあります。ですが、必勝判定関数が必勝だといえばそれは本当に必勝です。そういうふうに作ってあります。
ただし、みちづれのPP切れと相手がポケモンを交換する可能性を無視しているので厳密にそうとは言い切れません…そのためタイトルの必勝のあとに(?)をつけています

必勝判定関数は再帰によって実装されています。

手持ちのHPがすべて0なら「負け」と判定
相手のHPがすべて0なら「勝ち」と判定
場に出てる自分のポケモンのHPが0→交換。控えが2匹いる場合は両方試してどちらかで少なくとも一方で「勝ち」なら「勝ち」と判定
場に出てる相手のポケモンのHPが0→相手が次に出すポケモン最大2通りを両方試してどちらでも「勝ち」なら「勝ち」と判定
技を選択するか、交換をするか、両方試してどちらか少なくとも一方が「勝ち」なら勝ちと判定し、どちらもダメなら「負け」と判定
技を選択する場合
- 素早さが負けていたなら自分のHPを0にしたのち必勝判定
- 「みちづれ」を持っていて、相手が「どくどく」などの技を持っていない場合、自分と相手のHPを0にして必勝判定。「勝ち」になれば「勝ち」を返す。「負け」なら続けて他の選択肢を試す
- 一番相手にダメージを与えられる技を選ぶ。ただしダメージは最小乱数が出るものとして計算する。相手が「ゆうばく」などを持っていたら「負け」とする。
- 相手のHPが残っていて積み技などを持っていた場合「負け」とする。
- 攻撃した結果、相手のHPが残っている場合、自分のHPを0にする。その後必勝するか判定
交換する場合
- 相手が積み技などを持っている場合「負け」とする
- 交換先はすべて試す。すなわちどれか一つ「勝ち」が出れば「勝ち」。
- 「きあいのタスキ」を持っている場合を除いて、交換先は即倒れるものとする。

(こうして書いてみると何か穴があるような気がしてきたぞ…)

くわしくは次のcomplete_winnnable関数を見てください

https://github.com/oupo/factory-search/blob/master/factory-search.html#L129

2^32通りしらみつぶしをJavaScriptで書いて、しかもそれでうまくいったのだから、すごいですよね。

NPC消費によるズレを防ぐ

さきさんによるバトルステージ乱数調整ではポケモン選択画面でペラップを鳴かせることでNPC消費によるズレが起こらないようにされていました。しかしバトルファクトリーではそんなことできません。

そこでNPCのシミュレータを作ってNPCの消費が入らないような一定の長さ以上を持つ区間を探してそこを狙いました。
NPCのシミュレータを作るためにNPCの動きをわざわざエミュレータを使って解析しました。

またバトルファクトリーの外にもNPCがいるので、建物に入った瞬間に消費がいくつであったかをNPCの動きから特定します。

図で示したのがNPC消費の隙間です。婦人が上を向いた直後に隙間ができるのでそれを目印にします。

実現できていないこと

鉢巻ムクホークとタスキ道連れムウマージなどという超強いポケモンがいるので必勝seedを見つけられましたが、そういうポケモンが出ない周(つまりLv50では1〜6周目、オープンレベルでは1〜3周目)では難しいと思います。
どなたか挑戦してみませんか？

oupoの日記

線形合同法の周期を一般に決定する方法

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