#popcnt

Bit count algorithm

整数 \(i\) を2進数表記したとき、何個のビットが１であるかを数える。 整数 \(i\) が次のように書けたとする。変数が \(n\) ビットであるとき、ビットは全部で \(n\) 個ある。 以下では、整数 \(i\) を2進数表記したとき、左から \(m\) 番目 \((1 \leq m \leq n)\) のビットが \(b_m\) であるとき \[ i = b_1b_2 \cdots b_m \cdots b_n \] と表記する。 整数の2進数表記と対応させると、1の位が \(b_n\) 、 \(2^p\) の位が \(b_{n - p}\) である。 また、普通のIEEEでは整数ビットの個数は2の累乗にすることがほとんどなので、 \(n = 2^c\) とする。

この数列 \({b_i}\) のうち、1である個数を数えることを目的とする。 各桁ごとにビットを調べていくと、 \(n\) 回のループと \(n\) 回のビット演算が必要になる。

ここで工夫を施すと、隣り合ったビットのペア \(n/2\) 個について、 それぞれのペア中に1が何個あるかは3回のビット演算と1回の足し算でまとめて計算できる。 また同様の工夫を \(r\) 回施すと、隣り合ったペアのペア \(n/2^r\) 個それぞれのペアについて3回のビット演算と1回の足し算でまとめて計算できる。 これで必要なビット演算の回数を大きく減らすことを考える。 ビット演算の回数がおおよそ \(3\log (n)\) 回になることが期待される。

ビット列に対する操作

次のような整数で書けるフィルタを考える。

\(f(p)\) :0と1が \(p\) 個ずつ交互に並んだ2進数。たとえば \(f(1) = 010101\cdots\) 、 \(f(3) = 000111000111 \cdots\) 。

続いて、整数 \(i = {b_m}\) と \(f(p) = {f_m}\) に次の演算をする。

\(f(p)\) と \(i\) について、各桁ごとに論理積を取る。 すなわち \[ k_0 = f(p) \wedge i = (f_1 \wedge b_1)(f_2 \wedge b_2) \cdots (f_n \wedge b_n) \]

\(i\) を右に \(p\) 桁ずらし、また各桁ごとに論理和を取る。 すなわち \[ k_1 = f(p) \wedge (i/2^p) = (f_1 \wedge b_{1 + p})(f_2 \wedge b_{2 + p}) \cdots (f_n \wedge b_{n - p}) \]

\(k = k0 + k1\)

Cでは次のようにコーディングされる。

int k0 = f(p) & i; int k1 = f(p) & (i >> p); int k = k0 + k1;

この3つの演算をまとめて \(a(p)\) とする。すなわち、 \(a(p): i \mapsto k\) 。

アルゴリズム

まず最初の \(i\) に \(a(1)\) を作用させてみる。 \[ i = b_1b_2b_3b_4 \cdots b_{n-1}b_n \] \[ k_0 = 0b_20b_4 \cdots 0b_n \] \[ k_1 = 0b_10b_3 \cdots 0b_{n-1} \] \[ k = (0b_1 + 0b_2)(0b_3 + 0b_4) \cdots (0b_{n - 1} + 0b_n) \] つまり、隣り合った \(b_{2m}\) と \(b_{2m + 1}\) のビットを足して、その値を元あった２つのビットにおいている。見方を変えると、 \(2m\) と \(2m + 1\) のうちの１の数をそこに２進数でおいている。 \[ c_j = b_{2m} + b_{2m + 1} \] とし、２桁になるようにビットを０で埋める。とすると、 \(c_m\) は00、01、10のどれかである。これを使うと、 \[ k = c_1c_2 \cdots c_{n/2} \] とかける。これを \(i_1 = k\) としておく。

続いて、 \(i_1\) に \(a(2)\) を作用させる。 \[ f(2) = 001100110011 \cdots \] との論理積を各桁でとること、 \(c_m\) がそれぞれ２桁であることを考えると、 \[ k_0 = 00c_200c_4 \cdots 00c_{n/2} \] \[ k_1 = 00c_100c_3 \cdots 00c_{n/2 - 1} \] \[ k = (00c_1 + 00c_2)(00c_3 + 00c_4) \cdots (00c_{m/2 - 1} + 00c_{m/2}) \] となるが、 \(a(1)\) のときと同様に４桁のビットの塊 \[ d_m = c_{2m} + c_{2m + 1} \] となるように定義すると、 \(d_m\) は隣り合った \(4m\) 、 \(4m + 1\) 、 \(4m + 2\) 、 \(4m + 3\) のうち１の個数となり、0000、0001、0010、0011、0100のどれかである。 \[ i_2 = d_1d_2 \cdots d_{m/4} \] とかける。

結論

なんとなく一般化すると、 \[ a(2^k): i_{k - 1} \mapsto i_k \] \[ i_k = x_1x_2 \cdots x_{m/2^k} \] ここで \(x_m\) は \(2^k\) 桁の0もしくは1の列で、 \(b_{2^km} \cdots b_{2^{k + 1}m - 1}\) までの1の数を２進数表記したもので足りない桁を0で埋めたものである。

