今日から始めるWingBIT - 非再帰SegmentTreeのようなもの

再帰SegmentTreeと一緒にはしたくない()

経緯

トイレでなんとなく思いついたデータ構造がうまくいったので好きになっているだけの話です(論文があったというのは内緒で)

BIT(BinaryIndexedTree)を知っていますか?

このような感じでそれぞれの要素を持ち,[0,k]の和を求めるデータ構造です.

f:id:Kutimoti:20181005223207p:plain

O(logN)区間和が取得できる上,SegmentTreeより定数倍が早い,実装が軽いことで有名ですね.

しかし,,,

BITは和以外では使いみちがあまりわからない

RMQ(Range Minimum Query)などは処理することができません...悲しいね

そこで†WingBIT†なるものを紹介します.

WingBIT

f:id:Kutimoti:20181005223224p:plain

このように今まで空いていたスペースにBITをもう一枚入れる感じのデータ構造です.(名前の由来はここから来ています)

いまから具体的にクエリをどのように捌くか説明します.

点更新 - update

index = 3の場所を更新する例です.

だんだん高さをあげていき,updateをします.

まず,3の真上にあるBlueとRedの場所をマークします.

f:id:Kutimoti:20181005223240p:plain Blue(3)は一番下のノードなので,値を更新して,次のBlueにマークを移動させます.

移動は+depthをするとできます.

f:id:Kutimoti:20181005223254p:plain

次にdepth=2をみるとRed(6)があります. これを下の2つのノードを使って更新します. そして次のRedにマークを移動させます.

f:id:Kutimoti:20181005223308p:plain

depth=4をみるとBlue(4)があります 同じように下の2つのノードを使って更新します. そして次のBlueにマークを移動させます.

f:id:Kutimoti:20181005223319p:plain

これを繰り返し,depthが一番上に行くまで繰り返すことで,更新が出来ます.

区間取得 - get_inter

一例を見てみましょう.

f:id:Kutimoti:20181005223333p:plain

f:id:Kutimoti:20181005223344p:plain

区間は必ず,左側がRed,右側がBlueでできていませんか?これを使って実装していきます.(RedとBlueでWingが出来ています.)

先と同じように下から処理していきます.

左側をRedにマーク,右側をBlueにマークします.

f:id:Kutimoti:20181005223359p:plain

depth=1を見ると,Blue(7)があります.これは区間の一部なので,利用していきます.

次の位置は-depthをすることで移動できます.

f:id:Kutimoti:20181005223404p:plain

depth=2を見ると,Red(6),Blue(6)があります.これも区間の一部なので,利用します.

f:id:Kutimoti:20181005223407p:plain

RedとBlueが交差したので終了です.

f:id:Kutimoti:20181005223410p:plain

実装例

Nim Langで書いています.

shl...左シフト
shr...右シフト
seq...C++のvectorと思っていただければ
import algorithm
type
    WingBIT[Monoid] = object
        rw : seq[Monoid]
        lw : seq[Monoid]
        ide : Monoid
        sz : int
        update_func : proc(node : Monoid , x : Monoid) : Monoid
        f : proc(x : Monoid , y : Monoid) : Monoid
proc newWingBIT*[Monoid](n : int , ide : Monoid,
    update_func : proc(node : Monoid , x : Monoid) : Monoid , 
    f : proc(x : Monoid , y : Monoid) : Monoid) : WingBIT[Monoid] =
    var bit : WingBIT[Monoid]
    bit.rw = @[]
    bit.lw = @[]
    bit.sz = 1
    while bit.sz < n: bit.sz = bit.sz * 2
    bit.rw.setLen(bit.sz + 1)
    bit.lw.setLen(bit.sz + 1)
    bit.rw.fill(ide)
    bit.lw.fill(ide)
    bit.ide = ide
    bit.update_func = update_func
    bit.f = f
    return bit

proc update*[Monoid](bit : var WingBIT[Monoid] , k : int , x : Monoid) =
    var depth = 1
    var right = k
    var left = bit.sz + 1 - k
    if (right and depth) > 0:
        bit.rw[right] = bit.update_func(bit.rw[right],x)
        right = right + depth
    if (left and depth) > 0:
        bit.lw[left] = bit.update_func(bit.lw[left],x)
        left = left + depth
    depth = 2
    while depth <= bit.sz:
        var dd = depth shr 1
        if (left and depth) > 0:
            bit.lw[left] = bit.f(bit.rw[bit.sz - left + dd],bit.lw[left - dd])
            left = left + depth
        if (right and depth) > 0:
            bit.rw[right] = bit.f(bit.rw[right - dd],bit.lw[bit.sz - right + dd])
            right = right + depth
        depth = depth shl 1

proc get_inter*[Monoid](bit : WingBIT[Monoid] , left : int , right : int) : Monoid =
    var al = bit.ide
    var ar = bit.ide
    var depth = 1
    var l = bit.sz + 1 - left
    var r = right
    while bit.sz + 1 - l <= r:
        if (l != bit.sz) and ((l and depth) > 0):
            al = bit.f(al,bit.lw[l])
            l -= depth
        if (r and depth) > 0:
            ar = bit.f(bit.rw[r],ar)
            r -= depth
        depth = depth shl 1
    return bit.f(ar,al)

# verify arc008 - d

import strutils
import sequtils

var temp = stdin.readLine.split.map(parseInt)

var N = temp[0]
var M = temp[1]

var p : array[101010,int64]
var a : array[101010,float64]
var b : array[101010,float64]

var se : seq[int64] = @[]

for i in 0..<M:
    var temp2 = stdin.readLine.split
    p[i] = temp2[0].parseInt
    a[i] = temp2[1].parseFloat
    b[i] = temp2[2].parseFloat
    se.add(p[i])
sort(se,system.cmp)

type TT = tuple[a : float64 , b : float64]
proc ffunc(x : TT,y : TT) : TT=
    return (y.a * x.a , y.a * x.b + y.b)

proc uupdate(x : TT,y : TT) : TT=
    return y

var ran = newWingBIT(M + 1000,(1.0,0.0),uupdate,ffunc)

var mr : float64 = 1.0
var ir : float64 = 1.0
for i in 0..<M:
    p[i] = lowerBound(se,p[i]) + 1
    ran.update((int)p[i],(a[i],b[i]))
    var tup : TT = ran.rw[ran.sz]
    var rrr : float64 = tup.a + tup.b
    mr = max(mr,rrr)
    ir = min(ir,rrr)
echo ir.formatFloat
echo mr.formatFloat

利点

再帰SegmentTreeより早い.(無駄が無い)

実装がそんなに大変なわけでもない

名前がかっこいい