概要
Fenwick Tree がよくわからなかったので自分の言葉で整理する。
問題意識
配列 A の index 0 から i にそれぞれ対応する value の和を計算したい。
素朴な方法に、配列全ての value を見て足し合わせる、というものがある。このとき A の長さを N とするならば、この方法の時間計算量は $O(N)$ と見て良い。
もう一つの方法に、予め配列 A の index 0 から j までの和を計算しておいたものを、新たな配列 S に格納しておく、というものがある。つまり配列 S は
$$
S[j] = \sum_{i=0}^{j} A[i]
$$
ということ。この方法ならば、 S[i] の要素にアクセスするだけの時間だけを考えればいいから、時間計算量は $O(1)$ とみなせる。
ここでさらに、配列 A の要素をアップデートすることを考える。例えば、 A[i] の値がもともと 3 とかだったものを 100 で上書きすることを考える。
上に述べた2つの方法の内最初に述べた方法だと、 A の要素の素朴な書き換えだから、アップデートにかかる計算量は $O(1)$ である。もう一方の方は、 A の一つの要素を書き換えると、 S では書き換えられた要素以降の全ての要素を書き換えなければならないので、計算量は $O(N)$ とみる。
Fenwick Tree を用いればいずれも $O(\log N)$ でできる。
Fenwick Tree の構成方法
オリジナルの配列を A とする。 Fenwick Tree は A をもとに新たに作られたある配列である。 Fenwick Tree の配列を F と呼ぶ。
これから自然数と呼ぶときは 0 も含める。
Fenwick Tree を構成するために、まず lsb(least significant bit) 関数を定義する。 lsb は自然数から自然数への関数であり、 lsb(n) は「自然数 n を 2 進数展開したときに 1 が立っている最小の桁に対応する 10 進数の数」のこと。特に lsb(0) は 0 と約束する。例えば lsb(9) は、 9 を 2 進数展開すると 1001 で、ここで 1 が立ってる最小の桁は $2^0=1$ のところだから lsb(9) = 1. ここから奇数 n のときは lsb(n) = 1 だとわかる。他の例だと、 lsb(12) は、 12 を 2 進数展開すると 1100 で、ここで1が立ってる最小の桁は $2^2=4$ だから lsb(12) = 4.
この lsb を用いて以下の関数 T を考える:
$$
T: \mathbb{N} \ni n \longmapsto n - \mathrm{lsb(n)} \in \mathbb{N}.
$$
例えば、 T(12) = 12 - lsb(12) = 12 - 4 = 8。直感的には、 2 進数展開したときの一番右側にある 1 を 0 に置き換える、というもの。 この例では 1100 を 1000 にしている。
$T^0$ を恒等関数、すなわち $T^0 (n) = n$ とする。また $T$ の m 回の iteration を $T^m = T^{m-1} \circ T$ と表現する。
任意の自然数 n に対して、 $T^m (n)$ は、 m をどんどん大きくすると 0 に収束する。以降でそれを説明する。任意の自然数 n を固定して考える。 $T^m$ は m に関して減少関数である。これは lsb(n) が非負であることから明らか。 さらに $T(n) = n$ が成立するならば $lsb(n) = 0$ であり、すなわちそれは n = 0 であるから、その対偶を考えれば $n \neq 0$ なるときは狭義の減少関数となる( $n > T(n)$). $T^m (n) \geq 0$ であることを考えると、ある M が存在して M < m なる任意の m に対して $T^m (n) = 0$ となる。
さて、 Fenwick Tree F を以下のように構成する。
F の index 0 を F[0] = 0 とする。 これに本質的な意味はないが、実装の便宜上そうしている。
$i \neq 0$ なる i について、
$$
F[i] = A[i-1] + \sum_{m=0}^{M} F[T^m (i-1)]
$$
という風に、再帰的に定義する。ここで $M := \log_2 \mathrm{lsb} (i)$.
ここで注意だが、 F の長さは A の長さよりも 1 だけ長い。なぜなら F[0] が 0 であるため。そのため、 F と A は 1 ずつ index がずれている。上の式でもケアしているように、 F[i] と A[i-1] が対応している。
この式の例を出すと、例えば F[8] を構成するには、 A[7] と(上の注意を考慮せよ)、 F[7], F[6], F[4] を全部足し合わせれば OK. これを特に F の index について 2 進数で表現すると、そのパターンを悟ることができる:
8: 1000
---
7: 0111
6: 0110
4: 0100
同様に F[24] を考えると、それは A[23] と、 F[23], F[22], F[20] を全て足し合わせたものになる。これも 2 進数で考えると、
24: 11000
---
23: 10111
22: 10110
20: 10100
というかんじ。
Fenwick Tree を用いて和を求める
$$
\sum_{i=0}^{L-1} A[i] = \sum_{l=0}^{Q} F[T^l (L)]
$$,
特に $Q := \log_2 \mathrm{msb} (L)$ とする。 $\mathrm{msb} (n)$ は、 lsb 同様自然数から自然数に対する関数で、 lsb とは逆の「自然数 n を 2 進数展開したときに 1 が立っている最大の桁に対応する 10 進数の数」である。例えば msb(12) は 1100 であることから一番左の 1 の桁 $2^3 = 8$ に注目して msb(12) = 8. msb は most significant bit の略。
実は Q は理論上は特に意味がない。 $Q = \infty$ としても同じ意味に成る。なぜなら 任意の n に対して $\lim_{m \rightarrow \infty} T^m (n) = 0$ であるから。 Q を有限の数として表したのは、とはいえ実装上では気にする必要があるから。
例えば A[0] から A[31] までの和を求めたいときは、 F[32] だけを見れば OK. 他にも例えば A[0] から A[14] までの和は、 F[15], F[14], F[12], F[8] を全部足すだけで OK. これも以下のように 2 進数展開したものを比べれば、その雰囲気がわかる。1をどんどん右から減らした数に対応する index の F の value を足してるだけ。
15: 1111
14: 1110
12: 1100
8: 1000
計算量を考えると、 Fenwick Tree では足し算する回数が $Q = \log_2 \mathrm{msb} (L)$ なので $O( \log_2 L)$ となる。
実装して検証する
# %%
import random
from typing import List
from math import log2, floor
import time
# %%
N = 10000000
a = [random.randint(0, 100) for _ in range(N)]
# %%
# 普通にやるとどれくらい時間がかかるか。
# from 0 to i の index の和を求める。
print('N', N)
start_time = time.time()
result = 0
for i in a:
result += i
print(result)
end_time = time.time()
print('duration', end_time - start_time)
結果:
N 10000000
500195758
duration 0.7198359966278076
つまり一千万の配列の全ての和を素朴に計算すると、大体 0.7 秒くらいかかる。
以下で、 Fenwick Tree を構成して同様に一千万の配列の全ての和を計算する。
# %%
def lsb(i: int) -> int:
if i <= 0:
return 0
return i & -i
def func_t(n: int) -> int:
return n - lsb(n)
def initialize_fenwrick_tree(array: List) -> List:
return_array = [0]
for index in range(len(array)):
# 格納する value
node_value = 0
# 都合により index は 1 から始める。
index = index + 1
# これから足し算をする上限 M.
m = int(log2(lsb(index)))
# print('m', m)
# まず対象の index に対応する value を格納する。
# 今 index を 1 から数え直していることに注意する
node_value += array[index - 1]
# print('node_value', node_value)
# initialize の式の通り、必要な node の value を足していく。
cur_index = index - 1
# print('cur_index', cur_index)
for _ in range(m):
node_value += return_array[cur_index]
cur_index = func_t(cur_index)
# 格納する。
return_array.append(node_value)
return return_array
def sum_f_tree(f_tree: List, n: int) -> float:
cur_n = n
result = 0
while cur_n > 0:
result += f_tree[cur_n]
cur_n = func_t(cur_n)
return result
# やってみる
s = time.time()
f_tree_a = initialize_fenwrick_tree(a)
print('time', time.time() - s)
# %%
s = time.time()
print(sum_f_tree(f_tree_a, N))
print('time', time.time() - s)
結果:
time 9.028936862945557
500195758
time 0.0001239776611328125
つまり一千万の配列の全ての和を Fenwick Tree を用いて計算すると、大体 1 マイクロ秒以下。圧倒的に速い。
なんだけど、 initialize に 9 秒かかっているので、このままだと使い物にはならない。
Fenwick Tree の更新
A[j-1] の値を書き換えるとすると、 F では何を書き換えなければならないか？
それは以下のようになる。まず自然数から自然数への関数 L を left shift 関数とする。それは「2進数展開したやつを左に一度ずらす」。すなわち L(n) = 2n. F の更新は L を用いて以下のようにする。任意の自然数 m に対して、
$$
F[j - L^m \mathrm{lsb} (j) + L^{m+1} \mathrm{lsb} (j)]
$$
も更新する。
m は理論上無限個なんだけど、 left shift するのもあるところからずっと L(n) = 0 になるので、どこかで意味がなくなる。その有限の値はどれくらいかというと、 L の Operation が意味がなくなるとこだから、配列 A の長さを N にするなえば $\log_2 N$ 回になる。つまり、計算量は $O(\log_2 N)$ になる。
以下実装。
def update_f_tree(f_tree: List, org_list: List, updated_index: int, to_value: float) -> List:
# 足し算する diff を計算する。ちなみに updated_index は f_tree を基準にしているので、
# original list にアクセスするときは updated_index - 1 にする。
diff = to_value - org_list[updated_index - 1]
# iterate する index の一覧をとる。
cur_index = updated_index
updated_indice = []
while cur_index < len(f_tree):
updated_indice.append(cur_index)
cur_index = cur_index - lsb(cur_index) + (lsb(cur_index) << 1)
print(updated_indice)
for i in updated_indice:
f_tree[i] += diff
return f_tree
# %%
modified_f_tree = update_f_tree(f_tree, input_l, 9, 2)
modified_f_tree
# %%
sum_f_tree(modified_f_tree, 10)
# %%
s = time.time()
modified_f_tree_a = update_f_tree(f_tree_a, a, 500000, 900000000000000000)
print('duration', time.time() - s)
sum_f_tree(modified_f_tree_a, N)
結果:
[500000, 500032, 500096, 500224, 500736, 501760, 503808, 507904, 524288, 1048576, 2097152, 4194304, 8388608]
duration 0.00012683868408203125
900000180500195500