眠い
ひたすら眠い. 起きては寝てを繰り返してる.
ビット演算クイズ回答
問題はビット演算クイズ参照.
昔, 私が試行錯誤して見つけた方法を振り返る.
int calc(int a, int b, int s)
{
const int Q = 1 << s, Q2 = Q * 2, Q3 = Q * 3;
assert(0 <= s && s <= 16 && && 0 <= b && b < a && a < Q * 4);
int n = 0;
for (;;) {
if (a < Q2) {
n = n * 2;
} else if (b >= Q2) {
n = n * 2 + 1;
b -= Q2;
a -= Q2;
} else if (b >= Q && a < Q3) {
b -= Q;
a -= Q;
} else {
break;
}
b = b * 2;
a = a * 2 + 1;
}
return n;
}
まず関数の形を理解する.
sが決まると(a, b)は0 <= b < a < Q * 4の三角形Dの中を動く.
ifの中を見ると
- Dの上半分(a < Q2)のときは動かない.
- Dの右半分(b >= Q2)のときは左上に移動する.
- 真ん中の正方形のときは左上に移動する.
- それ以外は終了.
その後(a, b)は(a * 2 + 1, b * 2)に移動する.
そうやってループを抜けるまでDの中を移動した回数を返す.
テーブル引きをしようと思ってもテーブルサイズが大きそうだ.
次に具体的な数値を眺めてみる.
値は結構小さい. 三角形の斜辺近くに比較的大きな値が並んでいる感じがする. 斜辺の値(a, b = a-1)を取り出して並べてみよう.
a 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 n 0 0 1 0 2 1 3 0 4 2 5 1 6 3 7 0 8 4 9 2
奇数のところは値が半分. 2べきのところは0になっている. それ以外の偶数は?
何ビットシフト(s)したらnになるのかみてみる.
a 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 s 1 2 1 3 1 2 1 4 1 2 1 3 1 2 1 5 1 2 1 3
なんとなく規則性がありそうな, 無いような.
aを2進数表記して下から0が続く個数zを数えてみる.
a 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 1 0 0 1 1 0 0 0 0 1 1 1 1 0 0 0 0 0
0 1 0 1 0 0 1 1 0 0 1 1 0 0 0 0 1
0 1 0 1 0 1 0 1 0 0 1 1 0
0 1 0 1 0
z 0 1 0 2 0 1 0 3 0 1 0 2 0 1 0 4 0 1 0 2
s = z + 1っぽい. つまりaの下から0が続く個数 + 1だけ右シフトするとn.
- bsf(x) = xを2進数展開して下からみて初めて1になる場所
- bsr(x) = xを2進数展開して上からみて初めて1になる場所
とするとcalc(a, a - 1) = a >> (bsf(a) + 1).
b = a - 1以外のbではどうだろう.
頑張って眺めてみるとa ^ bの上からみて初めて1になる場所(bsr)だけずらせばよいようだ.
a >> (bsr(a ^ b) + 1).
しかし, 斜辺がbsfでそれ以外がbsrだと扱いづらい.
bsfとbsrの関係
x ^ (x - 1)が何を意味するかを思い出す.
xを2進数展開して下から何個か0が続き初めて1が立ち, そのあと?な状態として表す.
xは一般的にそのように書ける(連続する0が0個のときを含める).
x = ? ? ? ? 1 0 0 ... 0 x - 1 = ? ? ? ? 0 1 1 ... 1 x^(x-1)= 0 0 0 0 1 1 1 ... 1
xの下からみて初めて1になる場所bsf(x)はx^(x-1)の上からみて初めて1になる場所bsrに等しい. つまり
bsr(x ^ (x - 1)) = bsf(x).
斜辺はb = a - 1なのでbsf(a) = bsr(a ^ (a - 1)) = bsr(a ^ b)と同じ形になった.
結局
int calc2(int a, int b, int s)
{
return a >> (bsr(a ^ b) + 1);
}
という形になった(実際にこれが元の関数に一致する検証は省略 : 追記 @shinhさんによるロジカルな回答 ビット演算クイズを解いた時の話). 興味深い点が二つ.
- calcはsに依らない
- a >> (bsr(a ^ b) + 1) == b >> (bsr(a ^ b) + 1).
さて, ここまで簡単になるとSIMD化したくなる. しかしbsrのSIMD命令は無いのでどうするか.
bsr(x)は上からみて初めて1が立つ場所なのでceil(log2(x))に等しい.
xの浮動小数点数表記はfを1以上2未満の小数, eを整数としてx = f * 2eとかけるので 整数xをfloatに変換して指数部を取り出せばよい(codec-calc-ans.cpp).
static inline int bsr(int x)
{
assert(x > 0);
union {
int i;
float f;
} ui;
ui.f = (float)x;
return (ui.i >> 23) - 127;
}
これならcvtpi2psを使うことで容易にSIMD化できる.