module mod_unweighted_graph
  implicit none

  type t_list
    integer, allocatable :: at(:)
  end type

  type unweighted_graph
    type(t_list), allocatable :: list(:)
  contains
    procedure :: is_bipartite => is_bipartite
    procedure :: get_cycle => get_cycle
    procedure :: topological_sort => topological_sort
    procedure :: find_centroids => find_centroids
  end type

  type unweighted_graph_list
    type(unweighted_graph), allocatable :: graphs(:)
  end type

  interface unweighted_graph
    module procedure :: newug0
  end interface unweighted_graph

contains

  type(unweighted_graph) function newug0(n, m, a, b, undirected) result(res)
    integer, intent(in) :: n, m, a(:), b(:)
    logical, optional, intent(in) :: undirected
    logical :: u
    integer :: c(n), i
    u = merge(undirected, .true., present(undirected))
    c = 0
    do i = 1, m
      c(a(i)) = c(a(i)) + 1
      if (u) c(b(i)) = c(b(i)) + 1
    end do
    allocate(res%list(n))
    do i = 1, n
      allocate(res%list(i)%at(c(i)))
    end do
    do i = m, 1, -1
      res%list(a(i))%at(c(a(i))) = b(i)
      c(a(i)) = c(a(i)) - 1
      if (u) then
        res%list(b(i))%at(c(b(i))) = a(i)
        c(b(i)) = c(b(i)) - 1
      end if
    end do
  end

  logical function is_bipartite(this) result(res)
    class(unweighted_graph), intent(in) :: this
    integer(1) :: color(size(this%list))
    integer :: i
    color = 0_1
    color(1) = 1_1
    res = ibdfs(this, color, 1)
  end

  recursive logical function ibdfs(this, color, u) result(res)
    type(unweighted_graph), intent(in) :: this
    integer(1), intent(inout) :: color(:)
    integer, intent(in) :: u
    integer :: i, v
    do i = 1, size(this%list(u)%at)
      v = this%list(u)%at(i)
      if (color(v) == color(u)) then
        res = .false.; return
      end if
      if (color(v) == 0_1) then
        color(v) = -color(u)
        if (.not.ibdfs(this, color, v)) then
          res = .false.; return
        end if
      end if
    end do
    res = .true.
  end

  function get_cycle(this) result(res)
    class(unweighted_graph), intent(in) :: this
    integer :: res(size(this%list) + 1), idx, pos, i
    logical :: used(size(this%list)), finished(size(this%list))
    used = .false.
    finished = .false.
    idx = 0
    pos = -1
    call gcdfs(this, used, finished, 1, -1, res, idx, pos)
    do i = 1, idx
      if (res(i) == pos) then
        res(1:idx - i + 1) = res(i:idx)
        idx = idx - i + 1
      end if
    end do
    res(idx + 1:size(this%list) + 1) = -1
    res(size(this%list) + 1) = idx
  end

  recursive subroutine gcdfs(this, used, finished, u, p, stack, idx, pos)
    type(unweighted_graph), intent(in) :: this
    logical, intent(inout) :: used(:), finished(:)
    integer, intent(in) :: u, p
    integer, intent(inout) :: stack(:), idx, pos
    integer :: i, v
    used(u) = .true.
    idx = idx + 1
    stack(idx) = u
    do i = 1, size(this%list(u)%at)
      v = this%list(u)%at(i)
      if (v == p) cycle
      if (finished(v)) cycle
      if (used(v)) then
        pos = v; return
      end if
      call gcdfs(this, used, finished, v, u, stack, idx, pos)
      if (pos > 0) return
    end do
    idx = idx - 1
    finished(u) = .true.
  end

  function topological_sort(this) result(res)
    class(unweighted_graph), intent(in) :: this
    integer :: res(size(this%list)), idx, i
    logical :: used(size(this%list))
    used = .false.
    idx = size(this%list)
    do i = 1, size(this%list)
      if (used(i)) cycle
      call tsdfs(this, used, i, res, idx)
    end do
  end

  recursive subroutine tsdfs(this, used, u, sorted, idx)
    type(unweighted_graph), intent(in) :: this
    logical, intent(inout) :: used(:)
    integer, intent(in) :: u
    integer, intent(inout) :: sorted(:), idx
    integer :: i, v
    used(u) = .true.
    do i = 1, size(this%list(u)%at)
      v = this%list(u)%at(i)
      if (used(v)) cycle
      call tsdfs(this, used, v, sorted, idx)
    end do
    sorted(idx) = u
    idx = idx - 1
  end

  function find_centroids(this) result(res)
    class(unweighted_graph), intent(in) :: this
    integer :: sizes(size(this%list))
    integer :: res(3), idx, i
    res = -1
    idx = 0
    call fcdfs(this, sizes, 1, -1, res, idx)
    res(3) = idx
  end

  recursive subroutine fcdfs(this, sizes, u, p, centroids, idx)
    type(unweighted_graph), intent(in) :: this
    integer, intent(in) :: u, p
    integer, intent(inout) :: sizes(:), centroids(:), idx
    integer :: i, v
    logical :: found
    found = .true.
    sizes(u) = 1
    do i = 1, size(this%list(u)%at)
      v = this%list(u)%at(i)
      if (v == p) cycle
      call fcdfs(this, sizes, v, u, centroids, idx)
      if (sizes(v) > size(this%list) / 2) found = .false.
      sizes(u) = sizes(u) + sizes(v)
    end do
    if (size(this%list) - sizes(u) > size(this%list) / 2) found = .false.
    if (found) then
      idx = idx + 1
      centroids(idx) = u
    end if
  end

  type(unweighted_graph_list) function strongly_connected_components(n, m, a, b) result(res)
    integer, intent(in) :: n, m, a(:), b(:)
    type(unweighted_graph) :: g
    integer :: order(n), number(n), idx, i, j
    logical :: used(n)
    integer :: x(m, 2), y(m, 2), cnt(2)
    g = newug0(n, m, a, b, .false.)
    used = .false.
    idx = n
    do i = 1, n
      if (used(i)) cycle
      call sccdfs1(g, used, i, order, idx)
    end do
    g = newug0(n, m, b, a, .false.)
    number = -1
    idx = 0
    do i = 1, n
      if (number(order(i)) > 0) cycle
      idx = idx + 1
      call sccdfs2(g, order(i), number, idx)
    end do
    cnt = 0
    do i = 1, m
      j = merge(2, 1, number(a(i)) == number(b(i)))
      cnt(j) = cnt(j) + 1
      x(cnt(j), j) = a(i)
      y(cnt(j), j) = b(i)
    end do
    call remove_multiedge(cnt(1), x(:, 1), y(:, 1))
    allocate(res%graphs(3))
    do i = 1, 2
      res%graphs(i) = newug0(n, cnt(i), x(:, i), y(:, i), .false.)
    end do
    allocate(res%graphs(3)%list(1))
    allocate(res%graphs(3)%list(1)%at(n))
    res%graphs(3)%list(1)%at = number
  end

  recursive subroutine sccdfs1(this, used, u, order, idx)
    type(unweighted_graph), intent(in) :: this
    logical, intent(inout) :: used(:)
    integer, intent(in) :: u
    integer, intent(inout) :: order(:), idx
    integer :: i, v
    used(u) = .true.
    do i = 1, size(this%list(u)%at)
      v = this%list(u)%at(i)
      if (used(v)) cycle
      call sccdfs1(this, used, v, order, idx)
    end do
    idx = idx - 1
    order(idx) = u
  end

  recursive subroutine sccdfs2(this, u, number, idx)
    type(unweighted_graph), intent(in) :: this
    integer, intent(in) :: u, idx
    integer, intent(inout) :: number(:)
    integer :: i, v
    number(u) = idx
    do i = 1, size(this%list(u)%at)
      v = this%list(u)%at(i)
      if (number(v) > 0) cycle
      call sccdfs2(this, v, number, idx)
    end do
  end

  subroutine remove_multiedge(m, a, b)
    integer, intent(inout) :: m, a(:), b(:)
    integer :: i, x(2, size(a)), idx
    do i = 1, m
      x(1, i) = a(i)
      x(2, i) = b(i)
    end do
    call merge_sort(x(1:2, 1:m))
    a(1) = x(1, 1)
    b(1) = x(2, 1)
    idx = 1
    do i = 2, m
      if (a(idx) == x(1, i) .and. b(idx) == x(2, i)) cycle
      idx = idx + 1
      a(idx) = x(1, i)
      b(idx) = x(2, i)
    end do
    m = idx
  end

  subroutine merge_sort(a)
    integer, intent(inout) :: a(:, :)
    integer :: n, m, l, i, u, j
    n = size(a,1)
    m = n
    l = 1
    do while (m > 1)
      do i = 1, m / 2
        u = min(2 * i * l, n)
        j = 2 * (i - 1) * l + 1
        call merger(a(j:j + l - 1, :), a(j + l:u, :))
      end do
      l = 2 * l
      m = (m + 1) / 2
    end do
  end

  subroutine merger(a1, a2)
    integer, intent(inout) :: a1(:, :), a2(:, :)
    integer :: a(size(a1, 1) + size(a2, 1), size(a1, 2))
    integer :: i1, i2, n1, n2
    i1 = 1
    i2 = 1
    n1 = size(a1, 1)
    n2 = size(a2, 1)
    do while (i1 <= n1 .and. i2 <= n2)
      if (a1(i1, 1) < a2(i2, 1) .or. &
        (a1(i1, 1) == a2(i2, 1) .and. a1(i1, 2) <= a2(i2, 2))) then
        a(i1 + i2 - 1, :) = a1(i1, :)
        i1 = i1 + 1
      else
        a(i1 + i2 - 1, :) = a2(i2, :)
        i2 = i2 + 1
       end if
    end do
    if (i1 <= n1) then
      a(i1 + i2 - 1:n1 + n2, :) = a1(i1:n1, :)
    else
      a(i1 + i2 - 1:n1 + n2, :) = a2(i2:n2, :)
    end if
    a1 = a(1:n1, :)
    a2 = a(n1 + 1:n1 + n2, :)
  end

