prime's diary

そすうの日々を垂れ流しちゃうやつだよ

ビット列による部分集合表現 【ビット演算テクニック Advent Calendar 2016 1日目】

はじめに

この記事はビット演算テクニック Advent Calendar 2016

www.adventar.org

の1日目の記事です。

この記事ではビット演算を使って有限集合の部分集合を表現、操作するテクニックを紹介します。 本文中で説明を省略した部分については、詳しい説明のある参考文献を紹介しておきます。

表現方法

S={A,B,C,D}という4要素の集合があったとします。 各要素を次のように2進数の数字と対応付けます。 A: 0001 B: 0010 C: 0100 D: 1000 こうすると、Sの部分集合は存在する要素のビットごとの論理和を取ることで4ビットのビット列として一意に表すことができ、例えば {A,B,C}: 0111 {B,D}: 1010 のようになります。

以後、部分集合とそのビット列による表現をしばしば同一視します。

基本操作

空集合、全体集合

それぞれ0000, 1111が対応します。

和集合

和集合はビットごとの論理和(|)で書くことができます。

T = {A,B,C}: 0111, U = {B,D}: 1010

T∪U = {A,B,C,D}: 1111 = 0111 | 1010

共通部分

共通部分はビットごとの論理積(&)で書くことができます。

T = {A,B,C}: 0111, U = {B,D}: 1010

T∩U = {B}: 0010 = 0111 & 1010

補集合

補集合はビット反転(と必要に応じて&でマスクを取る操作)で書くことができます

T = {A,B,C}: 0111

Tc = {D}: 1000 = ~0111

(上位のビットが無視できない場合は 1000 = ~0111 & 1111)

部分集合全体を渡る

numを全体集合の要素数とすると、

for (int i = 0; i < (1 << num); ++i) {
  //ここに処理を書く
}

(本当はシフト演算子のところの括弧は不要ですが、見やすさの為付けています)

包含関係を判定する

補集合と共通部分を持つかどうかで判定できるので、 TがUに含まれるかどうかは~U&T==0で書けます。

素数を数える

1になっているビットの数を数える操作(population count, popcount)が相当します。 最近のIntelのCPUならずばりPOPCNT命令がありますし、そうでない場合も高速に処理する方法が知られています。 32要素の集合Sの部分集合Tの要素数は、

T = (T & 0x55555555) + ((T >>  1) & 0x55555555);
T = (T & 0x33333333) + ((T >>  2) & 0x33333333);
T = (T & 0x0F0F0F0F) + ((T >>  4) & 0x0F0F0F0F);
T = (T & 0x00FF00FF) + ((T >>  8) & 0x00FF00FF);
T = (T & 0x0000FFFF) + ((T >> 16) & 0x0000FFFF);

で計算できます。 ここで、 0x5=0b0101 0x3=0b0011 となっていて、 n行目で2nビットごとのpopcountを計算しています。 さらに最適化して、

T -= (T >> 1) & 0x55555555;
T = (T & 0x33333333) + ((T >> 2) & 0x33333333);
T = (T + (T >> 4)) & 0x0F0F0F0F;
T = (T * 0x01010101) >> 24;

と書くこともできます。

応用操作

ある部分集合の部分集合全体を渡る

T⊃UなるUを全部列挙します。 全体集合の部分集合全体を渡って、Tに含まれているかを判定すればできますが、要素数が少ない場合は無駄が多いですね。 実は、次のようにすることでTの要素数をNとしたとき2N回でループを抜けることができます。

for (int U = (1 << num) - 1; U >= 0; --U) {
  U &= T;
  // ここに処理を書く
}

UをTでマスクしているのがポイントです。 1を引く(デクリメントする)という操作は1になっている一番下位のビットを0にして、それより下のビットをすべて1にする、という操作なので(0のときは除く)、結果的にもれなくかぶりなく列挙できます。

ある集合を部分集合に持つ集合全体を渡る

T⊂UなるUを全部列挙します。

for (int i = T; i < (1 << num); i=(i+1)|T) {
  // ここに処理を書く
}

仕組みは部分集合列挙とほとんど同じです。

素数kの部分集合全体を渡る

int x, y;
for (int T = (1 << k) - 1; T < (1 << N); x = T & -T, y = T + x, T = (((T & ~y) / x) >> 1) | y) {
  // ここに処理を書く
}

