Fortran で skip list を作ってみた

動機

以前 Fortran で平衡二分探索木の一つである AA 木を実装したのですが、比較的実装が軽い AA 木でも実装は割と大変で、アルゴリズムも分かりづらいものでした。そこで、もっと実装が楽でアルゴリズムも分かりやすい同様のデータ構造はないかな、と平衡二分探索木の Wikipedia を見ていたところ、

木構造ではないが、同じような用途に使えるものとしてスキップリストがある。

という記述を見つけました。スキップリストについて調べてみると、アルゴリズムが分かりやすく面白かったので実装することにしました。

skip list とは

skip list は平衡二分探索木に似た用途に使われる乱択アルゴリズムのデータ構造です。構造自体は連結リスト（linked list）を拡張したものとなっています。

（William Pugh (1990), A Skip List Cookbook, Figure 1. より引用）

上の図のように、skip list は複数の linked list が層状に重なった構造となっており、各ノードの縦方向の高さは確率的に決められます。ノードを挿入する際はまず最下層に一つノードを置き、例えば「コインを振って裏が出るまで」上方向にノードを伸ばし続けます。つまり、最下層はソートされた全ての要素を含む linked list となっており、その上は最下層の要素を50%、25%、12.5%…程度含む linked list となります。

skip list におけるデータ探索アルゴリズムは、電車を想像すると分かりやすいです。上層ほど「特急」、下層ほど「鈍行」のような役割を担っています。目的のノードまでは出来るだけ「特急」を使って向かいます。例えば上の図で「21」のキーに格納されている値を取得したい場合は、左上をスタート位置として「スタート位置」→「6」→「下」→「下」→「9」→「17」→「下」→「19」→「21」のように遷移します。

また、それぞれのノードに次のノードまでの距離を保持することで指定した位置のノードへのランダムアクセスが比較的高速に行えます。詳しくは下記の参考を参照してください。

計算量の比較

一般的な平衡二分探索木と skip list で、要素数を $N$ としたときの計算量は以下の通りです。

操作	平衡二分探索木	skip list
要素の挿入	$O(\log N)$	$O(\log N)$
要素の探索	$O(\log N)$	$O(\log N)$
要素の削除	$O(\log N)$	$O(\log N)$
$k$ 番目の要素の取得	$O(N)$	$O(\log N)$
キーの位置の取得	$O(N)$	$O(\log N)$

環境

macOS Mojave (10.14.6)
GNU Fortran (GCC) 8.2.0

プログラム

実際に実装したソースコードは以下のとおりです。

skip_list_map.f08

module mod_skip_list_map
  implicit none
  integer, private, parameter :: key_kind = 4 ! key の kind
  integer, private, parameter :: val_kind = 4 ! val の kind
  real(8), private, parameter :: threshold = 0.5d0 ! ノードを上に伸ばすときの基準
  integer(key_kind), private, parameter :: infty = lshift(1_key_kind, 8 * key_kind - 2) ! infinity

  ! ノード
  type t_node
    private
    integer(key_kind) :: key
    integer(val_kind) :: val
    integer :: length = 0
    type(t_node), pointer :: prev => null(), next => null()
    type(t_node), pointer :: above => null(), below => null()
  end type

  ! skip list
  type t_skip_list_map
    integer :: size = 0
    integer :: level = 0
    integer(val_kind) :: default = -1
    type(t_node), private, pointer :: head => null()
    type(t_node), private, pointer :: tail => null()
  contains
    procedure :: search => search
    procedure :: contains => contain
    procedure :: insert => insert
    procedure :: index_of => index_of
    procedure :: get_key_at => get_key_at
    procedure :: remove => remove
    procedure :: remove_key_at => remove_key_at
    procedure :: first_key => first_key
    procedure :: poll_first => poll_first
    procedure :: last_key => last_key
    procedure :: poll_last => poll_last
    procedure :: floor_key => floor_key
    procedure :: lower_key => lower_key
    procedure :: ceiling_key => ceiling_key
    procedure :: higher_key => higher_key
    final :: finalize
  end type

  ! skip list のコンストラクタ
  interface skip_list_map
    module procedure :: newslm0, newslm1
  end interface

  private :: t_node, random_seed_clock, new_node, less, finalize, increase_level, search_node