end module mod_unweighted_graph
program urban_planning
  use mod_unweighted_graph
  implicit none
  integer(8), parameter :: modu = 1000000007
  integer :: n, m = 0, x, y, acc = 0, cnt = 0, i, j, k, l
  integer, dimension(5000) :: p = 0, a = 0, b = 0, c = 0
  integer(8), dimension(0:5000) :: fac = 1, pow = 1
  integer(8) :: dp(0:2, 0:5000) = 0, ans
  logical :: used(5000) = .false.
  type(unweighted_graph) :: g
  read(*,*) n
  read(*,*) p(1:n)
  do i = 1, n
    fac(i) = mod(fac(i - 1) * i, modu)
    pow(i) = mod(pow(i - 1) * (n - 1), modu)
    if (p(i) < 0) cycle
    m = m + 1
    a(m) = i
    b(m) = p(i)
  end do
  g = unweighted_graph(n, m, a, b)
  do i = 1, n
    if (used(i)) cycle
    x = 0
    y = 0
    call dfs(g, i)
    y = y / 2
    if (y == x - 1) then
      acc = acc + 1
      c(acc) = x
    else
      cnt = cnt + 1
    end if
  end do
  dp(0, 0) = 1
  do i = 1, acc
    do j = 0, i
      k = and(i, 1)
      l = 1 - mod(i, 2)
      dp(k, j) = mod(dp(l, j) * (n - 1) + dp(l, j - 1) * c(i), modu)
    end do
  end do
  ans = mod(pow(acc) * cnt, modu)
  k = mod(acc, 2)
  do i = 2, n
    ans = mod(ans + dp(k, i) * fac(i - 1), modu)
  end do
  do i = 1, acc
    ans = mod(ans + (c(i) - 1) * pow(acc - 1), modu)
  end do
  write(*,'(i0)') modulo(pow(acc) * n - ans, modu)
contains
  recursive subroutine dfs(this, v)
    type(unweighted_graph), intent(in) :: this
    integer, intent(in) :: v
    integer :: i, u
    used(v) = .true.
    x = x + 1
    do i = 1, size(this%list(v)%at)
      y = y + 1
      u = this%list(v)%at(i)
      if (used(u)) cycle
      call dfs(this, u)
    end do
  end
end program urban_planning

Measured memory_bytes: 3608