えびちゃんの日記

えびちゃん（競プロ）の日記です。

Re: 明日使えないすごいビット演算

たぶん今日も使えないと思うんですけど（名推理）．

タイトルの元ネタは これ．

上のスライドでは，ワードサイズ \(w\) に対して，以下の演算を \(O(\log w)\) time で求める方法が書かれています．

  • 立っている最上位のビットを取り出す：msb
  • 立っているビットの個数を取り出す：popcount

この記事では，このうち msb の添字を \(O(1)\) time で求める方法を紹介してみます． 今後，単に msb と言った場合に添字の方を指すことにします．

ソースコードは一番下にあります．

参考にした資料はいつもの CS166 の スライド．

以下では，ワードサイズの四則演算およびビット演算を一単位時間で計算できることを仮定します*1．すなわち，\(w\) bits の整数に関するこれらの演算を一単位時間で行えます．

また，扱う整数は符号なし整数とし，\(2^w\) を法として考えます（オーバーフローはよしなにするということ）．

基本方針

まず，\(w\) bits の整数を \(\sqrt{w}\) 個の \(\sqrt{w}\) bits の整数たちだと見なしてみます．この各々の小さい整数のことを指してブロックと呼んだりすることにします．

こうすることにより，\(\sqrt{w}\) 個の整数を（一単位時間で）並列処理することができます．

この並列処理によって定数時間でできる演算をいくつか紹介し，それらを組み合わせることで msb を求めます．

補助の演算たち

以下では，説明のために \(w = 16\) としてみますが，平方数ならなんでも適用できるはずです．平方数じゃない場合は適当な定数倍回行うと（ワードサイズより大きい最小の平方数をエミュレートする気持ちで）いい感じになると思います．

log2p1

\(\sqrt{w}\) bits の整数 \(x\) に対して，\(\lfloor\log_2 x\rfloor+1\) を求めます．これは，\(x\) の立っている最上位ビットの（下位ビットから 1-indexed での）添字に対応します．

すなわち，msb(x) は log2p1(x) - 1 に相当します（ただし 0 については適当にどうにかしたりしなかったりします）．

最上位ビットが立っていたら \(\sqrt{w}\) を返してしまいます．

整数 0010 の log2p1 を求めるのを例として挙げながら手順を書いていきます．答えが \(2\) になってくれるとうれしいです．

まず，互いに異なる位置に 1 を配置した整数を \(\sqrt{w}\) 個並べた整数を用意します．

1000 0100 0010 0001

これらの各要素である \(y_i = 2^i\) と \(x\) を比較し，\(x \ge y_i\) である \(i\) の個数であることがわかります．

\(x \ge y_i \iff x-y_i \ge 0\) でありますから，引き算をしてみましょう．

   0010 0010 0010 0010
-) 1000 0100 0010 0001
----------------------
   1001 1110 0000 0001

\(x-y_i \ge 0\) が成り立つというのは（2 の補数表現を考えるなどして）最上位ビットが立たないことに相当します．ここでいう最上位ビットというのは各ブロックごとの最上位ビットを指しています．

最上位ビットが立っている整数があったとき，キャリー（ボロー？）が発生して，それより上側のブロックにも影響が出ますね．具体的には，そのブロックのみを取り出して引き算したときよりも結果が 1 小さくなってしまいます．

なんですが，大きい数を上位側に並べていることにより，引きすぎるぶんには問題ないはずです．1 減ることで最上位ビットが変わるのは（独立に計算した場合の）結果が 0000 になるときですが，今の桁で最上位ビットが 1 になっていることから，そうなることはないため*2．

逆に，小さい数を上位側に並べていると以下のようになり，最上位ビットの結果が異なってしまいます．

   0010 0010 0010 0010
-) 0001 0010 0100 1000
----------------------
   0000 1111 1101 1010

で，結局，こうではなくて，さっきの順序で引き算をして，次のものが得られました．

   1001 1110 0000 0001

ほしいのは立っていない方なので，反転させておきます．

~) 1001 1110 0000 0001
----------------------
   0110 0001 1111 1110

