Fortranとセグメント木と抽象化された型

セグメント木なので競技プログラミングの話...と見せかけて**Fortranの抽象構造型**の話をします.

# セグメント木とは?
競技プログラミングで良く目にする**データ構造**です.
区間に対するある種の操作と配列の値の変更ができます.

## 例え話
配列 `[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]` の区間`2〜8`番目の和を求めるとします.
その場合, (区間`2〜2`の和, 区間`3〜4`の和, 区間`5〜8`の和)の和が分かれば計算できます.
つまり, あらかじめ`上手い感じに選ばれた区間`の和を計算しておくことで, `任意の区間`の和を計算することができます.
更に, $i$ 番目の値の変更も素早くできます.
何故なら, `上手い感じに選ばれた区間` には $i$ に関する区間が少ないからです.
(上の例だと, `2`番目に対応するのは `1〜16`, `1〜8`, `1〜4`, `1〜2`, `2〜2` の5個しかありません. これは $N = 16$ として, $\log(N) = 4$ くらいの値です.

詳細に知りたい場合はグーグルで検索すると図付きで解説されているページがたくさん出てきます.

## 条件
- 2項演算子がモノイド(結合則と単位元を持つ)
  + 要するに, 和なら $(a + b) + c == a + (b + c)$, $a + 0 = 0 + a = a$.
  + 最大値, 最小値でも似たものを作れるので試してみましょう.

だけですね.
この性質があればセグメント木に乗せられます.

# Fortran で実装するには?
一応, [私の実装](https://github.com/osada-yum/Fortran_competitive_library/blob/master/src/segment_tree_m.fypp)があります.
セグメント木の作り方は他の言語と同じなので省略します.
モノイドの実装方法は他の言語と違うと思います.

## モノイドをどのように表現するか
例えば, [AtCoderの実装](https://github.com/atcoder/ac-library) では`template` によって関数と単位元を渡しています.
Fortran では`template`は使えないので, 他の方法を考える必要があります.
私の実装では, [抽象構造型](https://www.nag-j.co.jp/fortran/fortran2003/Fortran2003_2_5.html) を使っています.

## [抽象構造型](https://www.nag-j.co.jp/fortran/fortran2003/Fortran2003_2_5.html)
`abstract interface` ([抽象引用仕様](https://www.nag-j.co.jp/fortran/fortran2003/Fortran2003_8_1.html)) を使うことで実装できます.

<details open><summary>monoid_op_m.f90</summary>

```fortran:monoid_op_m.f90
module monoid_op_m
  use, intrinsic :: iso_fortran_env
  implicit none
  type, abstract :: monoid_op_int32
     private
   contains
     procedure(identity_monoid_op_int32), nopass, deferred :: identity
     procedure(bin_op_monoid_op_int32), nopass, deferred :: bin_op
  end type monoid_op_int32
  abstract interface
     pure integer(int32) function identity_monoid_op_int32() result(res)
       import int32
     end function identity_monoid_op_int32
     pure integer(int32) function bin_op_monoid_op_int32(x, y) result(res)
       import int32
       integer(int32), intent(in) :: x, y
     end function bin_op_monoid_op_int32
  end interface
  public :: plus_monoid_op_int32
  type, extends(monoid_op_int32) :: plus_monoid_op_int32
     private
   contains
     procedure, nopass :: identity => identity_plus_monoid_op_int32
     procedure, nopass :: bin_op   => bin_op_plus_monoid_op_int32
  end type plus_monoid_op_int32
  public :: min_monoid_op_int32
  type, extends(monoid_op_int32) :: min_monoid_op_int32
     private
   contains
     procedure, nopass :: identity => identity_min_monoid_op_int32
     procedure, nopass :: bin_op   => bin_op_min_monoid_op_int32
  end type min_monoid_op_int32
contains
  pure integer(int32) function identity_plus_monoid_op_int32() result(res)
    res = 0_int32
  end function identity_plus_monoid_op_int32
  pure integer(int32) function bin_op_plus_monoid_op_int32(x, y) result(res)
    integer(int32), intent(in) :: x, y
    res = x + y
  end function bin_op_plus_monoid_op_int32
  pure integer(int32) function identity_min_monoid_op_int32() result(res)
    res = huge(0_int32)
  end function identity_min_monoid_op_int32
  pure integer(int32) function bin_op_min_monoid_op_int32(x, y) result(res)
    integer(int32), intent(in) :: x, y
    res = min(x, y)
  end function bin_op_min_monoid_op_int32
end module monoid_op_m
```

</details>

<details open><summary> test_monoid_op.f90 </summary>

```fortran:test_monoid_op.f90
program test_monoid_op
  use, intrinsic :: iso_fortran_env
  use monoid_op_m
  implicit none
  class(monoid_op_int32), allocatable :: myop
  allocate(myop, source = plus_monoid_op_int32())
  write(output_unit, '(i0)') myop%identity() ! 単位元.
  write(output_unit, '(i0)') myop%bin_op(myop%identity(), 100) ! 単位元.
  write(output_unit, '(i0)') myop%bin_op(myop%bin_op(3, 4), 5) ! 結合則.
  write(output_unit, '(i0)') myop%bin_op(3, myop%bin_op(4, 5)) ! 結合則.
  deallocate(myop)
  allocate(myop, source = min_monoid_op_int32())
  write(output_unit, '(i0)') myop%identity()
  write(output_unit, '(i0)') myop%bin_op(myop%identity(), 100)
  write(output_unit, '(i0)') myop%bin_op(myop%bin_op(3, 4), 5)
  write(output_unit, '(i0)') myop%bin_op(3, myop%bin_op(4, 5))
end program test_monoid_op
```

</details>

- 実行結果
```bash
$ gfortran --version
GNU Fortran (Ubuntu 12.3.0-1ubuntu1~22.04) 12.3.0
Copyright (C) 2022 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
$ gfortran monoid_op_m.f90 test_monoid_op.f90 && ./a.out
0
100
12
12
2147483647
100
3
3
```

単位元や結合則を満たしているように見えます.


# 蛇足

<details><summary> char_monoid.f90 </summary>

```fortran:char_monoid.f90
module vec_string_m
  use, intrinsic :: iso_fortran_env
  implicit none
  private
  integer(int32), public :: vec_string_iostat
  character(len=100), public :: vec_string_iomsg = ""
  public vec_string
  type :: vec_string
     private
     character, allocatable, public :: str_(:)
     integer(int32) :: size_ = 0_int32, capa_ = 0_int32
   contains
     procedure, pass :: init => init_vec_string
     procedure, pass :: with_capacity => with_capacity_vec_string
     procedure, pass :: resize => resize_vec_string
     procedure, pass :: push_back => push_back_vec_string
     procedure, pass :: pop_back => pop_back_vec_string
     procedure, pass :: back => back_vec_string
     procedure, pass :: size => size_vec_string
     procedure, pass, private :: assign => assign_vec_string
     generic :: assignment(=) => assign
     procedure, pass :: read => read_vec_string
     procedure, pass, private :: formatted_read => formatted_read_vec_string
     generic :: read(formatted) => formatted_read
     ! final :: destroy_vec_string
  end type vec_string
  interface vec_string
     module procedure :: construct_vec_string, construct_by_size_vec_string, construct_by_str_vec_string
  end interface vec_string
  public destroy_vec_string
  public :: compare, operator(<), operator(<=), operator(>), operator(>=), operator(==), operator(/=)
  interface compare
     module procedure :: compare_vec_string
  end interface compare
  interface operator(<)
     module procedure :: less_vec_string
  end interface operator(<)
  interface operator(<=)
     module procedure :: less_equal_vec_string
  end interface operator(<=)
  interface operator(>)
     module procedure :: greater_vec_string
  end interface operator(>)
  interface operator(>=)
     module procedure :: greater_equal_vec_string
  end interface operator(>=)
  interface operator(==)
     module procedure :: equal_vec_string
  end interface operator(==)
  interface operator(/=)
     module procedure :: not_equal_vec_string
  end interface operator(/=)
contains
  pure type(vec_string) function construct_vec_string() result(res)
    call res%with_capacity(1)
  end function construct_vec_string
  pure type(vec_string) function construct_by_size_vec_string(n) result(res)
    integer(int32), intent(in) :: n
    call res%with_capacity(n)
  end function construct_by_size_vec_string
  pure type(vec_string) function construct_by_str_vec_string(str) result(res)
    character(len=*), intent(in) :: str
    integer(int32) :: n, capa
    integer(int32) :: i
    n = len_trim(str)
    capa = 1_int32
    do while (capa < n)
       capa = capa * 2
    end do
    call res%with_capacity(capa)
    res%size_ = n
    do i = 1, n
       res%str_(i) = str(i:i)
    end do
  end function construct_by_str_vec_string

  pure subroutine destroy_vec_string(this)
    type(vec_string), intent(inout) :: this
    if (allocated(this%str_)) &
         deallocate(this%str_)
    this%size_ = 0
    this%capa_ = 0
  end subroutine destroy_vec_string

  pure subroutine assign_vec_string(lhs, rhs)
    class(vec_string), intent(inout) :: lhs
    class(vec_string), intent(in) :: rhs
    call destroy_vec_string(lhs)
    allocate(lhs%str_, source = rhs%str_(:))
    lhs%size_ = rhs%size_
    lhs%capa_ = rhs%capa_
  end subroutine assign_vec_string
  pure subroutine init_vec_string(this)
    class(vec_string), intent(inout) :: this
    call this%with_capacity(1)
  end subroutine init_vec_string
  pure subroutine with_capacity_vec_string(this, n)
    class(vec_string), intent(inout) :: this
    integer(int32), intent(in) :: n
    character, allocatable :: tmp(:)
    integer(int32) :: capa
    if (n < 1) return
    capa = 1_int32
    do while (capa < n)
       capa = capa * 2
    end do
    !> capa >= n.
    ! write(error_unit, '(a, i0, " < ", i0)') "hi", this%capa_, capa
    if (this%capa_ >= capa) return
    !> this%capa_ < capa
    this%capa_ = capa
    if (allocated(this%str_)) then
       allocate(tmp(capa))
       tmp(1:this%size()) = this%str_(1:this%size())
       call move_alloc(from = tmp, to = this%str_)
    else
       allocate(this%str_(capa))
    end if
  end subroutine with_capacity_vec_string
  pure subroutine resize_vec_string(this, n)
    class(vec_string), intent(inout) :: this
    integer(int32), intent(in) :: n
    if (n < 0) return
    call this%with_capacity(n)
    this%size_ = n
  end subroutine resize_vec_string
  pure subroutine push_back_vec_string(this, c)
    class(vec_string), intent(inout) :: this
    character, intent(in) :: c
    if (this%size() == this%capa_) &
         call this%with_capacity(max(1, 2 * this%capa_))
    this%size_ = this%size_ + 1
    this%str_(this%size()) = c
  end subroutine push_back_vec_string
  pure subroutine pop_back_vec_string(this)
    class(vec_string), intent(inout) :: this
    this%size_ = this%size_ - 1
  end subroutine pop_back_vec_string
  pure character function back_vec_string(this) result(res)
    class(vec_string), intent(in) :: this
    res = this%str_(this%size())
  end function back_vec_string
  pure integer(int32) function size_vec_string(this) result(res)
    class(vec_string), intent(in) :: this
    res = this%size_
  end function size_vec_string
  impure subroutine read_vec_string(this, unit)
    class(vec_string), intent(inout) :: this
    integer(int32), intent(in) :: unit
    character :: c
    ! write(error_unit, '(a)') this%str_(1:this%size())
    call this%resize(0)
    ! write(error_unit, '(a)') "hi"
    do
       read(unit, '(a1)', advance = "no", iostat = vec_string_iostat, iomsg = vec_string_iomsg) c
       ! write(error_unit, *) this%size(), ": ", c, ", ", this%str_(1:this%size())
       select case(vec_string_iostat)
       case(iostat_end)
          if (this%size() == 0) error stop "End of file in reading Symbol’s value as variable is void: vec_string."
          vec_string_iomsg = "End of file in reading Symbol’s value as variable is void: vec_string."
          exit
       case(iostat_eor)
          if (this%size() == 0) cycle
          vec_string_iostat = 0
          vec_string_iomsg = ""
          ! write(error_unit, '(a)') "End of record"
          exit
       case(0)
          if (c == " ") then
             if (this%size() == 0) cycle
             exit
          end if
       case default
          ! vec_string_iomsg = "Unknown iostat in reading Symbol’s value as variable is void: vec_string."
          return
       end select
       call this%push_back(c)
    end do
    ! write(error_unit, '(a)') this%str_(1:this%size())
  end subroutine read_vec_string
  impure subroutine formatted_read_vec_string(this, unit, iotype, vlist, iostat, iomsg)
    class(vec_string), intent(inout) :: this
    integer(int32), intent(in) :: unit
    character(len=*), intent(in) :: iotype
    integer(int32), intent(in) :: vlist(:)
    integer(int32), intent(out) :: iostat
    character(len=*), intent(inout) :: iomsg
    character :: c
    associate(dummy => iotype, dummy2 => vlist)
    end associate
    ! write(error_unit, '(a, 2(i0, 1x))') "loc: ", loc(this), loc(this%str_)
    ! write(error_unit, '(a, i0)') "loc: ", loc(iomsg)
    ! write(error_unit, '(*(a1))') "|", this%str_(1:this%size()), "|"
    ! call destroy_vec_string(this)
    ! this = vec_string(1)
    ! write(error_unit, '(a, *(1x, i0))') iotype, vlist(:)
    call this%resize(0)
    do
       read(unit, '(a1)', advance = "no", iostat = iostat) c
       ! write(error_unit, *) this%size(), ": ", c, ", ", this%str_(1:this%size())
       ! write(error_unit, '(2(i0, 1x), a, i0, 1x, a)') iostat_end, iostat_eor, ", ", iostat, "|"//c//"| "//trim(iomsg)
       select case(iostat)
       case(iostat_end)
          if (this%size() == 0) error stop "End of file in reading Symbol’s value as variable is void: vec_string."
          iomsg = "End of file in reading Symbol’s value as variable is void: vec_string."
          exit
       case(iostat_eor)
          if (this%size() == 0) cycle
          ! iostat = 0
          iomsg = "hi"
          ! backspace(unit)
          exit
       case(0)
          if (c == " ") then
             if (this%size() == 0) cycle
             exit
          end if
       case default
          iomsg = "Unknown iostat in reading Symbol’s value as variable is void: vec_string."
          return
       end select
       call this%push_back(c)
    end do
  end subroutine formatted_read_vec_string

  pure integer(int32) function compare_vec_string(lhs, rhs) result(res)
    type(vec_string), intent(in) :: lhs, rhs
    integer(int32) :: i
    associate(ls => lhs%size(), rs => rhs%size())
      do i = 1, min(ls, rs)
         associate(lc => lhs%str_(i), rc => rhs%str_(i))
           if (lc == rc) cycle
           if (lc < rc) then
              res = -1; return
           else
              res = 1; return
           end if
         end associate
      end do
      if (ls == rs) then
         res = 0
      else if (ls < rs) then
         res = -1
      else
         res = 1
      end if
    end associate
  end function compare_vec_string
  pure logical function less_vec_string(lhs, rhs) result(res)
    type(vec_string), intent(in) :: lhs, rhs
    res = compare(lhs, rhs) < 0
  end function less_vec_string
  pure logical function less_equal_vec_string(lhs, rhs) result(res)
    type(vec_string), intent(in) :: lhs, rhs
    res = compare(lhs, rhs) <= 0
  end function less_equal_vec_string
  pure logical function greater_vec_string(lhs, rhs) result(res)
    type(vec_string), intent(in) :: lhs, rhs
    res = compare(lhs, rhs) > 0
  end function greater_vec_string
  pure logical function greater_equal_vec_string(lhs, rhs) result(res)
    type(vec_string), intent(in) :: lhs, rhs
    res = compare(lhs, rhs) >= 0
  end function greater_equal_vec_string
  pure logical function equal_vec_string(lhs, rhs) result(res)
    type(vec_string), intent(in) :: lhs, rhs
    integer(int32) :: i
    res = .false.
    if (lhs%size() /= rhs%size()) return
    do i = 1, lhs%size()
       if (lhs%str_(i) /= rhs%str_(i)) return
    end do
    res = .true.
  end function equal_vec_string
  pure logical function not_equal_vec_string(lhs, rhs) result(res)
    type(vec_string), intent(in) :: lhs, rhs
    res = .not. (lhs == rhs)
  end function not_equal_vec_string
end module vec_string_m

module segment_tree_m
  use, intrinsic :: iso_fortran_env
  use vec_string_m
  implicit none
  private
  public :: segment_tree_vec_string
  public :: monoid_op_vec_string
  type :: segment_tree_vec_string
     private
     integer(int32) :: arr_size_, tree_size_, depth_
     type(vec_string), allocatable :: arr_(:)
     class(monoid_op_vec_string), allocatable :: monoid_
   contains
     procedure, pass :: init   => init_segment_tree_vec_string

     procedure, pass :: update => update_segment_tree_vec_string
     procedure, pass :: query => query_segment_tree_vec_string
  end type segment_tree_vec_string
  type, abstract :: monoid_op_vec_string
     private
   contains
     procedure(identity_vec_string), nopass, deferred :: identity
     procedure(bin_op_vec_string)  , nopass, deferred :: bin_op
  end type monoid_op_vec_string
  abstract interface

     pure type(vec_string) function identity_vec_string() result(res)
  import vec_string
     end function identity_vec_string
     pure type(vec_string) function bin_op_vec_string(x, y) result(res)
  import vec_string
       type(vec_string), intent(in) :: x, y
     end function bin_op_vec_string
  end interface
contains
  !> indices rage [1:2^a-1]
  subroutine init_segment_tree_vec_string (this, num_elems, monoid_)
    class(segment_tree_vec_string), intent(inout) :: this
    integer(int32), intent(in) :: num_elems
    class(monoid_op_vec_string), intent(in) :: monoid_
    integer(int32) :: tree_size
    allocate(this%monoid_, source = monoid_)
    tree_size = 1
    this%depth_ = 1
    do while (tree_size < num_elems)
       tree_size = tree_size * 2
       this%depth_ = this%depth_ + 1
    end do
    this%tree_size_ = tree_size
    this%arr_size_ = 2 * tree_size - 1
    allocate(this%arr_(this%arr_size_), source = this%monoid_%identity())
  end subroutine init_segment_tree_vec_string

  subroutine update_segment_tree_vec_string (this, idx, val)
    class(segment_tree_vec_string), intent(inout) :: this
    integer(int32), intent(in) :: idx
    type(vec_string), intent(in) :: val
    integer(int32) :: pos
    pos = this%tree_size_ + idx - 1
    this%arr_(pos) = val
    do while (pos > 1)
       pos = pos / 2
       this%arr_(pos) = this%monoid_%bin_op(this%arr_(2 * pos), this%arr_(2 * pos + 1))
    end do
  end subroutine update_segment_tree_vec_string
  !> close interval [a, b].
  pure type(vec_string) function query_segment_tree_vec_string (this, a, b) result(res)
    class(segment_tree_vec_string), intent(in) :: this
    integer(int32), intent(in) :: a, b
    type(vec_string) res_r
    integer(int32) :: p, r
    res = this%monoid_%identity()
    res_r = this%monoid_%identity()
    if (a > b) return
    p = this%tree_size_ + a - 1
    r = this%tree_size_ + b - 1
    do while (p <= r)
       if (iand(p, b'1') == 1) then !> right child.
          res = this%monoid_%bin_op(res, this%arr_(p))
          p = p + 1
       end if
       if (iand(r, b'1') == 0) then !> left child.
          res_r = this%monoid_%bin_op(this%arr_(r), res_r)
          r = r - 1
       end if
       p = p / 2
       r = r / 2
    end do
    res = this%monoid_%bin_op(res, res_r)
  end function query_segment_tree_vec_string
end module segment_tree_m

module mymonoid_m
  use, intrinsic :: iso_fortran_env
  use segment_tree_m
  use vec_string_m
  implicit none
  public plus_op_vec_string
  type, extends(monoid_op_vec_string) :: plus_op_vec_string
     private
   contains
     procedure, nopass :: identity => identity_plus_op_vec_string
     procedure, nopass :: bin_op => bin_op_plus_op_vec_string
  end type plus_op_vec_string
contains
  pure type(vec_string) function identity_plus_op_vec_string() result(res)
    res = vec_string("")
  end function identity_plus_op_vec_string
  pure type(vec_string) function bin_op_plus_op_vec_string(x, y) result(res)
    type(vec_string), intent(in) :: x, y
    integer(int32) :: i
    do i = 1, x%size()
       call res%push_back(x%str_(i))
    end do
    do i = 1, y%size()
       call res%push_back(y%str_(i))
    end do
  end function bin_op_plus_op_vec_string
end module mymonoid_m
```

</details>

```fortran:christmas.f90
program christmas
  use, intrinsic :: iso_fortran_env
  use vec_string_m
  use segment_tree_m
  use mymonoid_m
  implicit none
  character(len=*), parameter :: s = "abchristmashellomerry"
  type(segment_tree_vec_string) :: st
  integer(int32) :: i
  call st%init(len_trim(s), plus_op_vec_string())
  do i = 1, len_trim(s)
     call st%update(i, vec_string(s(i:i)))
  end do
  associate(res => st%query(12, len_trim(s)))
    write(output_unit, '(*(a1))') res%str_(1:res%size())
  end associate
  associate(res => st%query(3, 11))
    write(output_unit, '(*(a1))') res%str_(1:res%size())
  end associate
end program christmas
```

- 実行結果
```bash
 $ gfortran char_monoid.f90 christmas.f90 && ./a.out
hellomerry
christmas
```

というわけで, 文字列の結合はモノイドみたいです. (ちなみに計算量が増えているのでセグメント木に乗せる利点はありません.)
おあとがよろしいようで.


# 参考
## セグメント木
https://algo-logic.info/segment-tree/

https://github.com/atcoder/ac-library

https://github.com/osada-yum/Fortran_competitive_library/blob/master/src/segment_tree_m.fypp

## Fortran
https://www.nag-j.co.jp/fortran/fortran2003/Fortran2003_8_1.html

https://www.nag-j.co.jp/fortran/fortran2003/Fortran2003_2_5.html