Hatena::ブログ(Diary)

菊やんの雑記帳

2012-04-15 GCJ2012予選

[] GCJ予選 19:54  GCJ予選 - 菊やんの雑記帳 を含むブックマーク

証明できそうな問題がないかなと探したところA問題だけだった。BとCもできるかもしれないけどめんどくさいからいいや

2010-12-04 FFT

[] FFT書かないといけないんだけど 20:20  FFT書かないといけないんだけど - 菊やんの雑記帳 を含むブックマーク

世の中にはCooley-Tukeyの解説しかないんですよ。並列の時代なのでStockhamとか知りたいんですが…

で、難しい英語の解説とか読まなくてもサンプルコードがあれば実装できるので別に問題ないんだけど、やっぱりちゃんと理解したい。

FFTの解説ってΣをごにょごにょする系とクロネッカー積をごにょごにょする系しかないんですよ。

Σだらけの式とか読みたくないし、世の中全部行列で説明できるなんてものきもい。高階のテンソルなんかももちろん行列ですよ。行列の行列とか行列の分割とかそういった世界ですよわけがわかりません。

最初はもっと数学数学したテンソル代数を使うものではないかとも考えたのですが、それは間違いで、抽象的に考えると計算が効率的になったのかなんてさっぱり分からないし、テンソルのメモリ上の表現なんかきれいさっぱりなくなってこれはもうなにをやってるかわからない。

でも、行列だけにしぼるとメモリ上の配置とか計算量だとか見やすいかというとまったくそうでない。

テンソル的なものを表現するのはアインシュタイン先生がわりと記法を編み出してくれてるんだけど、あれはみんな4次元の話だからインデックスには何次元かの型情報がいるんですよ。あと普通の内積しか使わないからインデックスの上下なんて対して意味がないんです。

というわけで、テンソルはこう書くことにしよう

x[i:32][j:32]
y[k:32]
z

xは32x32の行列だ。yは32次元のベクトルだ。zはスカラーだ。インデックスは縮約するときに使うぞ。

Cのプログラムに変換する方法はまったく自明だ

complex x[32][32];
complex y[32];
complex z;

あとはテンソル積と縮約だ。ただ掛け算するのをテンソル積だとする。

w[k:32][i:32][j:32] = y[k:32]*x[i:32][j:32]

テンソル積の結果はただたんにインデックスを並べるだけだ。同じ名前のインデックスがあると縮約する

v[i:32] = y[j:32]*x[i:32][j:32]

これはベクトルと行列の積だ。この二つをプログラムで書くと

for (k) for (i) for (j) w[k][i][j] = 0;
for (k) for (i) for (j) w[k][i][j] += y[k]*x[i][j];
for (i) v[i] = 0;
for (i) for (j) v[i] += y[j]*x[i][j];

というわけで、積は機械的に総和の計算になる。

とか書いてると全部のインデックスに型をつけるのはめんどくさいので、一度型を書いたインデックスは後では書かなくてもいいことにしよう。あと先に型を宣言しておいてもいいことにしよう。まぎらわしくないときはコンマでつなげてもいいことにしよう。

w[k:32][i:32][j:32] = y[k]*x[i,j]

とか

i,j,k:32

v[i] = y[j]*x[i,j]

テンソルの要素にアクセスしたいときもあるので、それは丸括弧を使うことにしよう

omega[j:32][k:32]

omega(j,k) := exp(-2*pi*j*k*sqrt(-1)/32)

で定義するみたいに使う。

グループ化を行うテンソルが重要で、(N*M)要素ベクトルからNxM行列を作るような操作Gをテンソルで書く。

i:N, j:M, k:N*M
y[i,j] = G[i,j,k]*x[k]

のように使うとグループ分けができる。定義は

G(i,j,k) := if k = i*M+j then 1 else 0

グループ分けの解除も同じGでできる

x[k] = G[i,j,k]*y[i,j]

完全シャッフルテンソルPは一度グループ分けして逆方向に解除するとできる。トランプのシャッフルのこと。山を二つに分けて完全に交互に挿入するシャッフルFFTで重要。

i,j:2*N
y[i] = P[i,j]*x[j]
where
m:N, n:2
P[i][j] := G[m,n,i]*G[n,m,j]