contains

  ! 乱数のシード値を変えるサブルーチン
  subroutine random_seed_clock()
    integer :: nseed, clock
    integer, allocatable :: seed(:)
    call system_clock(clock)
    call random_seed(size = nseed)
    allocate(seed(nseed))
    seed = clock
    call random_seed(put = seed)
    deallocate(seed)
  end

  ! コインを振る
  logical function flip_coin() result(res)
    real(8) :: rand
    call random_number(rand)
    res = rand < threshold
  end

  ! ノードのコンストラクタ
  function new_node(key, val, length) result(res)
    integer(key_kind), intent(in) :: key
    integer(val_kind), intent(in) :: val
    integer, intent(in) :: length
    type(t_node), pointer :: res
    allocate(res)
    res%key = key
    res%val = val
    res%length = length
  end

  ! キーの comparator
  logical function less(key1, key2) result(res)
    integer(key_kind), intent(in) :: key1, key2
    res = key1 < key2
  end

  ! skip list のデストラクタ
  subroutine finalize(this)
    type(t_skip_list_map), intent(inout) :: this
    call clear(this)
  end

  ! skip list の全要素を削除する
  subroutine clear(this)
    class(t_skip_list_map), intent(inout) :: this
    type(t_node), pointer :: node, next, above
    if (.not.associated(this%head)) return
    node => this%head
    do while (associated(node%below))
      node => node%below
    end do
    do while (associated(node))
      next => node%next
      do while (associated(node))
        above => node%above
        deallocate(node)
        node => above
      end do
      node => next
    end do
    this%head => null()
    this%tail => null()
    this%size = 0
  end

  ! skip list のコンストラクタ
  type(t_skip_list_map) function newslm0() result(res)
    type(t_node), pointer :: head, tail
    call random_seed_clock()
    head => new_node(-infty, 0, 1)
    tail => new_node(infty, 0, infty)
    head%next => tail
    tail%next => head
    res%head => head
    res%tail => tail
  end

  ! skip list のコンストラクタ
  type(t_skip_list_map) function newslm1(default) result(res)
    integer(val_kind), intent(in) :: default
    res = newslm0()
    res%default = default
  end

  ! skip list の階層数を増やすサブルーチン
  subroutine increase_level(this, level)
    type(t_skip_list_map), intent(inout) :: this
    integer, intent(in) :: level
    type(t_node), pointer :: head, tail, habove, tabove
    integer :: i
    if (this%level >= level) return
    head => this%head
    tail => this%tail
    do i = 1, this%level
      head => head%above
      tail => tail%above
    end do
    do i = this%level + 1, level
      habove => new_node(-infty, 0, 1)
      head%above => habove
      habove%below => head
      tabove => new_node(infty, 0, infty)
      tail%above => tabove
      tabove%below => tail
      head => habove
      tail => tabove
      head%next => tail
      tail%prev => head
    end do
    this%level = level
  end

  ! node%key <= key となるノードの探索
  function search_node(this, key) result(res)
    type(t_skip_list_map), intent(in) :: this
    integer(key_kind), intent(in) :: key
    type(t_node), pointer :: res
    res => this%head
    do while (associated(res%above))
      res => res%above
    end do
    do
      do while (.not.less(key, res%next%key))
        res => res%next
      end do
      if (.not.associated(res%below)) exit
      res => res%below
    end do
  end

  ! node%key <= key となるノードの val を取得
  integer(val_kind) function search(this, key) result(res)
    class(t_skip_list_map), intent(in) :: this
    integer(key_kind), intent(in) :: key
    type(t_node), pointer :: node
    node => search_node(this, key)
    res = merge(node%val, this%default, node%key == key)
  end

  ! node%key == key となるノードが存在するかどうか
  logical function contain(this, key) result(res)
    class(t_skip_list_map), intent(in) :: this
    integer(key_kind), intent(in) :: key
    type(t_node), pointer :: node
    node => search_node(this, key)
    res = node%key == key
  end

  ! node%key == key となるノードが存在する場合は val を書き換え、
  ! 存在しない場合は新しいノードを挿入する
  subroutine insert(this, key, val)
    class(t_skip_list_map), intent(inout) :: this
    integer(key_kind), intent(in) :: key
    integer(val_kind), intent(in) :: val
    type(t_node), pointer :: node, prev, next, above
    integer :: i, level, length, prevlength
    prev => search_node(this, key)
    if (prev%key == key) then
      prev%val = val
      return
    end if
    this%size = this%size + 1
    node => new_node(key, val, 1)
    next => prev%next
    prev%next => node
    node%prev => prev
    node%next => next
    next%prev => node
    level = 0
    do while (flip_coin()) ! 裏が出るまでコインを振り続ける
      level = level + 1
    end do
    call increase_level(this, level)
    prevlength = 1
    length = 1
    do i = 1, level
      do while (.not.associated(prev%above))
        prev => prev%prev
        prevlength = prevlength + prev%length
      end do
      prev => prev%above
      prev%length = prevlength
      do while (.not.associated(next%above))
        length = length + next%length
        next => next%next
      end do
      next => next%above
      above => new_node(key, val, length)
      above%below => node
      node%above => above
      node => above
      prev%next => node
      node%prev => prev
      node%next => next
      next%prev => node
    end do
    do i = level + 1, this%level
      do while (.not.associated(prev%above))
        prev => prev%prev
      end do
      prev => prev%above
      prev%length = prev%length + 1
    end do
  end

  ! node%key == key となるノードが存在する場合は先頭から何番目か、
  ! 存在しない場合はそのキーを挿入するとき何番目になるかを返す
  integer function index_of(this, key) result(res)
    class(t_skip_list_map), intent(in) :: this
    integer(key_kind), intent(in) :: key
    type(t_node), pointer :: node
    res = 0
    node => this%head
    do while (associated(node%above))
      node => node%above
    end do
    do
      do while (.not.less(key, node%next%key))
        res = res + node%length
        node => node%next
      end do
      if (node%key == key) exit
      if (.not.associated(node%below)) exit
      node => node%below
    end do
    if (node%key /= key) res = -(res + 1)
  end

  ! idx 番目のノードのキーを取得
  integer(key_kind) function get_key_at(this, idx) result(res)
    class(t_skip_list_map), intent(in) :: this
    integer, intent(in) :: idx
    integer :: length
    type(t_node), pointer :: node
    if (idx < 1) then
      res = -infty
      return
    end if
    if (idx > this%size) then
      res = infty
      return
    end if
    length = 0
    node => this%head
    do while (associated(node%above))
      node => node%above
    end do
    do
      do while (length + node%length <= idx)
        length = length + node%length
        node => node%next
      end do
      if (length == idx) exit
      if (.not.associated(node%below)) exit
      node => node%below
    end do
    res = node%key
  end

  ! node%key == key となるノードが存在する場合はそのノードを削除、
  ! 存在しない場合は何もしない
  subroutine remove(this, key)
    class(t_skip_list_map), intent(inout) :: this
    integer(key_kind), intent(in) :: key
    type(t_node), pointer :: node, prev, next, above
    integer :: i, level
    node => search_node(this, key)
    if (node%key /= key) return
    this%size = this%size - 1
    level = 0
    prev => node%prev
    next => node%next
    prev%next => next
    next%prev => prev
    do
      above => node%above
      deallocate(node)
      node => above
      level = level + 1
      if (.not.associated(node)) exit
      do while (.not.associated(prev%above))
        prev => prev%prev
      end do
      prev => prev%above
      prev%length = prev%length + node%length - 1
      next => node%next
      prev%next => next
      next%prev => prev
    end do
    do i = level, this%level
      do while (.not.associated(prev%above))
        prev => prev%prev
      end do
      prev => prev%above
      prev%length = prev%length - 1
    end do
  end

  ! idx 番目のノードを削除
  subroutine remove_key_at(this, idx)
    class(t_skip_list_map), intent(inout) :: this
    integer, intent(in) :: idx
    integer :: key
    if (idx < 1 .or. idx > this%size) return
    key = get_key_at(this, idx)
    call remove(this, key)
  end

  ! 先頭のキーを取得
  integer(key_kind) function first_key(this) result(res)
    class(t_skip_list_map), intent(in) :: this
    res = merge(-infty, this%head%next%key, this%size == 0)
  end

  ! 先頭のキーを取得し、そのノードを削除する
  integer(key_kind) function poll_first(this) result(res)
    class(t_skip_list_map), intent(inout) :: this
    type(t_node), pointer :: node
    res = merge(-infty, this%head%next%key, this%size == 0)
    if (this%size > 0) call remove(this, res)
  end

  ! 末尾のキーを取得
  integer(key_kind) function last_key(this) result(res)
    class(t_skip_list_map), intent(in) :: this
    res = merge(infty, this%tail%prev%key, this%size == 0)
  end

  ! 末尾のキーを取得し、そのノードを削除する
  integer(key_kind) function poll_last(this) result(res)
    class(t_skip_list_map), intent(inout) :: this
    res = merge(infty, this%tail%prev%key, this%size == 0)
    if (this%size > 0) call remove(this, res)
  end

  ! key >= node%key となるノードのキーを取得
  integer(key_kind) function floor_key(this, key) result(res)
    class(t_skip_list_map), intent(in) :: this
    integer(key_kind), intent(in) :: key
    type(t_node), pointer :: node
    node => search_node(this, key)
    res = node%key
  end

  ! key > node%key となるノードのキーを取得
  integer(key_kind) function lower_key(this, key) result(res)
    class(t_skip_list_map), intent(in) :: this
    integer(key_kind), intent(in) :: key
    type(t_node), pointer :: node
    node => search_node(this, key)
    res = merge(node%prev%key, node%key, node%key == key)
  end

  ! key <= node%key となるノードのキーを取得
  integer(key_kind) function ceiling_key(this, key) result(res)
    class(t_skip_list_map), intent(in) :: this
    integer(key_kind), intent(in) :: key
    type(t_node), pointer :: node
    node => search_node(this, key)
    res = merge(node%key, node%next%key, node%key == key)
  end

  ! key < node%key となるノードのキーを取得
  integer(key_kind) function higher_key(this, key) result(res)
    class(t_skip_list_map), intent(in) :: this
    integer(key_kind), intent(in) :: key
    type(t_node), pointer :: node
    node => search_node(this, key)
    res = node%next%key
  end
end module mod_skip_list_map

動作が正しいかどうかは以下の問題で確認しました。
AtCoder Beginner Contest 177 F - I hate Shortest Path Problem
ハッシュマップのソースコードはこちらに置いてます。また、map を作れたらキーに対応する値の変数部分を削るだけで set も実装できるので、skip_list_set.f08 も作成しました。
set バージョンは以下の問題で動作確認しました。
AtCoder Regular Contest 033 C - データ構造

参考

スキップリスト（Wikipedia）
SkipLists.pdf
FPGA開発日記 - SkipListの勉強
A Skip List Cookbook