ここで，最上位ビットのみを取り出すのは簡単です．

   0110 0001 1111 1110
&) 1000 1000 1000 1000
-----------------------
   0000 0000 1000 1000

さて，あとはこの 1 の個数を数えればよいです．

作り方から，各ブロックにある 1 は \(0\) 個か \(1\) 個なので，各ブロックの値自体を 1 または 0 にしてしまい，ブロック全体の合計を求めてあげるとよいことになります．

    0000 0000 1000 1000
>>)                   3  // 4-1
-----------------------
    0000 0000 0001 0001

これを入力とする次のサブルーチンを考えましょう．

popcount

値が 0 または 1 であるブロック \(\sqrt{w}\) 個を並べた整数 \(x\) を受け取り，その 1 の個数を返します．

愚直には，次のようにシフトを行い，足すとよさそうです．

                  0000 0000 0001 0001 (x << 0)
             0000 0000 0001 0001      (x << 4)
        0000 0000 0001 0001           (x << 8)
+) 0000 0000 0001 0001                (x << 12)
-------------------------------------
   XXXX XXXX XXXX 0010 0010 0010 0001

XXXX はオーバーフローで無視される部分です．

ところで，これは乗算ですね．シフト幅が掛ける数で立っているビットに対応します．すなわち，0, 4, 8, 12 bit 目が立っている整数を掛けるとよいです．

   0000 0000 0001 0001
*) 0001 0001 0001 0001
----------------------
   0010 0010 0010 0001

同じになっていますね（適宜自分で計算してね）．

興味があるのは一番上のブロックのみですので，マスクをかぶせます．

   0010 0010 0010 0001
&) 1111 0000 0000 0000
----------------------
   0010 0000 0000 0000

あとはシフトして，答えを得ます．

    0010 0000 0000 0000
>>)                  12
-----------------------
    0000 0000 0000 0010

想定通り，\(2\) が得られました，うれしい．

summary

準備段階の最後です．やや複雑なのですが，これは次のような演算です．

\(w\) bits の整数を入力として受け取り，出力は \(\sqrt{w}\) bits の整数です．

入力 \(x\) を \(\sqrt{w}\) bits ごとに区切ります．これの \(i\) ブロック目が 0 であるとき出力 \(y\) の \(i\) bit 目は 0，そうでないとき 1 です．

たとえば，次のような感じです．

summary(0000 0010 1011 0000) = 0110

まず，ブロックの最上位ビットとそれ以外に分けます．

   0000 0010 1011 0000
&) 1000 1000 1000 1000
----------------------
   0000 0000 1000 0000
   0000 0010 1011 0000
&) 0111 0111 0111 0111
----------------------
   0000 0010 0011 0000

「最上位ビットのみが立っているもの」と「最上位ビットは常に立っていないもの」について分けて考えます． 各場合について summary を求め，それらの or を求めるとよいことがわかりますね．

最上位ビットのみが立ちうるもの

先ほどのようにシフトするとよいです．

シフトの結果，次のようになってくれるとうれしいです．

   a000 b000 c000 d000
*) ???? ???? ???? ????
----------------------
   abcd ???? ???? ????

たとえば次のようにするとよいです．

   000d 0000 0000 0000 (x << 9)
   00c0 00d0 0000 0000 (x << 6)
   0b00 0c00 0d00 0000 (x << 3)
   a000 b000 c000 d000 (x << 0)

各桁について，置かれるビットは高々ひとつなので，繰り上がりなどが発生することはありませんね．

   a000 b000 c000 d000
*) 0000 0010 0100 1001
----------------------
   abcd bcd0 cd00 d000

あとは，シフトして下位ビットを消しちゃいましょう．

    abcd bcd0 cd00 d000
>>)                  12
-----------------------
    0000 0000 0000 abcd

最上位ビットは立っていないもの

次のようになっていますね．

0??? 0??? 0??? 0???

とりあえず，最上位ビットのみが立ちうるパターンに帰着させられると楽そうなので，そういう方針で考えてみます． 0000 の場合に最上位ビットが 0，それ以外の場合に 1 になってほしいですね．

