module mod_wavelet_matrix
  implicit none
  integer, private, parameter :: intkind = 4
  type wavelet_matrix
    integer :: n, m
    integer, private, allocatable :: bit(:, :), 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 :: newwm0
  end interface
  private :: find, lower_bound
contains
  type(wavelet_matrix) function newwm0(a) result(res)
    integer(intkind), intent(in) :: a(:)
    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 :: flag
    maxa = maxval(a)
    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:n, 0:m))
    z = a
    do j = 0, m
      nx = 0
      ny = 0
      do i = 1, n
        flag = btest(z(i), m - j)
        res%bit(i, j) = res%bit(i - 1, j) + merge(1, 0, flag)
        if (flag) then
          ny = ny + 1
          y(ny) = z(i)
        else
          nx = nx + 1
          x(nx) = z(i)
        end if
      end do
      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(intkind) function access_to(this, idx) result(res)
    class(wavelet_matrix), intent(in) :: this
    integer, intent(in) :: idx
    integer :: i, j
    i = idx
    res = 0
    do j = 0, this%m
      if (this%bit(i, j) - this%bit(i - 1, j) == 1) then
        res = ibset(res, this%m - j)
        i = this%bit(i, j) + this%n - this%bit(this%n, j)
      else
        i = i - this%bit(i, j)
      end if
    end do
  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 function lower_bound(n, a, b) result(res)
    integer, intent(in) :: n, a(0:n), b
    integer :: l, r, m
    l = 0
    r = n + 1
    do while (r - l > 1)
      m = (l + r) / 2
      if (a(m) < b) then
        l = m
      else
        r = m
      end if
    end do
    res = r
  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
    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
      if (btest(val, this%m - j)) then
        i = this%bit(i, j) + this%n - this%bit(this%n, j)
      else
        i = i - this%bit(i, j)
      end if
    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
    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
      if (btest(val, this%m - j)) then
        i = lower_bound(this%n, this%bit(:, j), i - (this%n - this%bit(this%n, j)))
      else
        i = lower_bound(this%n, [(l - this%bit(l, j), l = 0, this%n)], i)
      end if
    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
    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) - (this%bit(r, j) - this%bit(l - 1, j))
      if (i > d) then
        res = ibset(res, this%m - j)
        l = this%bit(l - 1, j) + this%n - this%bit(this%n, j) + 1
        r = this%bit(r, j) + this%n - this%bit(this%n, j)
        i = i - d
      else
        l = l - this%bit(l - 1, j)
        r = r - this%bit(r, j)
      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
    integer :: l, r, i, j, d
    res = 0
    if (a > b) return
    l = a
    r = b
    do j = 0, this%m
      if (btest(val, this%m - j)) then
        res = res + (r - l + 1) - (this%bit(r, j) - this%bit(l - 1, j))
        l = this%bit(l - 1, j) + this%n - this%bit(this%n, j) + 1
        r = this%bit(r, j) + this%n - this%bit(this%n, j)
      else
        l = l - this%bit(l - 1, j)
        r = r - this%bit(r, j)
      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: 59552