prime's diary

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

Delta swapとその応用 【ビット演算テクニック Advent Calendar 2016 3日目】

ビット演算

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

www.adventar.org

の三日目の記事です。

ビット列の部分列を入れ替える操作を効率的に行えるDelta swapを紹介します。xorが好きなのでDelta swapも好きなテクニックの一つです。

行える操作

マスクとずらす量を自由に指定できます。

abcdefghijklmnop 元のビット列

0000011000011100 マスク
00000fg0000lmn00 マスクを取ったビット列
00fg0000lmn00000 3ビット左にずらしたビット列

0011000011100000 3ビット左にずらしたマスク
00cd0000ijk00000 マスクを取ったビット列
00000cd0000ijk00 3ビット右にずらしたビット列

ab00e00h000000op 入れ替えない残りのビット列

abfgecdhlmnijkop 結果のビット列

実装

まずは素朴な実装をしてみます。

uint32_t swap_bits_32(uint32_t bits, uint32_t mask, int delta) {
  return ((bits & mask) << delta) | ((bits >> delta) & mask) | (bits & ~(mask | (mask << delta)));
}

それぞれの部分列、そしてそれ以外の場所をそれぞれマスクを取り、シフトしてORで一つにまとめます。(最適化してシフトしてからマスク取ってたりしますがまあいいでしょう)

それがDelta swapでは次のように書けます。

uint32_t delta_swap_32(uint32_t bits, uint32_t mask, int delta) {
  uint32_t x = (bits ^ (bits >> delta)) & mask;
  return bits ^ x ^ (x << delta);
}

だいぶシンプルになりましたね。 ビット演算の演算回数も10回（andnotがあれば9回）から6回に減りました。

応用

Delta swapを複数回行ってビット列を任意に並べ替えることができます。

もちろん1ビットずつ入れ替えれば任意に並べ替えることができますが、 実際にはより効率的に並べ替えることができます。 32ビットなら次のようにdeltaは固定でマスクをうまく設定することで9回のDelta swapで任意の並べ替えを実現できます。(2nビットなら2n-1回(n>0))

bits = delta_swap_32(bits, mask1,  1);
bits = delta_swap_32(bits, mask2,  2);
bits = delta_swap_32(bits, mask3,  4);
bits = delta_swap_32(bits, mask4,  8);
bits = delta_swap_32(bits, mask5, 16);
bits = delta_swap_32(bits, mask6,  8);
bits = delta_swap_32(bits, mask7,  4);
bits = delta_swap_32(bits, mask8,  2);
bits = delta_swap_32(bits, mask9,  1);

証明は省略しますが、偶数ビット目と奇数ビット目に分けて再帰的に入れ替えを行って最後に偶奇を合わせる感じです。

bit matrix

16bit, 64bitのビット列をそれぞれ4x4、8x8のビットの行列(bit matrix)と思うと二次元のデータを効率的に扱えることがあります。

さて、Delta swapを使うと効率的にbit matrixを反転・回転させることができます。 たとえば対角線を軸とした反転は次のように書けます。

uint64_t flipDiagA8H1(uint64_t bits) {
  bits = delta_swap(bits, 0x00000000F0F0F0F0, 28);
  bits = delta_swap(bits, 0x0000CCCC0000CCCC, 14);
  return delta_swap(bits, 0x00AA00AA00AA00AA,  7);
}

他の例はhttps://chessprogramming.wikispaces.com/Flipping+Mirroring+and+Rotatingにあります。

参考文献

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

ビット演算 競技プログラミング

はじめに

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

www.adventar.org

の1日目の記事です。

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

2021/01/03: 高速ゼータ変換/メビウス変換/畳み込みに関する部分を全面的に書き換えました。

表現方法

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) {
  // ここに処理を書く
}

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

部分集合と和に関するテクニック

2021/01/03 高速ゼータ変換/メビウス変換/畳み込みに関する部分を全面的に書き換えました。

  • 部分集合について和を取るもの

  • {\displaystyle{ g(U) = \sum_{T \supset U} f(T)}} （Uを包含するTについて和を取る: sum over supersets）

  • {\displaystyle{ g(U) = \sum_{T \subset U} f(T)}} （Uに包含されるTについて和を取る: sum over subsets）

  • 逆変換、要素数ごとに符号が反転する

  • {\displaystyle{ g(U) = \sum_{T \supset U} (-1)^{|U-T|}f(T)}} （Uを包含するTについて和を取る）

  • {\displaystyle{ g(U) = \sum_{T \subset U} (-1)^{|U-T|}f(T)}} （Uに包含されるTについて和を取る）

これらの計算は、上述の「ある部分集合の部分集合全体を渡る」テクニックにより、O(3N)で計算することができます（計算量の証明は二項定理による）。 しかし、さらに効率よく計算するアルゴリズムが存在します。これらはしばしば高速ゼータ変換/高速メビウス変換と呼ばれています。

あらかじめdp[T]にf(T)の値を入れておくことにします。 1は

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)];
    }
  }
}

2は

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

3は

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)];
    }
  }
}

4は

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

というコードでO(N2N)でできます。((T >> i) & 1) が0か1か、Tにビットを付け足すのか取り除くのか、演算が+なのか-なのかが違うだけです。 また、ここでは加算演算を利用しましたが、可換で結合的で、単位元がある（つまり、可換モノイド）であれば1と2が、さらに逆元が存在（つまり、可換群（アーベル群））すれば3,4を適用できます。

畳み込み

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

fm = ranked_sum_over_subsets(f);
gm = ranked_sum_over_subsets(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 = inverse_ranked_sum_over_subsets(conv);

ランクkのsum over subsetsは{\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) == 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) == 1) {
        dp[k][T] -= dp[k][T ^ (1 << i)];
      }
    }
  }
}

でdp[popcount(T)][T]に結果が代入されるように計算できる。 また、ここでは加算と乗算を用いましたが、他の可換環でも同様の計算を行うことができます。 (max, +)を用いるトロピカル半環の場合は、元の集合が十分小さい場合にはexp(x)(底は適切に取る)に対して(+,×)の畳み込みを行えばできます。

参考文献

さいごに

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

www.adventar.org

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

CUDAで再帰するとき

GPGPU

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つです。C89でC89コンパイラを頒布するのです！

出落ち的アイデアですが、もう二度と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の決勝に出たけど記録を残せなかった【人権はない】

ISUCON

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

primenumber.hatenadiary.jp

を参照。

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

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

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

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

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