このまま \(n/2^k = 1\) つまり \(k = \log_2(n)\) になるまで繰り返すと、 \(i_c (c = \log_2(n))\) は \(i = {b_m}\) 中の1の数となる。ビット演算は \(3c\) 回である。

実用的には \(k

例

たとえば8byte（64bit）整数の場合。各桁をチェックしていくと64回のビット演算と64回のループ足し算となる。

\(i\) は64個の0もしくは1の列 \({b_m}\) とする。これに \(a(1)\) を作用させると、00、01、10のどれかが32個並んだ列 \(i_1 = {c_m}\) ができる。さらに \(a(2)\) を作用させると、

0000

0001

0010

0011

0100

のどれかが16個並んだ列 \(i_2 = {d_m}\) となる。 もういっちょ \(a(4)\) を作用させると、

00000000

00000001

00000010

00000011

00000100

00000101

00000110

00000111

00001000

のどれかが8個並んだ \[ i_3 = e_1e_2e_3e_4e_5e_6e_7e_8 \] となる。このとき \(e_m\) は \(i\) の中の \(8(m - 1)\) ～ \(8m - 1\) 桁の8個のビットのうちの1の個数を2進数で表記したものとなっている。 答えとなる \(i_c\) はこれら \(e_1\) から \(e_8\) を足すか、 さらに \(a(8)\) を作用させ、 \(i_4 = g_1g_2g_3g_4\) にして総和 \(i_c = g_1 + g_2 + g_3 + g_4\) を取る。 前者だとビット演算は17回で足し算が10回、後者だとビット演算16回足し算7回である。

\(a(32)\) まで作用させ \(i_6\) までまとめたとき、 \(n = 2^6\) なので \(i_6 = i_c\) は \(i\) の中の1の個数となるが、ビット演算は18回、足し算は6回である。

これくらいならループを書かずともよい。実際javaのbitcount関数は32ビット整数用であるが、 \(a(2^2)\) までの作用で実装され、ループが使われていない。

Fortran

Fortranで実装すると以下のようになる。

function bitcount(i) result(r) implicit none integer(8),intent(in) :: i integer :: r integer(8) :: i1,i2,i3,i4 ! ruler 1234567890123456 integer(8),parameter :: f1 = z'5555555555555555' integer(8),parameter :: f2 = z'3333333333333333' integer(8),parameter :: f4 = z'0F0F0F0F0F0F0F0F' integer(8),parameter :: f8 = z'00FF00FF00FF00FF' integer(8),parameter :: mr = z'00000000000000FF' ! a1: 01010101010101010101010101010101 ... = z'5555 ... ! a2: 00110011001100110011001100110011 ... = z'3333 ... ! a4: 00001111000011110000111100001111 ... = z'0f0f ... ! result < 64, mr = 256 is large enough i1 = iand(i,f1) + iand(ishft(i,-1),f1) i2 = iand(i1,f2) + iand(ishft(i1,-2),f2) i3 = iand(i2,f4) + iand(ishft(i2,-4),f4) i4 = iand(i3,f8) + iand(ishft(i3,-8),f8) r = i4 + ishft(i4,-16) + ishft(i4,-32) + ishft(i4,-48) r = iand(r,mr) end function bitcount

補足

C言語であれば/usr/lib/gcc/x86_64-linux-gnu/4.x/include/に各種intrin.hがあるので、population count関数をインクルードする。

#define __POPCNT__ #include long long unsigned int c,i; c = _mm_popcnt_u64(i)

でCPU拡張命令を利用したものが使える。 これはCPUが提供する拡張命令を利用しているので、環境によっては自分で関数を定義するよりはるかに速くなる。

Fortranであれば

integer :: n,p,i n = popcnt(i) !整数i中の1の数 p = poppar(i) !整数i中の1の数が偶数なら0、奇数なら1を返す

など各種ビット演算に関する関数が組み込まれている。

Fortranについて、性能調査を行った。1から1000000000まで繰り返し、かかった時間を比較する。

program test implicit none integer(8) :: n,i,c n = 10_8**9 c = 0 do i = 1,n c = c + function_of_bitcount(i) enddo write(*,*)c !最適化でループが消えてしまわないように stop end program test

このコード中のfunction_of_bitcountを、各種関数に置き換えて比較した。 コンパイルは、

gfortran -o hoge -O2 -msse4.2 hoge.f03 time ./hoge

とした。実行環境はAMD A10-6800K、Cygwin。 後ろでニコ動流したりBOINCが動いてたり各種クラウドサービスが動いてたりするのはご愛嬌。

関数 　　real(秒)　　 　 user(秒)　 　　sys(秒)　 実直ループ 142.098 122.538 0.389 a(8)まで 15.768 12.308 0.139 組み込み 3.832 2.495 0.109

実直ループ: 64個すべてのビットについてbtestでチェックし、trueなら1足すという素直な発想の関数。