メルセンヌツイスタの調律を行列で書く
擬似乱数生成器メルセンヌツイスタ(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)
ジョインアベニューの日替わり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個とかとっても少ないのに…。
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以下であるというのは
でかけます。
図を示すとこんな感じ。MAX = M/2としていて、黄色い領域が1つ先のseedが0以上MAX以下である範囲です。
これを使うと0以上MAX以下が2連続するseedというのは次の不等式を満たす格子点(s,x)に対するsの値となります。
このようにして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つの互いに平行な直線
のそれぞれの交点の間の距離が1であること、したがってどんな定数s_0に対しても
を満たす格子点(s,x)は1つしかないことからわかります。
■
重複なしにランダムに選んだ後ある特定のseed値になるようなseed値を求める その2 - oupoの日記に追記をしておいた。誰も興味がないようなネタだろうけど…。
確率計算
過去に書いたプログラムが、投稿したwebサービスから消えていたのでブログに残しておく。(2011年11月に作成)
ブラックシティとホワイトフォレストの最初に決まる13人の中にグレースがいる確率を求めるプログラム (だったと思う)。だだぢぢさんがTwitterで話題にしていて触発されて書いた。
各住人の重みのデータはktxadさんが公開されたものによる (URLは紛失)。
(追記2014/10/15: ktxadさんが公開されていたソースコードは現在消えていますが、dorara32さんがかわりに重みのデータを公開されています: http://note.chiebukuro.yahoo.co.jp/detail/n166836)
# encoding: utf-8 require "mathn" # 整数同士で割り算したら分数になるように # 重みのついたデータから重複しないように何件かランダムに選ぶとき、 # 特定の1件が選ばれる確率を計算する Candidate = Struct.new(:id, :weight) CANDIDATES = [ [ 0, 20], [ 1, 20], [ 2, 20], [ 3, 20], [ 4, 20], [ 5, 20], [ 6, 20], [ 7, 10], [ 8, 20], [ 9, 10], [10, 20], [11, 6], [12, 20], [13, 4], [14, 1], [15, 20], [16, 20], [17, 20], [18, 20], [19, 20], [20, 20], [21, 10], [22, 20], [23, 20], [24, 10], [25, 20], [26, 6], [27, 4], [28, 20], [29, 1] ].map{|(id,weight)| Candidate.new(id, weight) } TARGET_ID = 29 CALCLATORS = [ [:MonteCarlo, "モンテカルロ法による確率"], [:Naive, "ナイーブな方法による確率"], [:Fast, "高速化した方法による確率"], ] def main [1, 2, 3, 10, 13].each do |num_choice| puts "選ぶ個数: #{num_choice}" CALCLATORS.each do |(class_name, description)| klass = ProbabilityCalclator.const_get(class_name) if class_name == :Naive and num_choice >= 5 pr = "(時間がかかるため省略)" else pr = klass.new.calc(num_choice) end print " #{description}: #{pr}" if pr.kind_of?(Rational) print " ≒ #{pr.to_f}" end puts end end end def sumup_weight(candidates) candidates.map {|x| x.weight }.inject(:+) end ProbabilityCalclator = Module.new # モンテカルロ法 class ProbabilityCalclator::MonteCarlo TRY_TIMES = 10000 def calc(num_choice) succeeded_count = 0 TRY_TIMES.times do succeeded_count += 1 if try(num_choice) end succeeded_count.to_f / TRY_TIMES end def try(num_choice) candidates = CANDIDATES.dup num_choice.times do index = random_choose(candidates) candidate = candidates.delete_at(index) if candidate.id == TARGET_ID return true end end false end def random_choose(candidates) random_number = rand(sumup_weight(candidates)) sum = 0 candidates.each_with_index do |candidate, index| sum += candidate.weight if sum > random_number return index end end raise "unreachable" end end # ナイーブな方法 # CANDIDATES.size = 30, num_choice = 13 の場合 # 多分 29P13 回ぐらい探索しないといけなくて遅い class ProbabilityCalclator::Naive def calc(num_choice) @visited = [false] * CANDIDATES.size @total_sum = sumup_weight(CANDIDATES) recur(num_choice) end # 残りnum個での確率を計算 def recur(num) return 0 if num == 0 result = 0 CANDIDATES.each_with_index do |candidate, index| next if @visited[index] # 今すでに選択している候補はスキップ if candidate.id == TARGET_ID # TARGET_IDだった場合、確率は1。残りの分の探索はやらなくてOK pr = 1 else # candidateを選び、まだ決定していない残りnum - 1個分での確率を計算する @visited[index] = true @total_sum -= candidate.weight pr = recur(num - 1) @total_sum += candidate.weight @visited[index] = false end result += (candidate.weight / @total_sum) * pr end result end end # 同じ重みの値をまとめて扱い、さらにメモ化することで高速化 class ProbabilityCalclator::Fast def calc(num_choice) # 重みの値をキーにしてその重みの候補があと何個残っているかを記憶する @weight_to_count = CANDIDATES.each_with_object(Hash.new(0)) {|candidate, r| r[candidate.weight] += 1 } @total_sum = sumup_weight(CANDIDATES) @memo = Hash.new @weights = @weight_to_count.keys @target_weight = CANDIDATES.find{|x| x.id == TARGET_ID }.weight recur(num_choice) end def recur(times) # メモ化 # (@weight_to_countが同じ内容ならばrecur_internalの結果も同じだから # 記憶しておいて使いまわす) @memo[@weight_to_count.dup] ||= recur_internal(times) end def recur_internal(times) return 0 if times == 0 result = 0 @weights.each do |weight| count = @weight_to_count[weight] # 重みがweightの値であり未選択の候補の数 if weight == @target_weight pr = 1 / count # TARGET_ID を選んだときの確率 if count > 1 # TARGET_IDとweightは同じだがTARGET_IDではないものを選んだときの確率 pr += ((count - 1) / count) * choose(times, weight) end else pr = choose(times, weight) end result += pr * ((count * weight) / @total_sum) end result end # weightを選択した場合の残りtimes - 1個分での確率を計算 def choose(times, weight) if @weight_to_count[weight] > 0 @weight_to_count[weight] -= 1 @total_sum -= weight pr = recur(times - 1) @total_sum += weight @weight_to_count[weight] += 1 pr else 0 end end end if $0 == __FILE__ main() end
実行結果は次の通り。
なお「モンテカルロ法による確率」は不正確な結果で、「ナイーブな方法による確率」と「高速化した方法による確率」は完全に正確な結果である。
選ぶ個数: 1 モンテカルロ法による確率: 0.0025 ナイーブな方法による確率: 1/462 ≒ 0.0021645021645021645 高速化した方法による確率: 1/462 ≒ 0.0021645021645021645 選ぶ個数: 2 モンテカルロ法による確率: 0.0035 ナイーブな方法による確率: 9729246465/2204009196532 ≒ 0.004414340230661892 高速化した方法による確率: 9729246465/2204009196532 ≒ 0.004414340230661892 選ぶ個数: 3 モンテカルロ法による確率: 0.0064 ナイーブな方法による確率: 1568664824856411452425433501221/232169860937309598063431523052200 ≒ 0.006756539451431991 高速化した方法による確率: 1568664824856411452425433501221/232169860937309598063431523052200 ≒ 0.006756539451431991 選ぶ個数: 10 モンテカルロ法による確率: 0.0278 ナイーブな方法による確率: (時間がかかるため省略) 高速化した方法による確率: 65770250293419026078129194848669982508897034077533201092414085189931935552537702182324234499375672068027368327207661273376859495758981364818701358782794795202517/2475647499344346203202200274227924547084627755009408713100317804864210058755132452548664905544294615343129859802504836359196285761348970446409580976203759822039600 ≒ 0.026566888182117086 選ぶ個数: 13 モンテカルロ法による確率: 0.0396 ナイーブな方法による確率: (時間がかかるため省略) 高速化した方法による確率: 72920965323202474895401221315730728325203050221435766480570752583389314598872312680973006197795827752084921849225068647216327771173023300970901986908137886930833925903506944736597859681/1932080335186385797665285107673085006540033135470821753828553162094595243823360809973187589984927249351524831714776478574108021855863498572634867939515180025821224902545585179805950634000 ≒ 0.03774220149917724
C++による実装も実は作っていたのでした。
#include <iostream> #include <vector> #include <map> // 重みのついたデータから重複しないように何件かランダムに選ぶとき、 // 特定の1件が選ばれる確率を計算する struct Candidate { int id; int weight; }; const Candidate CANDIDATES[] = { { 0, 20}, { 1, 20}, { 2, 20}, { 3, 20}, { 4, 20}, { 5, 20}, { 6, 20}, { 7, 10}, { 8, 20}, { 9, 10}, {10, 20}, {11, 6}, {12, 20}, {13, 4}, {14, 1}, {15, 20}, {16, 20}, {17, 20}, {18, 20}, {19, 20}, {20, 20}, {21, 10}, {22, 20}, {23, 20}, {24, 10}, {25, 20}, {26, 6}, {27, 4}, {28, 20}, {29, 1} }; const int NUM_CANDIDATES = sizeof(CANDIDATES) / sizeof(Candidate); const int TARGET_ID = 29; double calc(int num_choice); void initialize(); double recur(int times); double recur_internal(int times); double choose(int times, int weight); int main(void) { for (int num_choice = 1; num_choice <= NUM_CANDIDATES; num_choice ++) { std::cout << "num_choice = " << num_choice << ": " << calc(num_choice) << std::endl; } } typedef std::map<int,int> WeightToCount; // 重みの値をキーにしてその重みの候補があと何個残っているかを記憶する WeightToCount g_weight_to_count; std::map<WeightToCount, double> g_memo; int g_total_sum; int g_target_weight; std::vector<int> g_weights; double calc(int num_choice) { initialize(); return recur(num_choice); } void initialize() { g_weight_to_count.clear(); g_weights.clear(); g_memo.clear(); g_total_sum = 0; for (int i = 0; i < (sizeof CANDIDATES / sizeof(Candidate)); i ++) { Candidate candidate = CANDIDATES[i]; int weight = candidate.weight; if (g_weight_to_count.find(weight) != g_weight_to_count.end()) { g_weight_to_count[weight] += 1; } else { g_weight_to_count[weight] = 1; g_weights.push_back(weight); } g_total_sum += weight; if (candidate.id == TARGET_ID) { g_target_weight = weight; } } } double recur(int times) { // メモ化 // (g_weight_to_countが同じ内容ならばrecur_internalの結果も同じだから // 記憶しておいて使いまわす) if (g_memo.find(g_weight_to_count) != g_memo.end()) { return g_memo[g_weight_to_count]; } else { return g_memo[g_weight_to_count] = recur_internal(times); } } double recur_internal(int times) { if (times == 0) return 0; double result = 0; for (int i = 0; i < g_weights.size(); i ++) { int weight = g_weights[i]; // 重みがweightの値であり未選択の候補の数 int count = g_weight_to_count[weight]; double pr; if (weight == g_target_weight) { pr = 1.0 / count; // TARGET_ID を選んだときの確率 if (count > 1) { // TARGET_IDとweightは同じだがTARGET_IDではないものを選んだときの確率 pr += ((double)(count - 1) / count) * choose(times, weight); } } else { pr = choose(times, weight); } result += pr * (count * weight) / g_total_sum; } return result; } // weightを選択した場合の残りtimes - 1個分での確率を計算 double choose(int times, int weight) { if (g_weight_to_count[weight] > 0) { g_weight_to_count[weight] -= 1; g_total_sum -= weight; double pr = recur(times - 1); g_total_sum += weight; g_weight_to_count[weight] += 1; return pr; } else { return 0; } }
asm.jsを試す
Firefoxで実装されているというasm.jsを試しました。これはJavaScriptのプログラムのうち、一部を非常に制限されたルールにしたがって書いてそこに"use asm"というディレクティブをつければ、その部分が非常に高速化されるというしろものです。
以前、「メルセンヌツイスタにおいて与えるseedを0から1000000まで変化させ、それぞれのseedで一番目の乱数を得る」というプログラムでいろんな言語のベンチマークをとりました。
これは要するに32bit整数の演算ばっかりをたくさんやるけれど、配列参照などはまったくしないプログラムです。
このベンチマークをasm.jsに対してもやってみました。
- http://oupo.github.io/tmp/mt-benchmark-use-asm.html
- http://oupo.github.io/tmp/mt-benchmark-no-use-asm.html
結果:
- use asm: 1773 msec
- no use asm: 987 msec
あれ？あれれれれ。"use asm"を無効にした方が高速になっています…。
しかも…
>cl.exe -Ox mt-benchmark.c
Microsoft(R) C/C++ Optimizing Compiler Version 18.00.21005.1 for x86
Copyright (C) Microsoft Corporation. All rights reserved.
mt-benchmark.c
Microsoft (R) Incremental Linker Version 12.00.21005.1
Copyright (C) Microsoft Corporation. All rights reserved.
/out:mt-benchmark.exe
mt-benchmark.obj
>type time.rb
t = Time.now
system *ARGV
puts "#{Time.now - t} sec."
>ruby time.rb .\mt-benchmark.exe
232881862
0.968796 secということなので、"use asm"なしはCと遜色のない速度で動いています。(C版の時間はプロセス起動と終了まで含まれてしまっているとはいえ)
これは二重に驚きの結果です。
環境: Firefox 27.0.1