かなりトリッキーですね。 詳しくは参考文献の解説を読んでください。

高速ゼータ変換

Sの部分集合を引数に取る関数fがあったとして、 {\displaystyle{ g(U) = \sum_{T \supset U} f(T)}} をSの各部分集合Uに対して求めたいとします。 各Uについて個別に求めれば二項定理により全体の計算量はO(3N)となりますが、Tの取りうる数が2N通りしかないので無駄が多いです。 これを効率よく計算する高速ゼータ変換というアルゴリズムがあって、 あらかじめdp[T]にf(T)の値を入れておくと、

for (int i = 0; i < N; ++i) {
  for (int T = 0; T < (1 << N); ++T) {
    if (((T >> i) & 1) == 0) {
      dp[T] += dp[T | (1 << i)];
    }
  }
}

というコードでO(N2N)でできます。

高速メビウス変換・逆変換

Sの部分集合を引数に取る関数fがあったとして、 [tex:{\displaystyle{ g(U) = \sum{T \subset U} f(T)}}] をメビウス変換(ゼータ変換と包含の向きが逆)、 [tex:{\displaystyle{ g(U) = \sum{T \subset U} (-1)^{|U-T|}f(T)}}] をメビウス逆変換といいます。(集合の引き算は差集合。) その名の通りメビウス変換して逆変換すると元に戻ります。

for (int i = 0; i < N; ++i) {
  for (int T = 0; T < (1 << N); ++T) {
    if ((T >> i) & 1) {
      dp[T] += dp[T ^ (1 << i)];
    }
  }
}

逆変換は差集合の要素数に依存する符号があって複雑そうですが、実はこちらも似たようなコードでできます。(符号が逆になっただけ)

for (int i = 0; i < N; ++i) {
  for (int T = 0; T < (1 << N); ++T) {
    if ((T >> i) & 1) {
      dp[T] -= dp[T ^ (1 << i)];
    }
  }
}

畳み込み

Sの部分集合を引数に取る関数f,gがあったとして、 {\displaystyle{ (f \circ g)(U) = \sum_{T \subset U} f(T)g(U-T)}} をSの各部分集合Uに対して求めたいとします。 (ランク付きの)高速メビウス変換・逆変換を用いると計算できます(多項式乗算におけるFFT/逆FFTと似た理屈)

fm = ranked_mobius(f);
gm = ranked_mobius(g);
for (int k = 0; k <= N; ++k) {
  for (int i = 0; i < (1 << N); ++i) {
    for (int j = 0; j <= k; ++j) {
      conv[k][i] += fm[j][i]*gm[k-j][i];
    }
  }
}
result = ranked_mobius_inv(conv);

ランクkのメビウス変換は{\displaystyle{ \hat{f}(k,U) = \sum_{T \subset U, |T| = k} f(T)}}

実際のコードはこんな感じ(dp[popcount(T)][T]にはf(T)の値を入れておく)

for (int k = 0; k <= N; ++k) {
  for (int i = 0; i < N; ++i) {
    for (int T = 0; T < (1 << N); ++T) {
      if ((T >> i) & 1) {
        dp[k][T] += dp[k][T ^ (1 << i)];
      }
    }
  }
}

メビウス逆変換は{\displaystyle{ f(U) = \sum_{T \subset U} (-1)^{|U-T|} \hat{f}(|U|,T)}}

dp[k][T]{\hat{f}(k,T)}の値を入れておくと、

for (int k = 0; k <= N; ++k) {
  for (int i = 0; i < N; ++i) {
    for (int T = 0; T < (1 << N); ++T) {
      if ((T >> i) & 1) {
        dp[k][T] -= dp[k][T ^ (1 << i)];
      }
    }
  }
}

でdp[popcount(T)][T]に結果が代入されるように計算できる。(実は畳み込みしなくてももともと結果が入っている)

参考文献

さいごに

ビット演算テクニック Advent Calendar 2016

www.adventar.org

みんな書いて!おねがい!

CUDAで再帰するとき