ん〜，1000 から引くとその逆になるので，とりあえずそれを求めてからビット反転させてみましょう．0000 から引くと繰り下がりを考えるのが大変そうです．

先ほどの例では，次のようになります．

   1000 1000 1000 1000
-) 0000 0010 0011 0000
----------------------
   1000 0110 0101 1000

反転させます．

~) 1000 0110 0101 1000
----------------------
   0111 1001 1010 0111

マスクをかぶせます．

   0111 1001 1010 0111
&) 1000 1000 1000 1000
----------------------
   0000 1000 1000 0000

最上位ビットのみが立ちうるパターンに帰着できましたので，あとは掛け算などでできます．

これらにより，summary を求めることができました．

access

下位から数えて 0-indexed で \(i\) 番目のブロックの整数を取得したいときはどうしましょうか． 適当にシフトしたりマスクをかぶせるといいですね．

\(2\) 番目のブロックを取得してみましょう．

    aaaa bbbb cccc dddd
>>)                   8  // 2*4
-----------------------
    0000 0000 aaaa bbbb
   0000 0000 aaaa bbbb
&) 0000 0000 0000 1111
----------------------
   0000 0000 0000 bbbb

本題

さて，準備が整いました．

確認すると，目標は，\(w\) bits の整数 x の msb を求めることでした．

これは，次のようにできます．

以下において，ブロックについての msb は log2p1 から求めます．

  • summary(x) の msb を求め，これを \(i\) とする
  • x の \(i\) 番目のブロックを取得し，これを b とする
  • b の msb を求め，これを \(j\) とする

このとき，msb(x) は \(i\sqrt{w}+j\) となります．

ソースコード

64-bit のもののつもりです．

constexpr uintmax_t operator ""_ju(unsigned long long n) {
  return n;
}

uintmax_t popcount(uintmax_t n) {
  constexpr uintmax_t multiplier = 0x0101010101010101_ju;
  return (n * multiplier) >> 56;
}

uintmax_t minimsb(uintmax_t n) {
  if (n >= 0x80) return 7;

  constexpr uintmax_t multiplier = 0x0101010101010101_ju;
  uintmax_t minuend = n * multiplier;

  constexpr uintmax_t subtrahend = 0x8040201008040201_ju;
  uintmax_t rest = minuend - subtrahend;

  constexpr uintmax_t mask = 0x8080808080808080_ju;
  return popcount((~rest & mask) >> 7) - 1;
}

uintmax_t summary(uintmax_t n) {
  constexpr uintmax_t hi_mask = 0x8080808080808080_ju;  
  constexpr uintmax_t lo_mask = 0x7F7F7F7F7F7F7F7F_ju;
  uintmax_t ones = (n | ~(hi_mask - (n & lo_mask))) & hi_mask;
  
  constexpr uintmax_t multiplier = 0x0002040810204081_ju;
  return (ones * multiplier) >> 56;
}

uintmax_t access(uintmax_t n, uintmax_t i) {
  constexpr uintmax_t mask = 0xFF_ju;
  return (n >> (i * 8)) & mask;
}

uintmax_t msb(uintmax_t n) {
  uintmax_t i = minimsb(summary(n));
  uintmax_t j = minimsb(access(n, i));
  return i * 8 + j;
}

定数時間アルゴリズムが速いとは言ってない系のやつですね．

実測

元スライドのものは添字を求めるものではないし，それに対応するようにする気力は今のえびちゃんにはないので，投げ出しました．ごめん．

1 以上の整数を std::mt19937_64 でランダムに \(10^7\) 個生成して，それらの msb の合計を求める時間を測りました*3．

__builtin_clzll に勝てませんでした． ナイーブに上から 1 を見つけるまでなめるのには勝ちました．

方法 時間 (ms)
built-in 113.4
定数時間くん 188.2
ナイーブ 219.6

おわりです．

*1:word-RAM model の仮定にのっていると思います

*2:これで不十分だったらごめん．

*3:合計にしたのは，最適化で関数呼び出しが捨てられたりするのが怖かったからです．