当たり前だが、Gを同じ方向に2回すると元に戻る。恒等テンソルをidで書くことにすると

n,s:N, m,t:M, j,k:N*M
G[n,m,j]*G[n,m,k] = id[j,k]
G[n,m,j]*G[s,t,j] = id[n,s]*id[m,t]

こんな感じで0と1しかないテンソルは代入だと理解したほうがいいかも?

n:N, m:M, j:N*M
G(n,m)[j]*T[j] = T(n*M+m)
G[n,m](j)*U[n,m] = U(j div M, j mod M)

とかやってればP(i,j)とかも計算できるかも

i,j:2*N
P(i,j) = G[m:N][n:2](i)*G[n][m](j)
= G(i mod 2, i div 2, j)
= id((i mod 2)*N+(i div 2), j)
= if iが偶数 then id(i div 2, j) else id(i div 2 + N, j)
= if j < N then id(i, j*2) else id(i, (j-N)*2+1)

というわけで、P(i)[j]はiの偶奇によりjに(i div 2)か(i div 2 + N)を代入し、P[i](j)はjとNの大小で、iに((j-N)*2+1)か(j*2)を代入するテンソルだと分かりました。

離散フーリエ変換を書くとどうなるかというと、まずomega[i:N][j:N]を

omega(i,j) := exp(-2pi*i*j*sqrt(-1)/N)

で定義して、ベクトルx[i:N]のフーリエ変換y[i:N]は

y[i] = omega[i,j]*x[j]

と書ける。サイズが偶数のときを考えて、omegaをインデックスの偶奇と大小で2分割する。iを分割するかjを分割するかの流儀がある。

iを偶奇にjを大小に分割するほうが周波数間引きで、

jを偶奇にiを大小に分割するほうが時間間引きというのでたぶんあってる。

iを偶奇で分割するときは