CUDAではDynamic Parallelism(GPU側から新たなカーネルを起動する)以外に普通の(シリアル実行の)再帰が可能です(Compute Capability 2.0以上)。 しかし、再帰的に関数を呼び出す場合、必要なスタックサイズが静的に割り出せないことがあり、このときはデフォルトのスタックサイズ(私の手元の環境では1024bytes/thread)が使用されます。 当然この状態でスタックサイズの制限を超えてスタックを消費すると不幸な目に合います(何も言わずにカーネルの実行が終了してしまうなど)。

解決方法はカーネル呼び出し前にホスト側でスタックサイズを指定してやること。

例えば1スレッドあたり4096バイトのスタックを割り当てたい場合は、ホスト側でカーネル呼び出し前に

cudaDeviceSetLimit(cudaLimitStackSize, 4096);

としてやればよい。

コミックマーケット89記念!ANSI C89コンパイラを作る【KMCアドベントカレンダー12日目】

この記事はKMC Advent Calendar 2015 - Adventar12日目の記事です。

宣伝

部誌、自作ゲーム、そして今回の本題であるコンパイラソースコードを置く予定です。 3日目(12月31日)東地区 "モ" 42bです!よろしくお願いします!!*1 

本題

今度のコミックマーケットは89回目です。コミックマーケット89は通称C89と呼ばれています。 C89といえばC言語の最初の標準規格ANSI-C:1989(通称C89, ANSI C89など)ですね!

こうなったらやることは1つです。C89C89コンパイラを頒布するのです!

出落ち的アイデアですが、もう二度とC89はやってこないので、このチャンスを逃すわけには行かないと部内でC89コンパイラを作るプロジェクトを立てました。

話し合いの結果コンパイラの名称はkmc89になりました。リポジトリはこちらになります。

github.com

規格書の入手

現行のC言語の標準規格ISO/IEC 9899:2011(通称C11)の規格書はISOのウェブサイトからPDFを購入できます。 一つ前のC言語の標準規格ISO/IEC 9899:1999(通称C99)は和訳(JIS X3010:2003)をJISCのウェブサイトから閲覧することが出来ます。

しかし、C89の規格書は上記の方法では購入・閲覧等ができません。 そこで、規格書のハードコピーをhttp://infostore.saiglobal.com/store/Details.aspx?productID=434093から注文して入手しました。

これは厳密にはANSI-C:1989ではなくISO/IEC 9899:1990(通称C90)の規格書ですが内容としては同じものです。*2

コンパイラを作る

コンパイラを作ると言っても、あまり時間はありません。 C89のすべてを実装するのは流石に無理があるということで、標準ライブラリ、リンカ、プリプロセッサ等は実装せず(既存のリンカ等を用いる)、狭義のコンパイラに絞って実装することにし、オブジェクトコードとしてLLVM IRを生成することにしました。

実際には規格書を注文してから届くまでにかなり時間がかかった(再送もあったので2ヶ月弱)ので届く前から少しずつ書き始めています。

最終的にセルフホスティングを目指したいのでコンパイラもC89で記述することになり、それにあたってC言語で貧弱な部分をライブラリとして実装する作業が行われています。 C++のstd::vectorやstd::stringに当たるものをC言語で実装しようとしています。KMCアドベントカレンダー9日目の記事template in C - KMC活動ブログも参照してください。

ある程度そろえばあとはflex、bison、LLVM等の便利なライブラリを使うことで12/31のコミックマーケット3日目に間に合うように実装できるはず!!ということで現在鋭意実装中です!応援よろしくお願いします!

このあたりのコンパイラの仕組みの話は今度のKMCの部誌、独習KMC vol.8に少し書いてあるので気になる人はぜひ読んでください。

今後の展望

今回はコンパイラのコア部分に絞って実装することでなんとか間に合わせる作戦ですが、C89の次のコミケはC90ですね。この時までにC89≒C90のコンパイラの周りの部分の実装を進めて完全版を配布したいと考えています。*3

また、C90の部誌ではコンパイラの実装や開発環境等についての苦労話や新たに得た知見なども書きたいと思っています。ご期待ください。

では早速この記事を書き上げたらコンパイラの実装を進めることにしたいと思います。応援よろしくお願いします!

github.com

コミケは3日目(12月31日)東地区 "モ" 42bです!よろしくお願いします!!

絶対C89にC89コンパイラを間に合わせます!間に合わなかったら木の下に埋めてもらって構わないよ!

次回予告

