日々drdrする人のメモ

今日もdrdr、明日もdrdr

Li Chao Treeのメモ

Li Chao (Segment) Treeについて理解をまとめた

以下の問題で使ったので、関連してきちんとまとめようと思った
smijake3.hatenablog.com

参考ページは以下
Convex hull trick and Li Chao tree - Competitive Programming Algorithms

目次

Li Chao Treeとは

f:id:smijake3:20180616042909p:plain

Convex Hull Trickとして解くアルゴリズムの1つで、直線をセグメント木で管理しながら、ある点における最小値(最大値)を計算する。
セグメント木のうちの区間 {[l, r)}では、 {[l, r)}区間で最小(最大)の値をとり得る直線 {f(x) = ax + b}が保持される。

そして、計算する座標xが {N}個あるとき、直線追加と最小値(最大値)計算が1回 {O(\log N)}で行える。

Li Chao Treeにおける最小値の計算

f:id:smijake3:20180616043112p:plain

ここでは、2つのクエリ処理、直線追加とある点xにおける最小値を計算する例について説明する。セグメント木は {2^k}の要素に対しての方が計算しやすいため、 {N = 2^k}とする。

今回、 {x_0, x_1, ..., x_{N-1} (x_0 < x_1 < ... < x_{N-1})} {N}個の座標xにおける最小値を計算するとする。(任意の区間に対しても動的セグメント木で計算できるが、今回は省略)
区間 {[l, r)}では、 {x_l \le x < x_r}の範囲で最小値をとり得る直線の式(つまり、 {f(x) = ax + b})を保持する。保持する情報は {ax + b} {a, b}である。全ての区間において、始めは {(a, b) = (0, \infty)}を持つとする。

直線追加

例えば、新たに直線 {f_A(x)}を追加するとする。
区間 {[0, N)}からトップダウンに探索を行い、その探索で訪れる各区間 {B = [l, r)}に既に存在する直線 {f_B(x)}と、現在持っている直線 {f_A(x)}とを比較し、交わる関係によっては {f_A} {f_B}を交換しながら更新を行う。

具体的に、区間 {B}に存在する関数 {f_B(x)}と現在持っている関数 {f_A(x)}について、 {m = \lfloor \frac{l + r}{2} \rfloor}とし、 {x_l, x_m, x_r}のそれぞれで2つの関数がとる値を比較し、各ケースに対応する。

値を比較した結果、以下のパターンが考えられる(赤線が {f_A(x)}青線が {f_B(x)}とする)
f:id:smijake3:20180616025041p:plain
(1)のように、 {\forall x \in [l, r). f_B(x) < f_A(x)}となる場合、 {f_A(x)}は最小値を取れないため、ここで {f_A(x)}を捨てて更新を終了する。
逆に(2)のように、 {\forall x \in [l, r). f_A(x) < f_B(x)}となる場合、 {f_B(x)}が最小値を取れないため、 {f_A(x)}区間 {B}の関数とし、 {f_B(x)}を捨てて更新を終了する。

f:id:smijake3:20180616030057p:plain
(3)もしくは(4)のように、半分の区間 {[l, m)}もしくは区間 {[m, r)} {f_A(x) < f_B(x)}となる場合、 {f_A(x)}が最小値を取れる区間((3)であれば {[l, m)}、(4)であれば {[m, r)})を更に探索する。最小値が取れない方の半分の区間は探索しない。

f:id:smijake3:20180616031114p:plain
(5)もしくは(6)のように、半分の区間以上で {f_A(x) < f_B(x)}となる場合、 {f_A(x)} {f_B(x)}を交換する。こうすることで、 {f_A(x)} {f_B(x)}が逆の関係になり、(3)もしくは(4)のパターンになる。こうすることで、半分の区間を探索するだけになる。

直線追加の実装

これらを踏まえ、直線を追加する処理をPythonで実装したものが以下である。
実際に計算する座標 {x_0, ..., x_{N-1}}の個数 {N} {2^k}の形にする必要はあるが、計算上存在しない座標 {x_N, ..., x_{2^k-1}}は計算上取り得ない大きい座標にしておけばよい。

# N: 座標の数
# N0 = 2**(N-1).bit_length()
# X[i]: i番目の座標x_i (|X| = N0にする必要あり)

# None: (a, b) = (0, ∞)を表す
data = [None]*(2*N0+1)

# 座標xにおける直線の値f(x)を計算
def f(line, x):
    p, q = line
    return p*x + q

# 区間B = [l, r)の探索
def _add_line(line, k, l, r):
    m = (l + r) // 2
    # 区間Bの直線が(a, b) = (0, ∞)の場合 -> (2)のパターン
    if data[k] is None:
        data[k] = line
        return

    lx = X[l]; mx = X[m]; rx = X[r-1]
    left = (f(line, lx) < f(data[k], lx))
    mid = (f(line, mx) < f(data[k], mx))
    right = (f(line, rx) < f(data[k], rx))

    # (2)のパターン
    if left and right:
        data[k] = line
        return
    # (1)のパターン
    if not left and not right:
        return
    if mid:
        # (5), (6)のパターン -> 区間Bの直線と交換して(3), (4)のパターンに変更
        data[k], line = line, data[k]
    if left != mid:
        # (3)のパターン: [l, m)の探索
        _add_line(line, 2*k+1, l, m)
    else:
        # (4)のパターン: [m, r)の探索
        _add_line(line, 2*k+2, m, r)

