Fortran で priority queue を作ってみた

動機

priority queue は Dijkstra 法などに用いられ、コレクション内の最小・最大の要素を高速に取り出すことができる便利なデータ構造です。しかしながら Fortran には実装されていないので、自分で実装します。今回は integer を要素として持ち、コレクション内の最小値を取り出せるものを作成します。

priority queue とは

priority queue は追加する要素に優先度をつけ、取り出す際には優先度の高いものから取り出すことのできるデータ構造です。実装にはヒープ構造を用いるものが多いようです。 priority queue に要素を追加するのは $O(\mathrm{log}N)$ 、priority queue から要素を取り出すのは $O(\mathrm{log}N)$ の時間で行うことができます。

環境

macOS Sierra (10.12.6)
GNU Fortran (GCC) 6.3.0

プログラム

実際に実装したのが以下のものです。

priority_queue.f08

module mod_priority_queue
  implicit none

  type t_priority_queue
    integer :: num = 0
    integer, pointer :: heap(:) => null()
  end type t_priority_queue

contains

  subroutine offer(pq,item)
    implicit none
    type(t_priority_queue), intent(inout) :: pq
    integer, intent(in) :: item
    integer :: n, i, t
    integer, allocatable :: tmp(:)

    if (.not.associated(pq%heap)) allocate(pq%heap(1))
    if (pq%num == size(pq%heap)) then
      allocate(tmp(pq%num))
      tmp = pq%heap
      deallocate(pq%heap)
      allocate(pq%heap(2*pq%num))
      pq%heap(1:pq%num) = tmp
      deallocate(tmp)
    end if

    pq%num = pq%num+1
    pq%heap(pq%num) = item

    n = pq%num
    do while (n > 1)
      i = n/2
      if (pq%heap(n) < pq%heap(i)) then
        t = pq%heap(n)
        pq%heap(n) = pq%heap(i)
        pq%heap(i) = t
      end if
      n = i
    end do
    return
  end subroutine offer

  subroutine clear(pq)
    implicit none
    type(t_priority_queue), intent(inout) :: pq

    if (associated(pq%heap)) deallocate(pq%heap)
    pq%num = 0
    return
  end subroutine clear

  function poll(pq) result(item)
    implicit none
    type(t_priority_queue), intent(inout) :: pq
    integer :: item, n, i, j, tmp

    n = pq%num
    item = pq%heap(1)
    pq%heap(1) = pq%heap(n)
    pq%num = pq%num-1

    i = 1
    do while (2*i < n)
      j = 2*i
      if (j+1 < n .and. pq%heap(j+1) < pq%heap(j)) j = j+1
      if (pq%heap(j) < pq%heap(i)) then
        tmp = pq%heap(j)
        pq%heap(j) = pq%heap(i)
        pq%heap(i) = tmp
      end if
      i = j
    end do
    return
  end function poll

  function peek(pq) result(item)
    implicit none
    type(t_priority_queue), intent(inout) :: pq
    integer :: item

    item = pq%heap(1)
    return
  end function peek

end module mod_priority_queue

また、実装できているか確認するために以下のプログラムも書きました。

test_priority_queue.f08

program test_priority_queue
  use mod_priority_queue
  implicit none
  integer, parameter :: n = 6
  type(t_priority_queue) :: pq
  integer :: x(n), i

  call random_seed_clock()
  call randint(x,10*n)

  do i = 1, n
    call offer(pq,x(i))
    write(*,'("offer: ",i0," size: ",i0)') x(i), pq%num
  end do

  do while (pq%num > 0)
    write(*,'("poll: ",i0," size: ",i0)') poll(pq), pq%num
  end do

  call randint(x,10*n)

  do i = 1, n
    call offer(pq,x(i))
    write(*,'("offer: ",i0," size: ",i0)') x(i), pq%num
  end do

  write(*,'("peek: ",i0," size: ",i0)') peek(pq), pq%num

  call clear(pq)
  write(*,'("clear, size: ",i0)') pq%num

  call randint(x,10*n)

  do i = 1, n
    call offer(pq,x(i))
    write(*,'("offer: ",i0," size: ",i0)') x(i), pq%num
  end do

  write(*,'("peek: ",i0," size: ",i0)') peek(pq), pq%num

  stop
contains

  subroutine random_seed_clock()
    implicit none
    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)
    return
  end subroutine random_seed_clock

  subroutine randint(x,maxint)
    implicit none
    integer, intent(out) :: x(:)
    integer, intent(in) :: maxint
    real(8) :: r(size(x))

    call random_number(r)
    x = int(maxint*r)
    return
  end subroutine randint

end program test_priority_queue

実行結果は以下のとおりです。

offer: 22 size: 1
offer: 19 size: 2
offer: 37 size: 3
offer: 40 size: 4
offer: 56 size: 5
offer: 50 size: 6
poll: 19 size: 5
poll: 22 size: 4
poll: 37 size: 3
poll: 40 size: 2
poll: 50 size: 1
poll: 56 size: 0
offer: 43 size: 1
offer: 57 size: 2
offer: 41 size: 3
offer: 23 size: 4
offer: 2 size: 5
offer: 50 size: 6
peek: 2 size: 6
clear, size: 0
offer: 8 size: 1
offer: 43 size: 2
offer: 6 size: 3
offer: 8 size: 4
offer: 29 size: 5
offer: 38 size: 6
peek: 6 size: 6

これで逐次的に最小値を取り出せるデータ構造が手に入りました。

追記 2019/08/11

不必要な部分を削り、よりオブジェクト指向な書き方に変更しました (ソースコード) 。