Fortranでinsert/eraseが早いvector作成記

　はじめに

　arrayListは配列より柔軟で実装も軽い便利なデータ構造だが、「末尾以外の要素に対する挿入/削除」の計算量が$O(N)$であり、ボトルネックである。
　それに対してlinkedListでは、このボトルネックは解消されており、末尾以外の要素に対する挿入/削除を$O(1)$で行うことができるが、arrayListで$O(1)$だった「i番目の値の取得」が$O(N)$と悪化している。

　主要なデータ構造に対して、基本的な操作を行った際の計算量を以下に示す。

	arrayList	linkedList	deque(arrayList)
i番目の値を検索する	$O(1)$	$O(N)$	$O(1)$
末尾の挿入/削除	$O(1)$	$O(1)$	$O(1)$
先頭への挿入/削除	$O(N)$	$O(1)$	$O(1)$
i番目の値の挿入/削除	$O(N)$	$O(1)$	$O(N)$

　この表から、どのデータ構造も、表に示した4つの操作の少なくとも1つの計算量が$O(N)$だとわかる。
　一方、今回紹介するデータ構造(vectree)に対する同様な表を以下に示す。

	vectree
i番目の値を検索する	$O(\log(N))$
末尾の挿入/削除	$O(\log(N))$
先頭への挿入/削除	$O(\log(N))$
i番目の値の挿入/削除	$O(\log(N))$

　この表から分かるように、vectreeは表に示した4つの操作に対して$O(\log(N))$で動作する。vectreeは「arrayList,linkedList,dequeのいいとこ取りをして、ちょっとデバフをかけた」ようなデータ構造だといえる。

　競技プログラミング的事情から見ても、$O(N\log N)$の処理は、おそらく$N = 10^5$~$10^6$くらいまで耐えられるため、実用的であろうと考えられる。

vectreeとは

　まず、vectreeという名前だが、これは正式なものではなく、この記事でつけられた名称であり、今回はこの名称を継承している。ここを参考にしたらしい。

　上の記事に図などを用いた丁寧な解説があるため、そちらを参考にしたほうが良いと思うが、vectreeとは、
* 平衡二分木によってデータを保持している。
* i番目の要素は、木を通りがけ順に見ていったときのi番目のノードに保存
* i番目の要素から見たi-1番目の要素のことをpredecessor(前任者)と呼ぶ
* i番目の要素から見たi+1番目の要素のことをsuccesor(後継者)と呼ぶ
* i番目への挿入では、predecessorの右の子ノードに挿入する
* i番目の要素の削除は、predecessorまたはsuccesorをi番目の位置へ移動させる
* i番目の要素を検索するために、各ノードが「自身を根とする部分木の要素数」を保持する
* 「自身を根とする部分木の要素数」は、挿入/削除の度に更新する。
　というものらしい。今回はこのvectreeをFortranで実装することにした。

実装

　平衡二分木にはAA木を用いた。リンク先のコードがFortranに似ていたため、簡易に実装することができた。
　実装は、
* 平衡二分木の構造をvectree_node_mod
* vectree_node_modのラッパーをvectree_mod
　というようにした。平衡二分木の実装は、vectreeの挿入の操作が普通の平衡二分木と異なるため、AA木の実装単体を分離することはできなかった。
　また、procedure名はc++に倣った。

環境

GNU Fortran (Ubuntu 7.5.0-3ubuntu1~18.04) 7.5.0
AtCoderのコードテスト

vectree_node_mod

　各procedureにはvtn_という接頭語をつけた。vectreeを使う時に、このmoduleの内部を気にする必要はない。

module vec_tree_node_mod
    use,intrinsic :: iso_fortran_env
    implicit none
    private
    public:: vtn_insert, vtn_push_back
    public:: vtn_erase
    public:: vtn_to_array
    public:: vtn_print, vtn_print_debug
    public:: vtn_at
    type,public:: vec_tree_node
        integer(int32):: level=1,tree_size=1
        integer(int32):: value
        type(vec_tree_node),pointer:: left_child => null()
        type(vec_tree_node),pointer:: right_child => null()
    end type