明日のKMC Advent Calendar 2015 - Adventar 13日目の記事はnonamea774(nona7)さんによる「[データ削除済]」です。お楽しみに。

www.adventar.org

*1:残念ながら当日私はいません

*2:C89の文章に章立てを追加したものになっている

*3:C99の時にC99コンパイラを作るかどうかは未定です。

高速文字列処理ライブラリを作った

この記事はポエムアドベントカレンダー4日目の記事です。

www.adventar.org

大量の文字列データを扱うことの多くなった現代において、文字列処理ライブラリの高速化は重要である。 しかしながら、個人レベルで汎用的かつ高速な文字列処理ライブラリを作成することは難しい。 今回は汎用性を少し下げることにより圧倒的な高速化をした文字列処理ライブラリ「A」を制作した。

ソースコード

gist.github.com

ライブラリの仕様

文字列の制約

  • すべての文字がAで構成されていること

制約を満たす文字列の例

AAAAAA
AAAAAAAAAA

制約を満たさない文字列の例

aaa
ABCDE

制約を満たさない文字列を用いた場合、正しくない結果を得る可能性がある。

利用方法

ライブラリ中では専用の効率的なデータ構造により文字列を管理する。そのため、利用するにはstd::stringから変換処理を行う必要がある。

A a(std::string(100, 'A'));

専用のデータ構造からstd::stringに変換するにはto_cppstringメンバ関数を用いる。

実装されているアルゴリズム

検索 find(A_1, A_2)

文字列A_1の部分文字列に正規表現A_2に一致するものが存在すればその位置(複数あればそのうち最も先頭に近いもの)を返す。 計算量:O(1)

置換 replace(A_1, A_2, A_3)

文字列A_1の部分文字列に正規表現A_2に一致するものが存在すればそれをすべてA_3で置き換える。 計算量:O(1)

編集距離 edit_distance(A_1, A_2)

2つの文字列の編集距離を計算する。 計算量:O(1)

最長共通部分列 longest_common_subsequence(A_1, A_2)

2つの文字列の共通部分列で最も長いものを返す。 計算量:O(1)

最長回文 longest_palindrome(A)

回文になっている部分文字列のうち最も長いものとその始点を返す。 計算量:O(1)

今後の課題

  • さらなるライブラリの充実
  • Bのみからなる文字列を高速に扱えるライブラリの開発

RubyからGoの関数をつかわなくても再帰をやめなくてもアルゴリズムを改善する→はやい

qiita.com

mattn.kaoriya.net

という記事を読んだのでアルゴリズムを改善して速くしました。

いわゆる行列累乗のテクニックを使うと(乗算部分を除けば)N番目のフィボナッチ数は{O(\log N)}で求まります。 実際にはフィボナッチ数は指数的に増大するので乗算にかかる時間が支配的になるのでもっと遅くなります(Karatsuba法なら{O(N^{1.585})}高速フーリエ変換を使えば{O(N \log N)})。 単純な足し算のループでは足し算の桁数を考えると { O(N^{2}) } なのでそれよりは随分と速くなります。

40番目ぐらいでは実行時間はほとんど0なので100万番目のフィボナッチ数を求めてみます。

こちらが元のソースコード

def fib(n)
  f0, f1, v = 0, 1, 0
  n.times do |x|
    if x <= 1
      v = x
    else
      v = f0 + f1
      f0, f1 = f1, v
    end
  end
  return f0 + f1
end

puts fib(1000000)

実行してみます。

$ ruby fib_naive.rb > /dev/null
ruby fib_naive.rb > /dev/null  9.63s user 0.87s system 100% cpu 10.499 total

行列累乗を用いたソースコード。行列を累乗するところで再帰を使っています。

def mat_mul(lhs, rhs)
  res = [[0,0],[0,0]]
  2.times {|i| 2.times {|j| 2.times {|k|
    res[i][j] += lhs[i][k] * rhs[k][j]
  }}}
  return res
end

def mat_pow(mat, index)
  if index == 1
    return mat
  else
    half = mat_pow(mat, index / 2)
    if (index % 2) == 0
      return mat_mul(half, half)
    else
      return mat_mul(half, mat_mul(half, mat))
    end
  end
end

def fib(n)
  return mat_pow([[1,1],[1,0]], n-1)[0][0]
end

