はじめに
整数二進対数とは
入力と出力がどちらも整数の、底を$2$とする対数です。式にすると以下の通りです。
$$
f(x) = \lfloor log_{2}{x} \rfloor
$$
これは整数を二進数で表現した際の、一番左に立っているビット位置（LSBを0として数えた場合）に等しいです。なお、対数なので$x=0$を未定義とします。
もう少し正確に表現すると、$x$の定義域は$1$以上の自然数です。
どうして整数二進対数？
整数二進対数は完全平衡二分木の階層と一致します。つまり、配列で表現された完全平衡二分木の要素インデックス$i$($i$はゼロオリジン)から、その要素が属する階層$h$を整数二進対数で得る事ができます。
$$
階層h = \lfloor log_{2}{(i + 1)} \rfloor
$$
完全平衡二分木の階層インデックスとして扱いたいので、なるべく計算速度が速い実装方法を採用したいです。
いくつかの実装と速度比較
本記事は、C言語で本実装する前の性能検証とプロトタイピングを目的として、pythonでいくつかの実装を試します。前提として、入力する整数は符号なし$32$bit整数までの幅を想定しています。
定義通りの実装
import numpy
def intlog2_math(v):
return int(numpy.floor(numpy.log2(v)))
まずは愚直に定義通りの実装です。浮動小数点数を使った数学関数なのでおそらく遅いでしょう。
上位ビットから線形探索
def intlog2_linsearch(v):
bit_mask = 0x80000000
for i in range(31, -1, -1):
if 0 < (v & bit_mask):
break
else:
bit_mask = bit_mask >> 1
return i
要は一番左の立っているビット位置がわかれば良いので、ビットマスクをシフトしつつ上位ビットから順番にビットが立っているかを探索していきます。
単純に考えれば計算量$O(N)$ですが、完全ランダムな値が入力される場合、上位ビット側の方がたった整数の方が下位ビット側しか立っていない整数よりも出現率が高いので、期待値は多分大体$2$回程度の探索回数です。
しかし、逆に実用上は下位ビット側しか使わないような小さな値しか扱わない可能性があるので、実用上はよりひどい性能になると予想します。
二分探索（二分法）
def intlog2_binsearch(v):
if 0 < (v & 0xffff0000):
if 0 < (v & 0xff000000):
if 0 < (v & 0xf0000000):
if 0 < (v & 0xc0000000):
if 0 < (v & 0x80000000):
return 31
else: # 0x40000000
return 30
else: # 0x30000000
if 0 < (v & 0x20000000):
return 29
else: # 0x10000000
return 28
else: # 0x0f000000
if 0 < (v & 0x0c000000):
if 0 < (v & 0x08000000):
return 27
else: # 0x04000000
return 26
else: # 0x03000000
if 0 < (v & 0x02000000):
return 25
else: # 0x01000000
return 24
else: # 0x00ff0000
if 0 < (v & 0x00f00000):
if 0 < (v & 0x00c00000):
if 0 < (v & 0x00800000):
return 23
else: # 0x00400000
return 22
else: # 0x00300000
if 0 < (v & 0x00200000):
return 21
else: # 0x00100000
return 20
else: # 0x000f0000
if 0 < (v & 0x000c0000):
if 0 < (v & 0x00080000):
return 19
else: # 0x00040000
return 18
else: # 0x00030000
if 0 < (v & 0x00020000):
return 17
else: # 0x00010000
return 16
else: # 0x0000ffff
if 0 < (v & 0x0000ff00):
if 0 < (v & 0x0000f000):
if 0 < (v & 0x0000c000):
if 0 < (v & 0x00008000):
return 15
else: # 0x00004000
return 14
else: # 0x00003000
if 0 < (v & 0x00002000):
return 13
else: # 0x00001000
return 12
else: # 0x00000f00
if 0 < (v & 0x00000c00):
if 0 < (v & 0x00000800):
return 11
else: # 0x00000400
return 10
else: # 0x00000300
if 0 < (v & 0x00000200):
return 9
else: # 0x00000100
return 8
else: # 0x000000ff
if 0 < (v & 0x000000f0):
if 0 < (v & 0x000000c0):
if 0 < (v & 0x00000080):
return 7
else: # 0x00000040
return 6
else: # 0x00000030
if 0 < (v & 0x00000020):
return 5
else: # 0x00000010
return 4
else: # 0x0000000f
if 0 < (v & 0x0000000c):
if 0 < (v & 0x00000008):
return 3
else: # 0x00000004
return 2
else: # 0x00000003
if 0 < (v & 0x00000002):
return 1
else: # 0x00000001
return 0
実装がひどい事になっていますが、愚直に上位ビットと下位ビットに分割しながら、ビット位置を二分探索しています。$32$bitを想定しているので全$32$通りのケースを愚直に条件分岐で実装しています。C言語だったらこれでも十分速くなりそうです。探索回数は常に$5$回です。
対数関数の解を二分探索で求めているので、実質的に整数版の二分法だと言えます。
二分探索と8bit幅テーブルのハイブリッド
bit8_table = [0 if v == 0 else intlog2_math(v) for v in range(256)]
def intlog2_8bit_table(v):
if 0 < (v & 0xffff0000):
if 0 < (v & 0xff000000):
return bit8_table[v >> 24] + 24
else: # 0x00ff0000
return bit8_table[v >> 16] + 16
else: # 0x0000ffff
if 0 < (v & 0x0000ff00):
return bit8_table[v >> 8] + 8
else: # 0x000000ff
return bit8_table[v]
途中まで二分探索を行い、探索範囲が$8$bit以下になったらあらかじめ用意していたテーブルを参照して答えを出す実装です。探索回数を$3$回節約できています。C言語の実装であればテーブルサイズは$256$byteで足りるはずです。(使わないですが$0$インデックスも含める想定です。)
とは言え、ROMサイズを消費する点と、場合によってはキャッシュフェッチのペナルティがあるかもしれません。
二分探索と16bit幅テーブルのハイブリッド
bit16_table = [0 if v == 0 else intlog2_math(v) for v in range(65536)]
def intlog2_16bit_table(v):
if 0 < (v & 0xffff0000):
return bit16_table[v >> 16] + 16
else: # 0x0000ffff
return bit16_table[v]
単にテーブル幅を$16$bitまで拡張しただけの実装です。空間計算量を犠牲にして時間計算量を減らす戦略です。と言ってもC言語であれば64KiB程度のROM/RAM消費で済むので、現代のパソコンであれば問題無いレベルです。探索回数を$4$回節約できており、比較$1$回、テーブル引きで二進対数を求める事が出来ます。
二分探索と32bit幅テーブルのハイブリッド?
これを実現できれば定数時間でほぼ最速だと思いますが、単純には$4$GiBも空間計算量を犠牲にするので、いくら現代的なパソコンといえどもあまり現実的では無いと考えます。なので試していません。
性能評価
それぞれの関数について、整数$1 \dots 999999$を入力して全て処理し終えるまでの総時間を計測して比較します。
時間計測コード
if __name__ == "__main__":
import time
n = 1000000
def _time(target_function):
st = time.time()
_ = [target_function(v) for v in range(1, n)]
et = time.time() - st
return (target_function.__name__, et, _)
results = []
results.append(_time(intlog2_math))
results.append(_time(intlog2_linsearch))
results.append(_time(intlog2_binsearch))
results.append(_time(intlog2_8bit_table))
results.append(_time(intlog2_16bit_table))
for v in results:
print(f"|{v[0]:20}|{v[1]:8.6f}|")
時間計測結果
| 実装 | 処理時間（秒） |
|---|---|
| intlog2_math | 2.375886 |
| intlog2_linsearch | 2.029245 |
| intlog2_binsearch | 0.454672 |
| intlog2_8bit_table | 0.326109 |
| intlog2_16bit_table | 0.258130 |
性能測定環境は以下の通りです。
Hardware Overview:
Model Name: MacBook Pro
Model Identifier: MacBookPro13,3
Processor Name: Quad-Core Intel Core i7
Processor Speed: 2.7 GHz
Number of Processors: 1
Total Number of Cores: 4
L2 Cache (per Core): 256 KB
L3 Cache: 8 MB
Hyper-Threading Technology: Enabled
Memory: 16 GB
Boot ROM Version: 428.0.0.0.0
SMC Version (system): 2.38f11
Pythonバージョンです。
Python 3.9.1 (default, Jul 15 2021, 23:56:30)
[Clang 12.0.0 (clang-1200.0.32.28)] on darwin
考察
Pythonだと$16$bit幅テーブルを用意して二分探索の途中からテーブル引きに切り替える実装が最速でした。一番遅いのは愚直に定義通りに実装した場合で、予想通りといえば予想通りでした。
愚直な二分探索はコーディングの労力虚しくそこまで速くなかったのが残念です。でもC言語で実装した場合に賢いコンパイラ様がきっと爆速に最適化してくれそうです。
テーブル引きはやはりキャッシュフェッチのペナルティが実行時間の揺れを大きくしたり、そもそもROM/RAM消費量が増えるので、候補としては愚直な二分探索、8bitテーブルのハイブリッドあたりが実装方針として妥当だと考えます。
おわりに
できればテーブル無しの定数時間計算量のアルゴリズムが欲しいです。探しても見つかりませんでした。
今回の検証コードはここにあります。
追加実装
ここからはコメントを頂いた後にさらに追加調査および実装した内容です。
文字列化して長さをはかる
def intlog2_py_string(v):
return len(bin(v)) - 3
非常にシンプルですが、いざC言語に持って行ったら文字列化するところのオーバーヘッドで遅い実装になりそうなのと、文字列長を数え上げる計算量が$O(N)$になりそうです。
Python組み込み機能のbit_length()を使う
def intlog2_py_bit_length(v):
return v.bit_length() - 1
整数オブジェクトの組み込み機能なのでおそらく高速でしょう。しかし組み込み機能で高速であってもC言語に移植できないと目的を果たせないので、実装を調べてみます。
_Py_bit_length関数がbit_length()の実装です。
#if (defined(__clang__) || defined(__GNUC__))
if (x != 0) {
// __builtin_clzl() is available since GCC 3.4.
// Undefined behavior for x == 0.
return (int)sizeof(unsigned long) * 8 - __builtin_clzl(x);
}
else {
return 0;
}
#elif defined(_MSC_VER)
ClangかGCCの場合は組み込み関数__builtin_clzlを使うように実装されていました。私の実行環境のPythonはClang 12.0.0でビルドした処理系なので、おそらくこの組み込み関数実装が有効になっているものと思います。
この組み込み関数__builtin_clzlはMSB(最上位ビット)から'0'が連続している数を返す関数です。CPUに専用命令があれば、コンパイラが専用命令を使うようにコンパイルしてくれるはずですが、専用命令がなかった場合の実装がどうなるのかよくわかりません。
そこで、組み込み関数を使わない実装の方をPython上で再現しました。
Python組み込み機能のbit_length()のC言語標準機能による実装
BIT_LENGTH_TABLE = [
0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5
]
def intlog2_py_bit_length_on_python(v):
msb = 0;
while (v >= 32):
msb +=6;
v >>= 6;
msb += BIT_LENGTH_TABLE[v]
return msb - 1
C言語の実装をわざわざPythonに移植してみるというのも本末転倒感がすごくありますが、折角なので実装して見て、実装内容や実測値を他の実装と比較してみます。
見たところ下位ビット側から線形探索していきますが、ある程度まとまったビット単位でテーブル参照をして比較回数を減らしています。線形探索なので理屈上は二分探索より計算量は大きくなるはずですが、入力する整数が小さい値に偏っている場合の探索回数は$1〜2$回程で済みそうなので、実用上は速い可能性があります。
処理時間を再計測する
| 実装 | 処理時間（秒） |
|---|---|
| intlog2_math | 2.426732 |
| intlog2_linsearch | 2.063314 |
| intlog2_binsearch | 0.456377 |
| intlog2_8bit_table | 0.333185 |
| intlog2_16bit_table | 0.257784 |
| intlog2_py_string | 0.286416 |
| intlog2_py_bit_length | 0.171990 |
| intlog2_py_bit_length_on_python | 0.527500 |
Python上であればbit_length()を使った実装が最速でした。おそらく組み込み関数実装になっているはずですが、一方で専用命令を持たないアーキテクチャ向けであったり、組み込み関数を持たないコンパイラ向けの実装も考える必要があります。ある程度処理の速さを犠牲にしても、なるべくC言語標準機能（特にC89）だけで実装できるアルゴリズムを採用したいです。
さらに調査
どうやらWikipediaの英語ページには、この辺りのビット操作に関するページが存在するようです。（日本語ページが存在しない……。）
ここの参考文献を辿ると、スタンフォード大学の以下のWebサイトにたどり着きました。
ここにも整数の二進対数を取得する方法について紹介されていたので、それらも実装してみました。
Find the log base 2 of an integer with a lookup table
def intlog2_stanford_with_lookup_table(v):
if (tt := v >> 16):
r = 24 + bit8_table[t] if (t := tt >> 8) else 16 + bit8_table[tt];
else:
r = 8 + bit8_table[t] if (t := v >> 8) else bit8_table[v];
return r
C言語のようなコードで説明されていたため、Pythonで実装するにあたり代入演算子を用いています。この実装は上下半分ずつに分けて２回探索し、その後テーブル参照して探索回数を節約しています。私のアイディアとほぼ同じですが、ビットマスクではなくビットシフトを使っている点が異なります。
と言っても計算量のオーダー的には二分探索に従うアルゴリズムです。
Find the log base 2 of an N-bit integer in O(lg(N)) operations
LG_N_BITMASK = [0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000]
LG_N_SHIFT = [1, 2, 4, 8, 16]
def intlog2_stanford_O_lg_N_operations_1(v):
r = 0
for i in range(4, -1, -1):
if (v & LG_N_BITMASK[i]):
v >>= LG_N_SHIFT[i];
r |= LG_N_SHIFT[i];
return r
このアルゴリズムは純粋な二分探索になっています。私の実装と異なるのは、ビットマスクとビットシフトを併用して予め用意しておくマスク値の種類を少なく抑えている点です。コードサイズはおそらく小さくなりそうですが、ループ時の条件判定が余計な命令数になっていそうで、もはやここまでくるとアルゴリズムというより特定のアーキテクチャ向け最適化の様相を呈してきています。
そもそも「条件分岐が遅いアーキテクチャ向けの実装」なる記述があったりするので、このあたりはPythonで実装して実行速度を測る事にほとんど意味はなさそうです。
条件分岐が遅いアーキテクチャ向け
def intlog2_stanford_O_lg_N_operations_2(v):
r = (v > 0xFFFF) << 4; v >>= r;
shift = (v > 0xFF ) << 3; v >>= shift; r |= shift;
shift = (v > 0xF ) << 2; v >>= shift; r |= shift;
shift = (v > 0x3 ) << 1; v >>= shift; r |= shift;
r |= (v >> 1);
return r
intlog2_stanford_O_lg_N_operations_1の亜種です。どういうアイディアで条件分岐がなくなっているのか全然解読できていないです……。
Find the log base 2 of an N-bit integer in O(lg(N)) operations with multiply and lookup
MultiplyDeBruijnBitPosition = [
0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31
]
def intlog2_stanford_O_lg_N_operations_with_multiply_and_lookup(v):
v |= v >> 1
v |= v >> 2
v |= v >> 4
v |= v >> 8
v |= v >> 16
return MultiplyDeBruijnBitPosition[((v * 0x07C4ACDD) % 0xffffffff) >> 27]
この実装も謎の実装です。右ビットシフトしつつ自身と論理和をとっているのは、一番左の立っているビットを利用して、その右側を全て$1$で埋める操作です。一番左の立っているビットが右側にシフトされて倍々に増えていきます。
そこまでは良いのですが、そのあとのオーバーフローによるラップアラウンド前提の掛け算と、その後のテーブル参照が意味不明です。ちょっと数学的に難しい概念を使っているようで、この実装が例え速いとしても採用に躊躇してしまいます。（原理のわからない実装をコピペで済ませてしまうのはまずい。）
ちなみに、C言語の符号なし整数型はオーバーフローするとラップアラウンドすることが言語規格で定義されており、大体のCPUアーキテクチャにおいても余計な処理無くラップアラウンドをオーバーフローで実現できていたりします。
一方、Pythonでそのような自然なラップアラウンドを実現する方法がわからなかったので、符号なし$32$bit整数の最大値で剰余計算しています。
全部の実装の処理時間を比較
| 実装 | 処理時間（秒） |
|---|---|
| intlog2_math | 2.426732 |
| intlog2_linsearch | 2.063314 |
| intlog2_binsearch | 0.456377 |
| intlog2_8bit_table | 0.333185 |
| intlog2_16bit_table | 0.257784 |
| intlog2_py_string | 0.286416 |
| intlog2_py_bit_length | 0.171990 |
| intlog2_py_bit_length_on_python | 0.527500 |
| intlog2_stanford_with_lookup_table | 0.286593 |
| intlog2_stanford_O_lg_N_operations_1 | 1.197447 |
| intlog2_stanford_O_lg_N_operations_2 | 0.817361 |
| intlog2_stanford_O_lg_N_operations_with_multiply_and_lookup | 0.819450 |
intlog2_binsearchとintlog2_stanford_O_lg_N_operations_1は計算量オーダーとしては同じはずですが、大きな差が生じてしまっています。単純に言語の特性であったりコーディングの最適化の度合いだったりと、何かアンフェアな要素があるのだと思います。そのあたりはPython上で試行錯誤しても仕方の無いところなので深くは追求しません。
同じように、intlog2_8bit_tableとintlog2_stanford_with_lookup_tableも同じアイディアの実装ですが、intlog2_stanford_with_lookup_tableの方が２割近く処理時間が短いです。これもなんとなくコードの最適化の度合い（ビットマスクを使うかビットシフトを使うか等）の差なのかなと予想します。（良くある定数倍遅い／速い実装）
やはり実装アイディアの理解のしやすさ、コードの単純さ、アルゴリズムの計算量オーダーの点から、愚直な二分探索か8bitテーブルとのハイブリッドあたりが無難そうです。
その他の実現手段
C++20で追加されたstd::countl_zeroが、まさにPythonのbit_length()実装で使われている__builtin_clzl組み込み関数と同じ機能です。
C++20で実装するのであれば標準ライブラリ関数であるstd::countl_zeroを使えば良いだけになりそうです。（標準ライブラリは各処理系、アーキテクチャ毎に極限まで最適化されている事を期待しています。）
また、以下のような提案もあるようです。
pow2ヘッダなるものを用意し、$2$の冪乗にまつわる整数処理に特化した機能を追加する提案です。この中に整数のlog2が存在しており、本記事でまさにやりたかった事が、標準ライブラリ機能として提案されています。
C89実装なんかにこだわらずに最初からC++20で実装するつもりでいれば本記事を書く事もなかったように思いますが後の祭りですね。
ここまでくると多項式近似やfloatの仮数部ハックなんかも試してみたくなりますが、脱線して戻ってこれなくなりそうなのでこの辺りで事前調査は一旦終了とします。