contains
    function vtn_has_left_child(n) result(ret)
        type(vec_tree_node),pointer:: n
        logical:: ret

        ret = associated(n%left_child)
    end function


    function vtn_has_right_child(n) result(ret)
        type(vec_tree_node),pointer:: n
        logical:: ret

        ret = associated(n%right_child)
    end function


    function vtn_is_leaf(n) result(ret)
        type(vec_tree_node),pointer:: n
        logical:: ret

        ret = .not. (vtn_has_left_child(n) .or. vtn_has_right_child(n))
    end function


    function vtn_predecessor(n) result(ret)
        !left -> rightend
        type(vec_tree_node),pointer:: n,ret

        ret => n%left_child
        do while(vtn_has_right_child(ret))
            ret => ret%right_child
        end do
    end function


    function vtn_succesor(n) result(ret)
        !right -> leftend
        type(vec_tree_node),pointer:: n,ret

        ret => n%right_child
        do while(vtn_has_left_child(ret))
            ret => ret%left_child
        end do
    end function


    subroutine vtn_size_update(n)
        type(vec_tree_node),pointer:: n
        integer(int32):: old_size, new_size

        old_size = n%tree_size
        new_size = vtn_left_size_of(n)+vtn_right_size_of(n) + 1
        n%tree_size = new_size
    end subroutine


    function vtn_left_size_of(n) result(ret)
        type(vec_tree_node),pointer:: n
        integer(int32):: ret

        ret = 0
        if (vtn_has_left_child(n)) ret = n%left_child%tree_size
    end function


    function vtn_right_size_of(n) result(ret)
        type(vec_tree_node),pointer:: n
        integer(int32):: ret

        ret = 0
        if (vtn_has_right_child(n)) ret = n%right_child%tree_size
    end function


    function vtn_ind_at(n) result(ret)
        type(vec_tree_node),pointer:: n
        integer(int32):: ret

        ret = vtn_left_size_of(n)+1
    end function


    function vtn_skew(n) result(ret)
        type(vec_tree_node),pointer:: n, ret, sub_n

        if (.not. associated(n)) then
            allocate(ret)
            ret => null()
            return
        else if (.not. vtn_has_left_child(n)) then
            ret => n
        else if (n%level == n%left_child%level) then
            sub_n => n%left_child
            n%left_child => sub_n%right_child
            sub_n%right_child => n
            ret => sub_n
        else
            ret => n
        end if
        call vtn_size_update(n)
        call vtn_size_update(ret)
    end function


    function vtn_split(n) result(ret)
        type(vec_tree_node),pointer:: n,ret,sub_n

        if (.not. associated(n)) then
            allocate(ret)
            ret => null()
            return
        else if (.not. vtn_has_right_child(n) .or. .not. vtn_has_right_child(n%right_child)) then
            ret => n
        else if (n%level == n%right_child%right_child%level) then
            sub_n => n%right_child
            n%right_child => sub_n%left_child
            sub_n%left_child => n
            sub_n%level = sub_n%level+1
            ret => sub_n
        else
            ret => n
        end if
        call vtn_size_update(n)
        call vtn_size_update(ret)
    end function


    recursive function vtn_push_back(n,i,val) result(ret)
        type(vec_tree_node),pointer:: n,ret
        integer(int32):: i
        integer(int32):: val

        if (.not. associated(n)) then
            allocate(ret)
            ret = vec_tree_node(value=val)
            return
        else
            n%right_child => vtn_push_back(n%right_child,i,val)
        end if
        n => vtn_skew(n)
        n => vtn_split(n)
        ret => n
    end function


    function vtn_level_of(n) result(ret)
        type(vec_tree_node),pointer:: n
        integer(int32):: ret

        ret = 0
        if (associated(n)) ret = n%level
    end function


    function vtn_decrease_level(n) result(ret)
        type(vec_tree_node),pointer:: n, ret
        integer(int32):: should_be

        should_be = min(vtn_level_of(n%left_child), vtn_level_of(n%right_child)) + 1
        if (should_be < n%level) then
            n%level = should_be
            if (should_be < vtn_level_of(n%right_child)) then
                n%right_child%level = should_be
            end if
        end if
        ret => n
    end function


    recursive function vtn_insert(n,i,val) result(ret)
        type(vec_tree_node),pointer:: n,ret
        integer(int32):: i
        integer(int32):: val

        if (.not. associated(n)) then
            allocate(ret)
            ret = vec_tree_node(value=val,level=1)
            return
        end if
        if (i == vtn_ind_at(n)) then
            n%left_child => vtn_push_back(n%left_child,i,val)
        else if (i > vtn_ind_at(n)) then
            n%right_child => vtn_insert(n%right_child, i-(vtn_ind_at(n)),val)
        else if (i < vtn_ind_at(n)) then
            n%left_child => vtn_insert(n%left_child,i,val)
        end if
        n => vtn_skew(n)
        n => vtn_split(n)
        ret => n
    end function


    recursive function vtn_erase(n,i,ev) result(ret)
        type(vec_tree_node),pointer:: n,ret,nn
        integer(int32):: i
        integer(int32),optional:: ev

        if (.not. associated(n)) then
            allocate(ret)
            ret => null()
            return
        else if (vtn_ind_at(n) > i) then
            n%left_child => vtn_erase(n%left_child, i,ev)
        else if (vtn_ind_at(n) < i) then
            n%right_child => vtn_erase(n%right_child, i-vtn_ind_at(n),ev)
        else
            if (present(ev)) ev = n%value
            if (vtn_is_leaf(n)) then
                allocate(ret)
                ret => null()
                return
            else if (vtn_has_left_child(n)) then
                nn => vtn_predecessor(n)
                n%left_child => vtn_erase(n%left_child,nn%tree_size)
                n%value = nn%value
            else
                nn => vtn_succesor(n)
                n%right_child => vtn_erase(n%right_child,nn%tree_size)
                n%value = nn%value
            end if
        end if
        n => vtn_decrease_level(n)
        n => vtn_skew(n)
        n%right_child => vtn_skew(n%right_child)
        if (vtn_has_left_child(n)) n%right_child%right_child => vtn_skew(n%right_child%right_child)
        n => vtn_split(n)
        n%right_child => vtn_split(n%right_child)
        ret => n
    end function


    recursive subroutine vtn_print_debug(n,i,d)
        type(vec_tree_node),pointer:: n
        integer(int32):: i,j
        character(1):: d

        if (.not. associated(n)) return
        call vtn_print(n%right_child,i+1,"r")
        do j=1,i; write(*,'(a)',advance="no") '  '; end do
        if (d == "r") then
            write(*,'(a)',advance="no") "┌ "
        else if (d == "l") then
            write(*,'(a)',advance="no") "└ "
        else
            write(*,'(a)',advance="no") "  "
        end if
        print'("[ ",3(i0,1x),"]")', n%value, n%level, n%tree_size
        call vtn_print(n%left_child,i+1,"l")
    end subroutine


    recursive subroutine vtn_print(n,i,d)
        type(vec_tree_node),pointer:: n
        integer(int32):: i,j
        character(1):: d

        if (.not. associated(n)) return
        call vtn_print(n%right_child,i+1,"r")
        do j=1,i; write(*,'(a)',advance="no") '  '; end do
        if (d == "r") then
            write(*,'(a)',advance="no") "┌ "
        else if (d == "l") then
            write(*,'(a)',advance="no") "└ "
        else
            write(*,'(a)',advance="no") "  "
        end if
        print'(i0)', n%value
        call vtn_print(n%left_child,i+1,"l")
    end subroutine


    recursive subroutine vtn_to_array(n,ind,ar)
        type(vec_tree_node),pointer:: n
        integer(int32),intent(inout):: ind
        integer(int32),intent(inout):: ar(:)
        integer(int32):: i

        if (.not. associated(n)) return
        i=ind
        call vtn_to_array(n%left_child,i,ar)
        ar(i) = n%value
        i=i+1
        call vtn_to_array(n%right_child,i,ar)
        ind=i
    end subroutine


    recursive function vtn_at(n,ind) result(ret)
        type(vec_tree_node),pointer:: n
        integer(int32):: ind
        integer(int32):: ret

        ret = 0
        if (.not. associated(n)) return
        if (vtn_ind_at(n) > ind) then
            ret = vtn_at(n%left_child, ind)
        else if (vtn_ind_at(n) < ind) then
            ret = vtn_at(n%right_child, ind-vtn_ind_at(n))
        else ! vtn_ind_at(n) == ind
            ret = n%value
        end if
    end function