def add_line(line):
    return _add_line(line, 0, 0, N0)

最小値計算

ある座標 {x}における最小値を計算する際は、座標 {x}を含む区間 {B}が持つ各関数について {f_B(x)}を計算し、それらの値のうちの最小値を求めればよい。

例えば、 {x_3}の最小値を計算するのであれば、以下のように区間 {[3, 4), [2, 4), [0, 4), [0, 8)}が持つ関数の値 {f_B(x_3)}を計算して最小値を求める。
f:id:smijake3:20180616140044p:plain

最小値計算の実装

Pythonにおける最小値計算の実装は以下になる。

def query(k):
    k += N0-1
    x = X[k]
    s = 1e30
    while k >= 0:
        if data[k]:
            s = min(s, f(data[k], x))
        k = (k - 1) // 2
    return s

線分の最小値計算

f:id:smijake3:20180616042544p:plain

Li Chao Treeにおける直線追加処理を変更することで、線分の最小値(最大値)計算もできるようになる。具体的には、直線追加では全体の区間 {[0, N}に対して更新をかけていたが、これを区間 {[a, b)}に制限し、部分的に更新するようにすればよい。
例えば、下の例のように区間 {[3, 6)}の範囲のみに含まれる( {x_3 \le x < x_6}に存在する)線分を更新する場合、区間 {[3, 4), [4, 6)}に直線を追加する処理を行えばよい。
f:id:smijake3:20180616040350p:plain

 {O(\log N)}かかる直線追加処理を、高々 {\log N}個の部分区間へ行うため、1回の線分追加処理は {O(\log^2 N)}になる。

線分追加の実装

従来の直線追加処理を少し変更して線分追加処理に変更したPython実装は以下のようになる。
トップダウンに探索を行い、追加対象の区間 {[l, r) \subseteq [a, b)}が見つかった場合にその区間への直線追加処理を行う。

# N: 座標の数
# N0 = 2**(N-1).bit_length()
# X[i]: i番目の座標x_i (|X| = N0にする必要あり)

data = [None]*(2*N0+1)

def f(line, x):
    p, q = line
    return p*x + q

# 区間B = [l, r)の探索
def _add_line(line, a, b, k, l, r):
    if r <= a or b <= l:
        # 区間[a, b)に到達できない場合は探索を打ち切る
        return
    m = (l + r) // 2
    if not (a <= l and r <= b):
        # 区間[l, r)の部分区間として[a, b)が含まれる
        # => [l, m) と [m, r)に分けてそれぞれ探索を行う
        _add_line(line, a, b, 2*k+1, l, m)
        _add_line(line, a, b, 2*k+2, m, r)
        return
    # ここ以降は、[a, b)の区間に[l, r)が含まれている場合の処理
    # (直線追加処理と同じ)

    if data[k] is None:
        data[k] = line
        return

    lx = X[l]; mx = X[m]; rx = X[r-1]
    left = (f(line, lx) < f(data[k], lx))
    mid = (f(line, mx) < f(data[k], mx))
    right = (f(line, rx) < f(data[k], rx))

    if left and right:
        data[k] = line
        return
    if not left and not right:
        return
    if mid:
        data[k], line = line, data[k]
    if left != mid:
        _add_line(line, a, b, 2*k+1, l, m)
    else:
        _add_line(line, a, b, 2*k+2, m, r)

def add_line(line, a, b):
    return _add_line(line, a, b, 0, 0, N0)

線分追加の実装(高速化)

高速化するのであれば、トップダウン探索で追加対象の区間 {[l, r) \subseteq [a, b)}へ移動するのではなく、区間 {[l, l + 2^k) \subseteq [a, b)}となる区間を列挙して、直接直線を追加するとよい。(Pythonだとこれが結構効く)

# N: 座標の数
# N0 = 2**(N-1).bit_length()
# X[i]: i番目の座標x_i (|X| = N0にする必要あり)

data = [None]*(2*N0+1)

def f(line, x):
    p, q = line
    return p*x + q

# _add_lineは直線追加処理から変更なし
def _add_line(line, k, l, r):
    m = (l + r) // 2
    if data[k] is None:
        data[k] = line
        return

    lx = X[l]; mx = X[m]; rx = X[r-1]
    left = (f(line, lx) < f(data[k], lx))
    mid = (f(line, mx) < f(data[k], mx))
    right = (f(line, rx) < f(data[k], rx))

    if left and right:
        data[k] = line
        return
    if not left and not right:
        return
    if mid:
        data[k], line = line, data[k]
    if left != mid:
        _add_line(line, 2*k+1, l, m)
    else:
        _add_line(line, 2*k+2, m, r)

def add_line(line, a, b):
    # 追加する区間[l, l+2^k)を列挙して、区間[l, l+2^k)に直線として追加
    L = a + N0; R = b + N0
    a0 = a; b0 = b
    sz = 1
    while L < R:
        if R & 1:
            R -= 1
            b0 -= sz
            _add_line(line, R-1, b0, b0+sz)
        if L & 1:
            _add_line(line, L-1, a0, a0+sz)
            L += 1
            a0 += sz
        L >>= 1; R >>= 1
        sz <<= 1