puts fib(1000000)
$ ruby fib.rb > /dev/null
ruby fib.rb > /dev/null  0.10s user 0.01s system 97% cpu 0.106 total

さらに100倍ぐらい速くなりました!すごい!

【failに】ISUCON5の決勝に出たけど記録を残せなかった【人権はない】

お金が欲しいので(直球)、優勝賞金100万円のISUCONの決勝に出ました。予選については

primenumber.hatenadiary.jp

を参照。

当日起きれるか不安だったので前日入りしてカプセルホテルに泊まった。 人生初のカプセルホテルだったが、あまり寝られなかった。お腹の調子が悪かったのもあると思う。 少し早めに出てゆったり朝食を取って会場入りした。

コンテスト始まってとりあえず問題の内容を理解したりベンチ回してみたり実装を読んだりプロファイルを取ったりしたりしていた。 なぜかrack-lineprofがSEGVかなんかで落ちるのが辛かった。 とりあえず手元で実験して/dataが外部のAPIを叩くところが遅いことはわかって、どうにかしてキャッシュ出来ないかとか並列にリクエスト送れないかとかあれこれやっていた。 しかし、15時頃に一番遅いtenkiのAPIが単純にキャッシュしただけではダメになったので、八方塞がりになってしまって「なんか自分の知らないすごいからくりを使うとうまく行くみたいなやつでは」という気分になっていた。

結局15:30ぐらいまで全然わからず1000点ぐらいで止まっていたが、そこからどうせ殆どの時間外部からのレスポンスを待って暇しているのなら、すごい並列度を上げても捌けるのでは? ということに気づいてnginxとunicornのworkerを20, 40とかに増やしてみたら8000点ぐらいに乗って、なるほどなぁと思ってそのあとworkerの数を微調整して、並列にリクエストを投げる(ただし、httpsAPIを同じトークンのものを2つ同時に投げるとToo many connectionsとかで怒られるのでそこだけ同期を取って実行した)ようにしたりして13000点まで上がった。 「外部APIのレスポンスが遅い時は、並列度を上げてスループットを上げる」という戦略、私には自明ではなかったけど、多くの他のチームは早いうちから1万点を超えていたのでこういうことはWebサービスを作る人々にとっては自明なことで、この辺が経験の差なのだろうと思った。

この時点で1台しか使っていなかったので頑張って3台で動かそうとしていたけど何故かうまく行かず、最後もなんかうまく行かないかと試しながらやっていたら時間切れ、一度も成功していない3台に分散する方で提出したので当然failして0点だった。 学生賞は1位しか無いのでall or nothingの精神で賭けに出たが、結局全チームfailだったので「安全策で行っていればなぁ」と後悔が残るけど後から言っても仕方無い。

多くの知見を得られた、よいコンテストでした。運営の皆さん、お疲れ様でした。

ISUCON5の予選に出て学生枠1位で通過した(チーム名: 古典論理の犬)

お金が欲しいので(直球)、優勝賞金100万円のISUCONの予選に参加して、13609点、二日目全体で11位(?)学生枠1位で通過した(学生枠で登録した中では一般枠で通過したチーム学生自治についで2位)。

isucon.net

チーム

チーム名: 古典論理の犬

チームメンバー: @lastcat_ @nonylene @_primenumber

チーム名の由来

KMC内でISUCONに出るぞ!!ってなって6人集まったのでnotaで働いてる人(@nonamea774 @pastak @uiureo, チーム学生自治)とそうでない人でそれぞれチームを作った。 not notaなのでnot not a ⇒ aってなって二重否定の除去マンだ、古典論理の犬だ!!ってなって決まった。

流れ

予選まで

8月の末ぐらいからチームで何回か練習をした。

最初にISUCON1からやろうとしてAMIとかもないし流石に古すぎてそもそもアプリ起動まで辿りつけなかったので失敗。結局練習で物理的に集まれたのはこの一回だけだった。

反省を活かしてそれ以降はISUCON3/ISUCON4の予選問題を解いたり、各種テクニックを各自勉強したりしてた。

ISUCON4の予選問題を解いた時は実装をRubyからGoにしただけでかなりスコア上がったので、Goの勉強をしたり、GCPになれるためにGCP上で個人練習したりしながら予選本番を迎えた。

当日