end module

vectree_mod

各procedureにはvt_という接頭語をつけた。vectreeを用いるときに使用する関数とその使用方法についてはコードの後に示す。

module vec_tree_mod
    use,intrinsic :: iso_fortran_env
    use vec_tree_node_mod
    implicit none
    private
    public:: vt_print, vt_print_debug
    public:: vt_size
    public:: vt_to_array
    type,public:: vectree
        type(vec_tree_node),pointer,private:: head => null()
        integer(int32),private:: len = 0
    contains
        procedure,public:: insert => vt_insert, erase => vt_erase
        procedure,public:: push_back => vt_push_back, push_front => vt_push_front
        procedure,public:: pop => vt_pop, pop_back => vt_pop_back, pop_front => vt_pop_front
        procedure,public:: at => vt_value_at
    end type
contains
    subroutine vt_insert(vt,i,val)
        class(vectree):: vt
        type(vec_tree_node),pointer:: new
        integer(int32):: i
        integer(int32):: val

        if (i < 1 .or. vt%len+1 < i) then
            print'(a)', "inserting index is too big"
            print'(a)', "index, vectree length"
            print'(*(i0,1x))', i, vt%len
            stop
        else if (.not. associated(vt%head)) then
            allocate(new)
            new = vec_tree_node(value=val)
            vt%head => new
        else if (i == vt%len+1) then
            vt%head => vtn_push_back(vt%head,i,val)
        else
            vt%head => vtn_insert(vt%head,i,val)
        end if
        vt%len=vt%len+1
    end subroutine


    subroutine vt_push_back(vt,val)
        class(vectree):: vt
        integer(int32):: val

        call vt_insert(vt,vt%len+1,val)
    end subroutine


    subroutine vt_push_front(vt,val)
        class(vectree):: vt
        integer(int32):: val

        call vt_insert(vt,1,val)
    end subroutine


    function vt_pop(vt,i) result(erased_value)
        class(vectree):: vt
        integer(int32):: i,erased_value

        if (i < 1 .or. vt%len < i) then
            print'(a)', "remove index is too big"
            print'(a)', "index, vectree length"
            print'(*(i0,1x))', i, vt%len
            stop
        end if
        vt%head => vtn_erase(vt%head,i,erased_value)
        vt%len=vt%len-1
    end function


    subroutine vt_erase(vt,i)
        class(vectree):: vt
        integer(int32):: i,tmp

        tmp = vt%pop(i)
    end subroutine


    function vt_pop_back(vt) result(ret)
        class(vectree):: vt
        integer(int32):: ret

        ret = vt%pop(vt%len)
    end function


    function vt_pop_front(vt) result(ret)
        class(vectree):: vt
        integer(int32):: ret

        ret = vt%pop(1)
    end function


    function vt_value_at(vt,i) result(ret)
        class(vectree):: vt
        integer(int32):: i
        integer(int32):: ret

        ret = vtn_at(vt%head,i)
    end function


    subroutine vt_print(vt)
        type(vectree):: vt

        call vtn_print(vt%head,0,' ')
    end subroutine


    subroutine vt_print_debug(vt)
        type(vectree):: vt

        print*, "==[value, vtn_level_of, tree_size]=="
        call vtn_print_debug(vt%head,0,' ')
    end subroutine


    function vt_to_array(vt) result(ret)
        type(vectree):: vt
        integer(int32):: ret(vt%len),i
        i=1
        call vtn_to_array(vt%head,i,ret)
    end function


    function vt_size(vt) result(ret)
        type(vectree):: vt
        integer(int32):: ret

        ret = vt%len
    end function
