module mod_wavelet_matrix
  implicit none
  integer, private, parameter :: intkind = 4
  integer, private, parameter :: bitsize = 64
  integer(8), private, parameter :: mask = -1
  type space_saving_array
    integer :: n, size, offset
    integer(8), allocatable :: bit(:)
    integer, allocatable :: cnt(:)
  end type
  type wavelet_matrix
    integer :: n, m
    type(space_saving_array), allocatable :: bit(:)
    integer, private, allocatable :: order(:), start(:)
    integer(intkind), private, allocatable :: vals(:)
  contains
    procedure :: access => access_to
    procedure :: rank => rank_of
    procedure :: select => select_from
    procedure :: quantile => quantile
    procedure :: range_less => range_less
    procedure :: range_freq => range_freq
  end type
  interface wavelet_matrix
    module procedure :: new_wm
  end interface
  private :: space_saving_array, new_ssa, access_ssa, rank_ssa, select_ssa
  private :: find
contains
  type(space_saving_array) function new_ssa(a) result(res)
    logical, intent(in) :: a(:)
    integer :: n, m, i, j, k
    n = size(a)
    m = (n + bitsize - 1) / bitsize
    res%n = n
    res%size = m
    allocate(res%bit(m), res%cnt(0:m))
    res%bit = 0
    j = 1
    k = 0
    do i = 1, n
      if (a(i)) res%bit(j) = ibset(res%bit(j), k)
      k = k + 1
      if (k == bitsize) then
        j = j + 1
        k = 0
      end if
    end do
    res%cnt(0) = 0
    do i = 1, m
      res%cnt(i) = res%cnt(i - 1) + popcnt(res%bit(i))
    end do
    res%offset = n - res%cnt(m)
  end
  logical function access_ssa(this, idx) result(res)
    type(space_saving_array), intent(in) :: this
    integer, intent(in) :: idx
    integer :: i, j
    i = (idx - 1) / bitsize + 1
    j = idx - 1 - (i - 1) * bitsize
    res = btest(this%bit(i), j)
  end
  integer function rank_ssa(this, idx, bool) result(res)
    type(space_saving_array), intent(in) :: this
    integer, intent(in) :: idx
    logical, intent(in) :: bool
    integer :: i, j
    i = (idx - 1) / bitsize + 1
    j = idx - (i - 1) * bitsize
    ! res = this%cnt(i - 1) + popcnt(and(not(lshift(mask, j)), this%bit(i)))
    res = this%cnt(i) - popcnt(and(lshift(mask, j), this%bit(i)))
    if (bool) return
    res = idx - res
  end
  integer function select_ssa(this, ord, bool) result(res)
    type(space_saving_array), intent(in) :: this
    integer, intent(in) :: ord
    logical, intent(in) :: bool
    integer :: l, r, m, s
    res = -1
    s = merge(this%cnt(this%size), this%n - this%cnt(this%size), bool)
    if (s < ord) return
    l = 0
    r = this%size + 1
    do while (r - l > 1)
      m = (l + r) / 2
      s = merge(this%cnt(m), bitsize * m - this%cnt(m), bool)
      if (s < ord) then
        l = m
      else
        r = m
      end if
    end do
    l = bitsize * (r - 1)
    r = bitsize * r + 1
    do while (r - l > 1)
      m = (l + r) / 2
      if (rank_ssa(this, m, bool) < ord) then
        l = m
      else
        r = m
      end if
    end do
    res = r
  end
  type(wavelet_matrix) function new_wm(a, sigma) result(res)
    integer(intkind), intent(in) :: a(:)
    integer(intkind), optional, intent(in) :: sigma
    integer(intkind) :: maxa, x(size(a)), y(size(a)), z(size(a))
    integer, allocatable :: s(:), t(:), u(:)
    integer :: n, m, nx, ny, i, j
    logical :: bool(size(a))
    maxa = merge(sigma, maxval(a), present(sigma))
    n = size(a)
    m = 8 * intkind - 1
    do while (.not.btest(maxa, m))
      m = m - 1
    end do
    res%n = n
    res%m = m
    allocate(res%bit(0:m))
    z = a
    do j = 0, m
      nx = 0
      ny = 0
      do i = 1, n
        bool(i) = btest(z(i), m - j)
        if (bool(i)) then
          ny = ny + 1
          y(ny) = z(i)
        else
          nx = nx + 1
          x(nx) = z(i)
        end if
      end do
      res%bit(j) = new_ssa(bool)
      z(1:nx) = x(1:nx)
      z(nx + 1:n) = y(1:ny)
    end do
    nx = n
    n = 1
    do i = 2, nx
      if (z(i) /= z(i - 1)) n = n + 1
    end do
    allocate(s(n), t(n), u(n), res%order(n), res%start(n), res%vals(n))
    u = [(i, i = 1, n)]
    n = 1
    res%start(1) = 1
    do i = 2, nx
      if (z(i) /= z(n)) then
        n = n + 1
        z(n) = z(i)
        res%start(n) = i
      end if
    end do
    do j = 0, m
      nx = 0
      ny = 0
      do i = 1, n
        if (btest(z(i), j)) then
          ny = ny + 1
          y(ny) = z(i)
          t(ny) = u(i)
        else
          nx = nx + 1
          x(nx) = z(i)
          s(nx) = u(i)
        end if
      end do
      z(1:nx) = x(1:nx)
      z(nx + 1:n) = y(1:ny)
      u(1:nx) = s(1:nx)
      u(nx + 1:n) = t(1:ny)
    end do
    res%order = u
    res%vals = z(1:n)
  end
  integer function find(this, val) result(res)
    type(wavelet_matrix), intent(in) :: this
    integer(intkind) :: val
    integer :: l, r
    l = 0
    r = size(this%vals) + 1
    do while (r - l > 1)
      res = (l + r) / 2
      if (this%vals(res) < val) then
        l = res
      else if (this%vals(res) > val) then
        r = res
      else
        return
      end if
    end do
    res = -1
  end
  integer(intkind) function access_to(this, idx) result(res)
    class(wavelet_matrix), intent(in) :: this
    integer, intent(in) :: idx
    logical :: bool
    integer :: i, j
    i = idx
    res = 0
    do j = 0, this%m
      bool = access_ssa(this%bit(j), i)
      i = rank_ssa(this%bit(j), i, bool)
      if (bool) then
        res = ibset(res, this%m - j)
        i = i + this%bit(j)%offset
      end if
    end do
  end
  integer function rank_of(this, idx, val) result(res)
    class(wavelet_matrix), intent(in) :: this
    integer, intent(in) :: idx
    integer(intkind), intent(in) :: val
    logical :: bool
    integer :: i, j, k
    res = 0
    k = find(this, val)
    if (k < 0) return
    k = this%order(k)
    i = idx
    do j = 0, this%m
      bool = btest(val, this%m - j)
      i = rank_ssa(this%bit(j), i, bool)
      if (bool) i = i + this%bit(j)%offset
    end do
    res = i - this%start(k) + 1
  end
  integer function select_from(this, ord, val) result(res)
    class(wavelet_matrix), intent(in) :: this
    integer, intent(in) :: ord
    integer(intkind), intent(in) :: val
    logical :: bool
    integer :: i, j, k, l
    res = -1
    k = find(this, val)
    if (k < 0) return
    k = this%order(k)
    i = this%start(k) + ord - 1
    if (k == size(this%start)) then
      if (i > this%n) return
    else
      if (i >= this%start(k + 1)) return
    end if
    do j = this%m, 0, -1
      bool = btest(val, this%m - j)
      if (bool) i = i - this%bit(j)%offset
      i = select_ssa(this%bit(j), i, bool)
    end do
    res = i
  end
  integer(intkind) function quantile(this, a, b, ord) result(res)
    class(wavelet_matrix), intent(in) :: this
    integer, intent(in) :: a, b, ord
    logical :: bool
    integer :: l, r, i, j, d
    res = -1
    if (b - a + 1 < ord) return
    l = a
    r = b
    i = ord
    res = 0
    do j = 0, this%m
      d = (r - l + 1) - (rank_ssa(this%bit(j), r, .true.) - rank_ssa(this%bit(j), l - 1, .true.))
      bool = i > d
      l = rank_ssa(this%bit(j), l - 1, bool) + 1
      r = rank_ssa(this%bit(j), r, bool)
      if (bool) then
        res = ibset(res, this%m - j)
        i = i - d
        l = l + this%bit(j)%offset
        r = r + this%bit(j)%offset
      end if
    end do
  end
  integer function range_less(this, a, b, val) result(res)
    class(wavelet_matrix), intent(in) :: this
    integer, intent(in) :: a, b
    integer(intkind), intent(in) :: val
    logical :: bool
    integer :: l, r, i, j, d
    res = 0
    if (a > b) return
    l = a
    r = b
    do j = 0, this%m
      d = (r - l + 1) - (rank_ssa(this%bit(j), r, .true.) - rank_ssa(this%bit(j), l - 1, .true.))
      bool = btest(val, this%m - j)
      l = rank_ssa(this%bit(j), l - 1, bool) + 1
      r = rank_ssa(this%bit(j), r, bool)
      if (bool) then
        res = res + d
        l = l + this%bit(j)%offset
        r = r + this%bit(j)%offset
      end if
    end do
  end
  integer function range_freq(this, a, b, x, y) result(res)
    class(wavelet_matrix), intent(in) :: this
    integer, intent(in) :: a, b
    integer(intkind), intent(in) :: x, y
    res = 0
    if (x > y) return
    res = range_less(this, a, b, y + 1) - range_less(this, a, b, x)
  end
end module mod_wavelet_matrix
program range_set_query
  use mod_wavelet_matrix
  implicit none
  integer :: n, q, l, r, i
  integer, dimension(500000) :: c = 0, a = 0, b = 0
  type(wavelet_matrix) :: wm
  read(*,*) n, q
  read(*,*) c(1:n)
  b(1:n) = n + 1
  do i = 1, n
    a(i) = b(c(i))
    b(c(i)) = i
  end do
  wm = wavelet_matrix(a(1:n))
  do i = 1, q
    read(*,*) l, r
    write(*,'(i0)') r - l + 1 - wm%range_freq(l, r, l, r)
  end do
end program range_set_query

Measured memory_bytes: 26144