予選二日目の方に参加した。一日目は開始が1時間遅れたりidobataが落ちたりして大変そうだったけど、二日目はそこまで大きな問題はなかった。 学生枠で出ることを考えると一日目と2日目を合わせた順位だけで本戦出場が決まるので、予選の日程によって参加者数とかが違うから有利不利がないし、1日多く練習出来たり運営の人たちが慣れたりして安定する二日目のほうが有利な気はする。

予選開始するとなんとGoの参考実装にバグがあるという話で、しかもコードを見ると800行弱もあって(Rubyだと400行弱)読む気を失ったので早々にRuby実装で行くことに決定。

前回の予選までと違って手元でベンチを回せない、チームあたりベンチキューに積めるのは1つまで、キューの消化が最初はかなり遅かったのもあり結構手間取る。 見た感じとにかく/が遅く、1クエリに2.5sかかってるSQLクエリとかあっておうおうみたいな感じ。 しかし、オンメモリにするにはあまりにも巨大だし、メモリもそんな無いし、真面目にクエリいい感じにしたいけどむずいなこれ…となってなかなかスコアが上がらない。

そうこうしてるうちに1時ぐらいになって、今のままだと絶望だし気分転換にみんなでコンビニ行くか、となって昼飯を書いにコンビニに行く。 ここで3人でいろいろ話して、全員で問題点や疑問点を共有したり出来たのは良かった。 もちろんそれ以前からクエリが遅いとか話してたけど、3人でちゃんと具体的にどういうクエリがどういう理由で遅いのか、今のところこういう解決策がありそう、とかをしっかり共有できたのが大きな前進だった。

いろいろ解決策はありそうだけど実装力がなくてなかなか進まず、16:30頃まではスコアが2000程度で絶望、って感じだったけど、そこからようやくfriendを取ってくるクエリとかも改善して、必要なところにインデックス貼ったりして3000点台→10000点台まで上がった。

その後は細かい調整をして11000点台になったあと、再起動試験をしてクエリキャッシュを有効にしたりカーネルパラメータをいじって13000点台になったところであと3分とかになって、まあ流石に学生枠ではいけるでしょ、お疲れ様でした、と言う感じになった。

その後はチーム学生自治、チームCAMPHOR-の人たちとお疲れ様でした打ち上げをした。

自分のやった作業

  • 都度スロークエリログを眺める
  • 自分の投稿へのコメントを取ってくるクエリはとりあえず自分の投稿のentry_id全部をサブクエリで取ってきてWHERE INで取ってくる実装に変更
  • footprintsのuser_id、relationsのone, anotherのそれぞれにインデックスを張った(oneに張る必要はなかった気がする、あとanotherは見なくてよかったというのを終わってから知った)
  • commentsのentry_id, created_atで複合インデックスを張ったけど効果があったかは不明
  • 「あなたの友達の日記エントリ」を取ってくる実装をfriend全部を取ってきたのを使いまわしてWHERE INの中に突っ込んで直接取ってこれるようにした
  • 「あなたの友達のコメント」のほうは自分に見れない記事へのコメントを弾くのをクエリに入れる方法を考える時間がなかったので、友達のコメント20件取ってきてそのうち自分の見れるものをRuby側で10件取ってくる場当たり的な実装になってしまった(もともとコメント1000件取ってきてRuby側で条件に合う10件取ってくるみたいなのだったので、同じようなもんだし妥協ということにした)

反省

  • 改善策が出てから実装できるまでにかなり時間がかかって、まともなスコアが出るのがかなり遅くなった
    • 実装力がまだまだ足りない、MySQLとかにももっと慣れないといけない
  • 結局Rubyでプロファイル取れてないので、改善すべき点がまだ結構残っている
  • 実はビューの方にもSQLクエリが書いてあったらしいけどビューとか全く見てなかった
  • 問題点とか疑問点はしっかりチームメンバーで共有すべき

感想

ボリュームがあるし面白い問題だった。 チームとしては結果として本選に行けたけど本選でまともに戦うにはまだまだ足りないのでちゃんと練習したい。 物理的に3人でホワイトボードのあるところに集まるとかなり捗ることがわかったので今後チーム練習するときはちゃんと物理的に集まって練習したい。

追記

本番で使ったリポジトリを公開しました。

github.com

mysqlの設定とかは置いてないです