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

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

ポケモン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関数を見てください

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

NPC消費によるズレを防ぐ

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

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

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

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

実現できていないこと

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

メルセンヌツイスタの調律を行列で書く

擬似乱数生成器メルセンヌツイスタ(MT)には以下のような調律(tempering)と呼ばれる関数が登場します。

u32 tempering(u32 x) {
    x ^= (x >> 11);
    x ^= (x << 7) & 0x9d2c5680;
    x ^= (x << 15) & 0xefc60000;
    x ^= (x >> 18);
    return x;
}

MTの作者自身の講義ノートでtemperingはある行列をかけていることに相当することが書いてあります。

(引用した講義ノートでは横ベクトルに行列を右から掛けていますが、今回の記事では縦ベクトルに左から行列を掛けることにします)


実際にこの行列を求めてみました。行列Tは:

さらにこの行列の逆行列T^{-1}は

見やすさのために0はドットにしておきました。

説明図:

32-bitの整数は集合{0,1}の要素を成分とする32次元の(列)ベクトルと同一視できます。
{0,1}という集合に、通常の算術の和と積を2で割ったあまりとして和と積を定義しましょう。すると有理数全体や実数全体と同じような四則演算の法則が成り立ちます。つまり{0,1}はこの演算で体になるということです。この体をF_2と呼びます。
このとき整数のxorはベクトルの和に対応します。

そうすればx -> x xor (x >> 11)という関数がこんな行列に対応することが分かるはずです。

temperingのほかの行も同様に行列で書けます。それらの積をとれば最初の行列Tが得られます。

F_2が体となることが効いて、有理数を成分とする行列や実数を成分とする行列と同じように線形代数の理論が使えて逆行列有理数などのときと同じ方法で求めることができます。
(ただ今回の場合はまず有理数逆行列を求めてからそれをF_2の元に変換しました)

追記 (2014/10/23)

参考までにtemperingの逆関数は次のようになります。どうぞご利用ください。

u32 tempering_inv(u32 x) {
    x ^= (x >> 18);
    x ^= ((x << 15) & 0xefc60000);
    x ^= ((x << 7) & 0x9d2c5680) ^ ((x << 14) & 0x94284000) ^ ((x << 21) & 0x14200000) ^ ((x << 28) & 0x10000000);
    x ^= (x >> 11) ^ (x >> 22);
    return x;
}

追記 (2014/10/25)

プラスルツールのぷらすさんがtemperingの逆関数についてとっても詳しく考察してくれました。皆さん見ましょう。

ジョインアベニューの日替わりseedの周期

ポケモンBW2ジョインアベニューの日替わりseedで使われている漸化式がちょっとおもしろいです。

s[n+1] = ((s[n] * 0x5d583d6d6c078979 + 0x26a693) mod (2^64)) >> 32

まず、0がこの漸化式の不動点になることがわかりますが、実際のゲーム中では0になったときは全く新しい値で初期化されるようになっています。

この数列は大部分の初項については途中から循環が始まる準周期的数列になります。図で書くとσの形になるだとかカタツムリの形になるとかいってもいいかもしれません。

そこで循環部分がいくつ存在して、それらの周期はいくつなのか、ひたすら計算させて調べてみました。

周期 循環部分から1つ代表させたseed
83458 79f18008
123305 2134bc83
31719 86e648d5
2598 12719a8c
2466 6cf5a895
773 0e50bafa
6 6a3db7e0
25 04cbaf9f
15 1c905da4
1 00000000
1 eae02079

このように循環部分は11個あり、そのすべてが比較的短いことがわかりました。

さらに0x00000000に至るseed全体をグラフにしてみました。

面白いですね。

周期が6の循環部分に至るseed全体のグラフも書きたかったのですが、6a3db7e0で初めて循環部分にいたるようなseedが84256個もあってうまくいかず。他の5つのseedは13個とか8個とか0個とかとっても少ないのに…。

近況

乱数講座

乱数講座というものに参加しました。 (9月7日)

LCG Ring

LCG Ringというプログラムを書きました。 (9月24日)

ポケモンの乱数調整ユーザーなら見ただけで何を意味しているかは分かるはず。

ジョインアベニューのくじの解析

ポケモンBW2ジョインアベニューのくじの解析をしました。 (9月29日〜10月10日)

以上です。

LCGにおけるseedの検索を離散幾何に帰着させる

LCG(線形合同法)でseedを検索することを離散幾何の問題に帰着させるということを考えてみました。

まず線形合同法の漸化式で使われる関数を(Ax + B) mod Mとしておきます。
seedを1増やしたとき、その一つ先のseedというのはA増えるか、A-M増えるかのどちらかです。
seedをsとして0からsまで1ずつ増やしていった間に一つ先のseedがA-M増えた回数をxとすると、1つ先のseedが0以上MAX以下であるというのは

 0 \leq B + As - Mx \leq \rm{MAX}

でかけます。

図を示すとこんな感じ。MAX = M/2としていて、黄色い領域が1つ先のseedが0以上MAX以下である範囲です。

これを使うと0以上MAX以下が2連続するseedというのは次の不等式を満たす格子点(s,x)に対するsの値となります。


 0 \leq s \leq \rm{MAX}
 0 \leq x
 0 \leq B + As - Mx \leq \rm{MAX}

このようにしてLCGのseedを探す問題を離散幾何の問題に帰着できました。
変数を3つに増やせばある範囲の値が3連続するseedというものも離散幾何の問題に書き直せます。4連続以上でも同様です。

では、既存の離散幾何のプログラムに、この帰着させた問題を解かせてみましょう。
線型不等式を満たす格子点の個数を求める「LattE」というプログラムがあります。

これを試してみたところ、M = 2^64, MAX = 2^48-1で2連続のseedの個数は一瞬 (0.01秒)で求まりました。
M, MAXは同じまま3連続にすると4.2秒。4連続にすると、いつ計算が終わるのかわからない感じだったので途中で中断。

これは3連続のときの画像。

もう一つ残念なことに、LattEは格子点の個数を求めることはできるものの、格子点を全て列挙することができない模様。
線型関数でコストを与えて、コストが最小な点を1個出力することはできるのですが、3連続の場合にそれをしてみても少し時間が経った後途中で落ちてしまいました。

LattEに与える入力ファイルを生成するプログラムはこちら

補足

不等式を満たす格子点(s,x)に対するsの値は解の候補であって、必ず解になるとは限らないのじゃないかと思われれるかもしれません。ところが実は必ず解になるのです。

それは(s,x)平面のs軸に垂直な直線lに対して2つの互いに平行な直線


 B + As - Mx = 0
 B + As - Mx = M

のそれぞれの交点の間の距離が1であること、したがってどんな定数s_0に対しても

 s = s_0
 0 \leq B + As - Mx \lt M

を満たす格子点(s,x)は1つしかないことからわかります。

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