プログラミング
sr_8e
、


　FFTする時に使うビット逆転インクリメントについて複数の実装を説明しつつ比較してみる。
　(一応断っておくと、本記事内ではビット逆転インクリメントと呼ぶものの、これが一般に普及している名称ではなさそうである。需要があまりないので統一されていない)

・ビット逆転インクリメントってなーに

　ビットの並びが前後逆転しているようなインクリメントのことである。当ブログでもFFTの記事では一度出てきたことがある。
　具体的に書けば、普通のインクリメントだと 0000 -> 0001 -> 0010 -> 0011 -> 0100 ... と進むところを、ビットの並びの前後をひっくり返して 0000 -> 1000 -> 0100 -> 1100 -> 0010 ... と進んでいくようなインクリメントである。

　以下、一般に \(N\) ビットの数を考える。コードはPythonで書いてあるが、一般のビット演算が定義されている言語であればほぼ同じように書けると思う。

・naiveな実装

　まずは凄まじく単純な実装を考えてみよう。
　下から \(j\) 番目のビットを見て、1が立っていたら上から \(j\) 番目(=下から \(N-j-1\) 番目)に1を詰めるという、サルでも思いつきそうな実装である。

def bit_reverse_naive(bits):
    table = [0 for _ in range(1 << bits)]
    for i in range(1 << bits):
        s = 0
        for j in range(bits):
            if bits > 2 * j:
                s += (i & (1 << j)) << (bits - 2 * j - 1)
            else:
                s += (i & (1 << j)) >> (2 * j - bits + 1)
        table[i] = s
    return table


def bit_reverse_naive2(bits):
    table = [0 for _ in range(1 << bits)]
    for i in range(1 << bits):
        s = 0
        for j in range(bits):
            if i & (1 << j):
                s += 1 << (bits - j - 1)
        table[i] = s
    return table

　適当に2つほど実装のバリエーションを考えてみた。特に説明もいらないだろう。計算量は(配列の要素の数 \(M\equiv2^N\) について)どちらも \(O(M\log M)\) である。

・再帰的な実装

　次は再帰的な方法でビット逆転列のテーブルを構成する方法である。2進数では(というか、位取り記数法であれば何進数でもそうなのだが)、すべての桁はずっとそれぞれ同じパターンを繰り返している。
　例えば一番下の桁は 0 -> 1 -> 0 -> 1 ->... を繰り返しているし、その次の桁は 0 -> 0 -> 1 -> 1 ->...を繰り返している。だから、\(N-1\) 桁まで (\(2^{N-1}\) 以下) のテーブルがすべて出来上がっているなら、そのテーブルの \(N\) 桁目にすべて1を立てたものを続けてやるだけで \(2^N\) 以下のテーブルが作れることになる。
具体例を示すと、例えば2桁までのテーブル 00 -> 10 -> 01 -> 11 が既に出来上がっているとすると、それぞれの数字の右端に1を立てた 001 -> 101 -> 011 -> 111 をくっつけてやれば、000 -> 100 -> 010 -> 110 -> 001 -> 101 -> 011 -> 111 となって3桁までのテーブルができるといった寸法である。

実装を示す。

def bit_reverse_recursive(bits):
    table = [0 for _ in range(1 << bits)]
    for i in range(bits):
        for j in range(1 << i):
            table[(1 << i) + j] = table[j] + (1 << (bits - i - 1))
    return table

計算量は \(O(M)\) である。事前に \(2^N\) 個のすべての数についてビット逆転列を計算する必要があるときに最適であろう。

・XORを用いた実装

　2進数の桁上がりについて考えよう。例えば 1011 という数に1を足すと 1100 になる。1001 に1を足すと 1010 になる。1000 に1を足すと 1001 になる。これらの例から、足す前の数の末尾に1が \(k\) 個 (\(k\geq0\)) 並ぶとき、末尾から \(k+1\) 個の数字が反転する(01...1 が 10...0になる)ことがわかる。したがって足した後の数の末尾には0が \(k\) 個並ぶことになる。
　したがって、ビット逆転列においては、末尾でなく先頭から数えて \(k+1\) 個の数字を反転させればよい。これは1が先頭から \(k+1\) 個並び、残りがすべて0であるような数(ビットマスク)とのXORをとることで実現できる(XORの性質について詳しくない場合は、各自で調べることをお勧めする。面白いので)。
　では「1が先頭から \(k+1\) 個並び、残りがすべて0であるような数=ビットマスク」をどのように求めればよいだろうか。これは　\(2^N-2^{N-k-1}\) であることは簡単に確かめられよう。例えば、\(N=8, k=3\)だったとすると$$2^8-2^{8-3-1}=100000000_{(2)}-00010000_{(2)}=11110000_{(2)}$$で確かに先頭から1が(3+1)個立った数が求められた。