omega'[i':N][i0:2][j0:2][j':N]=G[i', i0, i]*G[j0, j', j]*omega[i,j]

これを計算すると

omega'(i', 0, j0, j')
=G(i', 0)[i]*G(j0, j')[j]*omega[i,j]
=omega(i'*2, j0*N+j')
=omega(i', j')
=omega(i')[k:N]*id[k](j')

omega'(i', 1, j0, j')
=G(i', 1)[i]*G(j0, j')[j]*omega[i, j]
=omega(i'*2+1, j0*N+j')
=omega(i', j')*omega(1, j')*(-1)^j0
=omega(i')[k:N]*id[k](j')*omega(1, j')*(-1)^j0

というわけで、バタフライ演算B[i0:2][j0:2][i:N][j:N]を

B(0, j0, i, j) := id(i, j)
B(1, j0, i, j) := id(i, j)*omega(1, j)*(-1)^j0

で定めると

omega'[i':N][i0:2][j0:2][j':N] = omega[i':N][k:N]*B[i0:2][j0:2][k:N][j':N]

と書けました。両辺にGを2個掛けてomega[i:2*N][j:2*N]に戻すと

i0,j0:2, i',j':N

omega[i, j]
= G[i', i0, i]*G[j0, j', j]*omega'[i', i0, j0, j']
= G[i', i0, i] * omega[i', k] * B[i0, j0, k, j'] * G[j0, j', j]

というわけで、インデックスがjのベクトルフーリエ変換するには、jを大小で2分割して新しいインデックスを j0, j' としてバタフライ演算を行い、インデックスが i0, k になって、各 i0 に小さいフーリエ変換を行い、インデックスが i0, i' になり、最後に偶奇分割を解除すればできることがわかりました。めでたしめでたし。計算は機械的にできるようになったけどインデックスつけるのがめんどくせえ。

iに関する分割解除をする前に大小で分割解除して、さらに大小に分割するようにすると

s0:2, s':N, s:2*N
omega[i, j]
= G[s', s0, i] * G[s0, s', s] * G[i0, i', s] * omega[i', k] * B[i0, j0, k, j'] * G[j0, j', j]
= P[i,s] * G[i0, i', s] * omega[i', k] * B[i0, j0, k, j'] * G[j0, j', j]
= P[i,s] * F[s,j]

とかいった感じで進めてればCooley-Tukeyとか簡単に理解できるのではないだろうか、続くかも?

2010-06-21 ICFPC終了

[] ICFPC2010参戦記 21:26  ICFPC2010参戦記 - 菊やんの雑記帳 を含むブックマーク

問題設定とかは http://d.hatena.ne.jp/ku-ma-me/ を見ればいいとして、

金曜日21:00〜23:00

問題読みタイム。さっぱりわからないが、一番下に最初にやるべきことが書いてあったので、0を送ればいいっぽい、サーバが死に掛けてて全く送れない。やっと反応があったときにはすでに登録されてるよって出る。

それはいいとして、次は回路を理解せよということだったので、まず回路のフォーマットを推測して、適当に作った回路を送り始める。「22022022022022022」を返してくる回路を見つけてこれがヒントになるかと思ったが、時間切れで電車に乗る

土曜日

帰宅して、適当に回路を作って送りつける作業を続ける。結論として、次のような回路を考えればゲートの動作を推測できることが分かった。

  X +---+
  | |   |
+-*-*-+ |
|     | |
+-*-*-+ |
  | |   |
  | +---+
  A
  | +---+
  | |   |
+-*-*-+ |
|     | |
+-*-*-+ |
  | |   |
  | +---+
  B
  | +---+
  | |   |
+-*-*-+ |
|     | |
+-*-*-+ |
  | |   |
  | +---+
  C

最初に1ゲートのを送ってAの出力を調べて、次に2ゲートのを送ってBの出力を調べて、を繰り返せば入力と出力のペアはいくらでも得られる。これを左右反転した回路に対しても行って、

左入出力
02120112100002120
00100202221110100
01221000201200221
02110011002121022
00122120010101102
01211021122222000
02101101211111211
00112002101201012
右入出力
22022022022022022
21202210221202210
20221022020221022
22102202212020221
21121022111202211
20202202021120202
22120221202112022
21112020221102210

これをじっと眺めてれば、左に2を入力したときに左からは0,1,2の全通りが出てるので、

in  | out
2 ? | 0 X
2 ? | 1 ?
2 X | 2 ?

よく見ると、「2->0」の直後に「2->2」があるので、Xが等しいことが分かる。途中で計算間違いを莫大にしたので時間がかかったけど、記録によると1:30ごろにゲートは確定したようだ。

回路作成

というわけで、今度は任意の出力を出せる回路を作る方法を考えないといけないわけだが、最初に考えたのは常に0を出す基本回路を作るんじゃねとか思ってみたりしたんだけど、どう頭をひねっても想像もできないので、やめてもっと直接にやる方法を思いついた。

とりあえず、回路全体をこの形に固定することにする

 inp
  | +---+
  | |Y  |
+-*-*-+ |
|     | |
+-*-*-+ |
  | |X  |
  | +---+------X
  |     |         残りの回路
 out    +------Y'

この回路はinpとoutのtrit列が与えられたとき、 inp[0] = out[0] ならば解があって、

Y = 0::Y'
inp[0] = out[0]
X[0] = inp[0]×0-1
Y[k] = inp[k]-out[k]
X[k] = inp[k]×Y[k]-1

のように、XとY'のtrit列を完全に決定できる。しかも、Y'の列の長さはoutの長さより1短い。

よって、inpからoutを出す回路を作る問題はinp[0]=out[0]のときは、1短い入出力の問題に単純化できた。

同様に inp[0]-2=out[0]のときは、ちょっと改造して Y = 2::Y' にすればいいので、

    +---+ +-------+
    |   | |       |
    | +-*-*-+     |
    | |     |     |
    | +-*-*-+     |
    |   | |       |
    +---+ |       |
          |       |
        +-+       |
        |         |
        | +---+   |
        | |   |   |
      +-*-*-+ |   |
      |     | |   |
      +-*-*-+ |   |
        | |   |   |
        | +---+---+
 inp    |     |
  | +---+     |
  | |Y        |
+-*-*-+       Y'
|     |
+-*-*-+
  | |X
  | +-----X
  |
 out

とつなぐと、Y = 2::f(Y') となる。このfが全単射になるように繋がないといけないので右出力を使えないのです。

これでkey prefixは通った。朝の7時半ごろ。そのまま寝て起きて、夜になったら車フォーマットを研究することにする。

日曜日

午前2時ぐらいにやっと車のフォーマットがほぼわかった。引き続き燃料の解析。午前3時半に最初の燃料が通過する。この時点で簡単な車にがんがん給油してればもっと点数も高かったんだけど、メインの問題を考え始めてしまった。一人でやってるとしかたがないよね。

このへんで、サンプルの車を解こうとするが、さっぱり理解できない。こんな車解けるのか?ランダムに行列を生成したり、ちょっとした変異を加えても全く答えに近づいているような気配がない。9時半まで悩んで全く解決しなかったので寝る。

21時ごろに起きて、大衆車に給油する仕事をすませる。そのあとサンプル車の解き方を考えるのだが、こういう問題の基礎知識がないのでさっぱりわからない。非線形の最適化問題とかやったことないよ。

23時にあきらめてWORKINGを見たりする。さすがに時間切れ。

2010-04-08 ferNANDoでQuine(3)

[] なんも説明もしてなかった気がするので 18:37  なんも説明もしてなかった気がするので - 菊やんの雑記帳 を含むブックマーク

とりあえず少し書いておくと、まず第1回で書いたようにROMは作れた。

A0 0 0
A1 0 0
A2 0 0
...
Ae 0 0
Af 1 1
## こんな感じにADDR <- 0x7FFF みたいに、アドレスをセットして
## ROMのコードに突っ込むと

## BEGIN ROM
if (ADDR-- == 0) {
  DATA <- 'H';
}
if (ADDR-- == 0) {
  DATA <- 'e';
}
if (ADDR-- == 0) {
  DATA <- 'l';
}
...
if (ADDR-- == 0) {
  DATA <- '\n';
}
## END ROM

## ここで DATA に対応するバイトが読まれてる

次にCPUを作る。

PC <- 0x0000 # プログラムカウンタを0に
# リセット完了
LOOP1:
LOOP2:
ADDR <- PC   # プログラムカウンタをアドレスレジスタにコピー
ROMを実行する
DATAをデコードする
PC++
なんか実行する
if (出力命令でもHLT命令でもない) goto LOOP1;
出力する
if (HLT命令じゃない) goto LOOP2;

大枠はこんな感じ。あとは命令セットを決めて実装すればよい。Quineを作るだけなので簡単な命令だけでよくて

// 8bitレジスタ  A
// 16bitレジスタ HL
// 16bitレジスタ SAVED_PC
op = read_rom(PC);
if (op & 0x80) {
  if (op & 0x40) {
    oprand, eom = read_rom(HL); // オペランドをフェッチする。HLがメモリの終端を指してたらeom=1
    if (op & 0x20) {
      // INC_HL_AND_JUMP_UNLESS_EOM
      HL++;
      if (!eom) {
        PC = SAVED_PC;
      }
    } else {
      // LD A, (HL)
      A = oprand;
    }
  } else {
    if (op & 0x20) {
      // SHIFT_PUTC 命令
      putchar('0' + (A & 1));
      A >>= 1;
    } else {
      if (op & 0x10) {
        // HLT 命令
        hlt();
      } else {
        // SAVE_PC命令
        SAVED_PC = PC;
      }
    }
  }
} else {
  // PUTC 命令
  putchar(op & 0x7F);
}

というルールでデコードして実行する。オペランドをフェッチするところをCPUに足すのは簡単にできるからこれで完成。

定数文字を出力する命令が一番多いので、7ビット目が立ってない命令はそのまま出力にしてある。

2010-03-06 ferNANDoでQuine(2)

[] やっとおわた 16:37  やっとおわた - 菊やんの雑記帳 を含むブックマーク

kik@shadowgate:~/work/sandbox/fernando$ ruby quinecode.rb >quine.mem
kik@shadowgate:~/work/sandbox/fernando$ ruby quinegen.rb quine.mem >q.nand
kik@shadowgate:~/work/sandbox/fernando$ time ./ferc q.nand >q.nand.out

real    2055m44.082s
user    2054m54.533s
sys     0m7.412s
kik@shadowgate:~/work/sandbox/fernando$ cmp q.nand q.nand.out
kik@shadowgate:~/work/sandbox/fernando$

とりあえず、ここに置いときます。http://as305.dyndns.org/~kik/p/q.nand.gz