end module

実装したprocedureとその使用方法

　以下のコードで、vectreeを操作するためのすべてのprocedureについて紹介する。10個作成した。作りすぎた気がする。

program main
    use,intrinsic :: iso_fortran_env
    use vec_tree_mod
    implicit none
    type(vectree):: vt
    integer(int32):: i,n
    integer(int32),allocatable:: ar(:)

    n=5
    print"(a)", "===insert, push==="
    do i=1,2*n
        v = i*2
        call vt%insert(i,v) ! i番目にvを挿入 O(log N)
        v = i*3
        call vt%push_front(v) ! 先頭にvを挿入 O(log N)
        v = i*4
        call vt%push_back(v) ! 末尾にvを挿入 O(log N)
    end do
    print"(a)", "===vt_print==="
    call vt_print(vt) ! 平衡二分木を可愛らしく表示する O(N)
    print"(a)", "===erase, pop==="
    do i=1,n/2
        call vt%erase(i) ! i番目の要素を削除 O(log N)
        print*, vt%pop(i) ! i番目の要素を削除し、その値が返る O(log N)
        print*, vt%pop_front() ! 先頭の要素を削除し、その値が返る O(log N)
        print*, vt%pop_back() ! 末尾の要素を削除し、その値が返る O(log N)
    end do

    do i=1,n/2
        print*, vt%at(i) ! i番目の要素を参照 O(log N)
    end do

    n = vt_size(vt) ! vtの要素数を返す O(1)
    allocate(ar(n))
    ar(1:n) = vt_to_array(vt) ! vtの要素全体を配列で返す O(N)
    print"(a)", "===to_array==="
    print'(*(i0,1x))', ar(:)
end program main

実行結果

===insert, push===
===vt_print===
            ┌ 40
          ┌ 36
        ┌ 32
          └ 28
      ┌ 24
        └ 20
    ┌ 16
        ┌ 12
      └ 8
        └ 4
  ┌ 2
        ┌ 4
      ┌ 6
        └ 8
    └ 10
          ┌ 12
        ┌ 14
          └ 16
      └ 18
        └ 20
  3
      ┌ 6
    ┌ 9
      └ 12
  └ 15
        ┌ 18
      ┌ 21
        └ 24
    └ 27
      └ 30
===erase, pop===
          27
          24
          40
          15
          21
          36
          12
           9
===to_array===
12 9 6 3 20 18 16 14 12 10 8 6 4 2 4 8 12 16 20 24 28 32

バグなど

ランダムなinsert/eraseのクエリ($10^5$程度)をpythonで作成し、その出力をpythonのlist,fortranの自作vector、vectreeの3つでdiffを取って問題ない事を確認済み。また、AtCoderの「コードテスト」上でも正常に動くことを確認。(自作vectorのinsert/eraseの速度がpythonのリストに負けていて泣きそうになった。)

もしバグをみつけたら教えてください。

まとめ

Fortranでvectreeを作成したよ。