　また、\(2^{N-k-1}=2^{N-1}/2^k\) だから、\(2^k\), すなわち末尾から \(k\) 個が0で、\(k+1\) 番目に1が立つような数さえわかればビットマスクを計算できる。

　では \(2^k\) をどうやって求めるのかという話なのだが、これは割と有名な話で、Find first set(ffs)ないしFind first one(→最初の1を見つける)とか、Number of trailing zeros(NTZ)ないしCount trailing zeros(CTZ)(→末尾の0の個数)とかいう名前が付けられている。ここでは名前がわかりやすいNTZを採用する。
　\(\mathrm{NTZ}(x)\)を、\(x\) を2進表記したときの末尾に0がいくつ続くかを返す関数としよう。すると \(\mathrm{NTZ}(x)\) は \(\log(\verb|-x & x|)\) という至極単純な形で表される。ただしここで \(\log\) の底は2で、-x は x の2の補数である。そしてこの節の冒頭の段落を再度見直してもらうと、この \(\mathrm{NTZ}(x)\) というのは \(x-1\) に1を足して \(x\) にするという操作における \(k\) の値にほかならないことがわかるだろう。したがって
$$\begin{eqnarray*}2^k &=& 2^{\mathrm{NTZ}(x)} \\
&=& 2^{\log(\verb|-x & x|)} \\
&=& \verb|-x & x|
\end{eqnarray*}
$$
あるいは、NTZにビットを反転した ~x を与えることで、末尾に1がいくつ続くかを返す NTO (Number of trailing ones. 名前は適当) なる関数とすることもできる。この場合は\(x\) に1を足して \(x+1\) にするという操作における \(k\) の値が求まることになる。あえてビット演算を展開するならば、-x = ~x + 1 (2の補数)であることから
-(~x)& ~x = (~(~x) + 1) & ~x = (x + 1) & ~x であるので、\(\mathrm{NTO}(x)\) を \(\log(\verb|(x + 1) & ~x|)\) で定義することもできる。

以上をまとめると

  • NTZを使って、末尾に0が続く数 \(k\) に対し \(2^k\) を求める
  • \(2^{N-1}\)を \(2^k\) で割ったものを\(2^N\)から引き、ビットマスクを求める
  • ビットマスクを一つ前のビット逆転列とXORして、次のビット逆転列を求める

という順序で計算ができる。この方法では一つ前のビット逆転列の値さえピンポイントで分かっていれば次のビット逆転列を導けるので、本当の意味での「インクリメント」といえるかもしれない。

実装を示す。

def bit_reverse_xor(bits):
    table = [0 for _ in range(1 << bits)]
    s = 0
    m = 1 << bits
    for i in range(1 << bits):
        table[i] = s
        s ^= m - m // (2 * (~i & (i + 1)))
    return table


def bit_reverse_xor2(bits):
    table = [0 for _ in range(1 << bits)]
    s = 0
    m = 1 << bits
    for i in range(1, 1 << bits):
        s ^= m - m // (2 * (i & -i))
        table[i] = s
    return table

bit_reverse_xorがNTOを使った実装で、bit_reverse_xor2がNTZを使った実装である。計算量は\(O(M)\)である。

・パフォーマンスバトル

　それでは各実装でどれくらい差があるか見てみよう。以下の様に関数を渡して時間を測定するメソッドを定義した。

def test_performance(bits, *func):
    for f in func:
        st = time.perf_counter()
        revt = f(bits)
        en = time.perf_counter()

        print(f"{f.__name__} time: {en-st:.6f} sec.")

まずは20 bitのテーブルを生成してみる。


bit_reverse_naive time: 3.045853 sec.
bit_reverse_naive2 time: 1.935180 sec.
bit_reverse_recursive time: 0.164151 sec.
bit_reverse_xor time: 0.207048 sec.
bit_reverse_xor2 time: 0.198910 sec.

既にnaiveな実装が10 ~ 15倍ほど遅いことが判明してしまった。計算量のオーダー的には \(\log M = N\) 倍なのでまあ妥当な結果だろう。

\(O(M)\) の皆さんにはもっと大きなテーブルを生成してもらう。25 bitではどうだろうか。


bit_reverse_recursive time: 5.226015 sec.
bit_reverse_xor time: 6.664798 sec.
bit_reverse_xor2 time: 6.280106 sec.

先ほどもちょっと早かったが、やはり再帰的な実装がやや早い。XORを使った実装ではNTOよりNTZを使った実装のほうが僅かに早いようである。NTOのほうがビット演算の回数が1回多いので、その影響だろう。

・結論

　大体の用途では再帰的な実装で全部のテーブルを生成してしまうのが手っ取り早いだろう。

この記事で掲げたコードはgistに置いてあるので、弄って遊びたい方は見に行ってみてください。

ビット逆転インクリメント

投稿ナビゲーション

コメントを残す

メールアドレスが公開されることはありません。