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
module module_modint
  use, intrinsic :: iso_fortran_env
  implicit none
  integer(int64) :: modu = 1000000007
  type modint
    integer(int64), private :: val = 0
  contains
    procedure, private :: asgn => asgn
    procedure, private :: add => add
    procedure, private :: sub => sub
    procedure, private :: mul => mul
    procedure, private :: div => div
    procedure, private :: pow => pow
    procedure, public :: inv => inv
    procedure, public :: get => get
    generic, public :: assignment(=) => asgn
    generic, public :: operator(+) => add
    generic, public :: operator(-) => sub
    generic, public :: operator(*) => mul
    generic, public :: operator(/) => div
    generic, public :: operator(**) => pow
  end type
  interface modint
    module procedure :: new_modint, new_modint_with_num
  end interface modint
  interface operator(+)
    module procedure :: add_r
  end interface
  interface operator(-)
    module procedure :: sub_r
  end interface
  interface operator(*)
    module procedure :: mul_r
  end interface
  interface operator(/)
    module procedure :: div_r
  end interface
  interface dot_product
    module procedure :: dot_prod
  end interface dot_product
  interface matmul
    module procedure :: matmul1, matmul2, matmul3
  end interface matmul
contains
  pure elemental integer(int64) function to_int64(num) result(res)
    class(*), intent(in) :: num
    select type (num)
    type is (integer(int64))
      res = num
    type is (integer(int32))
      res = int(num, int64)
    type is (integer(int16))
      res = int(num, int64)
    type is (integer(int8))
      res = int(num, int64)
    class default
      res = 0_8
      res = 1_8 / res
    end select
  end
  pure elemental type(modint) function new_modint()
  end
  pure elemental type(modint) function new_modint_with_num(num) result(res)
    class(*), intent(in) :: num
    res%val = to_int64(num)
  end
  pure elemental subroutine asgn(this, other)
    class(modint), intent(inout) :: this
    class(*), intent(in) :: other
    integer(int64) :: val
    select type (other)
    type is (modint)
      this%val = other%val
    class default
      val = to_int64(other)
      if (0_8 <= val .and. val < modu) then
        this%val = val
      else
        this%val = modulo(val, modu)
      end if
    end select
  end
  pure elemental type(modint) function add(this, other) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    integer(int64) :: val
    select type (other)
    type is (modint)
      if (0_8 <= this%val + other%val .and. this%val + other%val < modu) then
        res%val = this%val + other%val
      else
        res%val = this%val + other%val - modu
      end if
    class default
      val = to_int64(other)
      if (0_8 <= this%val + val .and. this%val + val < modu) then
        res%val = this%val + val
      else
        res%val = modulo(this%val + val, modu)
      end if
    end select
  end
  pure elemental type(modint) function sub(this, other) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    integer(int64) :: val
    select type (other)
    type is (modint)
      if (0_8 <= this%val - other%val .and. this%val - other%val < modu) then
        res%val = this%val - other%val
      else
        res%val = this%val - other%val + modu
      end if
    class default
      val = to_int64(other)
      if (0_8 <= this%val - val .and. this%val - val < modu) then
        res%val = this%val - val
      else
        res%val = modulo(this%val - val, modu)
      end if
    end select
  end
  pure elemental type(modint) function mul(this, other) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    integer(int64) :: val
    select type (other)
    type is (modint)
      val = other%val
    class default
      val = to_int64(other)
    end select
    if (0_8 <= this%val * val .and. this%val * val < modu) then
      res%val = this%val * val
    else
      res%val = modulo(this%val * val, modu)
    end if
  end
  pure elemental type(modint) function div(this,other) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    type(modint) :: tmp
    call asgn(tmp, other)
    call asgn(res, mul(this, inv(tmp)))
  end
  pure elemental integer(int64) function get(this) result(res)
    class(modint), intent(in) :: this
    res = this%val
  end
  pure elemental type(modint) function inv(this) result(res)
    class(modint), intent(in) :: this
    integer(int64) :: a, b, c, n
    integer(int64) :: x, y, z, m
    a = this%val
    b = modu
    c = 0_8
    n = 1_8
    do while (b /= 0_8)
      x = b
      y = mod(a,b)
      z = n-a/b*c
      m = c
      a = x
      b = y
      c = z
      n = m
    end do
    call asgn(res,n)
  end
  pure elemental type(modint) function pow(this, other) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    integer(int64) :: n
    type(modint) :: i
    res%val = 1
    i%val = this%val
    n = to_int64(other)
    do while (n > 0_8)
      if (btest(n, 0)) call asgn(res, mul(res, i))
      call asgn(i, mul(i, i))
      n = rshift(n, 1)
    end do
  end
  pure elemental type(modint) function add_r(other, this) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    call asgn(res, add(this, other))
  end
  pure elemental type(modint) function sub_r(other, this) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    call asgn(res, sub(this, other))
  end
  pure elemental type(modint) function mul_r(other, this) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    call asgn(res, mul(this, other))
  end
  pure elemental type(modint) function div_r(other, this) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    call asgn(res, div(this, other))
  end
  pure type(modint) function dot_prod(x, y) result(res)
    type(modint), intent(in) :: x(:), y(:)
    integer :: i
    if (size(x, 1) /= size(y, 1)) i = to_int64('')
    do i = 1, size(x,1)
      call asgn(res, add(res, mul(x(i), y(i))))
    end do
  end
  pure function matmul1(x, y) result(res)
    type(modint), intent(in) :: x(:, :), y(:)
    type(modint) :: res(size(x, 1))
    integer :: i
    do i = 1, size(x, 1)
      call asgn(res(i), dot_prod(x(i, :), y))
    end do
  end
  pure function matmul2(x, y) result(res)
    type(modint), intent(in) :: x(:), y(:, :)
    type(modint) :: res(size(y, 2))
    integer :: i
    do i = 1, size(y, 2)
      call asgn(res(i), dot_prod(x, y(:, i)))
    end do
  end
  pure function matmul3(x, y) result(res)
    type(modint), intent(in) :: x(:, :), y(:, :)
    type(modint) :: res(size(x, 1), size(y, 2))
    integer :: i
    do i = 1, size(x, 1)
      call asgn(res(i, :), matmul2(x(i, :), y))
    end do
  end
end module module_modint
module modulo_utils
  use module_modint
  implicit none
  type(modint), allocatable :: f(:), invf(:)
contains
  subroutine init(n)
    integer, intent(in) :: n
    integer :: i
    if (allocated(f)) deallocate(f)
    if (allocated(invf)) deallocate(invf)
    allocate(f(0:n), invf(0:n))
    f(0) = 1
    do i = 1, n
      f(i) = f(i - 1) * i
    end do
    invf(n) = f(n)%inv()
    do i = n, 1, -1
      invf(i - 1) = invf(i) * i
    end do
  end
  pure elemental type(modint) function perm(n, k) result(res)
    integer, intent(in) :: n, k
    if (k > n .or. n < 0 .or. k < 0) return
    res = f(n) * invf(n - k)
  end
  pure elemental type(modint) function comb(n, k) result(res)
    integer, intent(in) :: n, k
    res = perm(n, k) * invf(k)
  end
  pure elemental type(modint) function homo(n, k) result(res)
    integer, intent(in) :: n, k
    res = 1
    if (n == 0 .and. k == 0) return
    res = comb(n + k - 1, k)
  end
end module modulo_utils
program squirrel_migration
  use mod_unweighted_graph
  use module_modint
  use modulo_utils
  implicit none
  integer :: n, i, c(3), s = 0, a, b, p = 0, v, k
  integer, dimension(5000) :: x = 0, y = 0, r = 0
  type(unweighted_graph) :: g
  type(modint) :: m, dp(0:5000, 0:5000)
  read(*,*) n
  call init(n)
  do i = 1, n - 1
    read(*,*) x(i), y(i)
  end do
  g = unweighted_graph(n, n - 1, x, y)
  c = g%find_centroids()
  if (c(3) == 2) then
    m = f(n / 2) * f(n / 2)
    write(*,'(i0)') m%get()
    stop
  end if
  call dfs(g, c(1), -1)
  dp(0, 0) = 1
  k = size(g%list(c(1))%at)
  do i = 1, k
    v = g%list(c(1))%at(i)
    do a = 0, p
      do b = 0, r(v)
        dp(i, a + b) = dp(i, a + b) + dp(i - 1, a) * coef(r(v), b)
      end do
    end do
    p = p + r(v)
  end do
  do i = 0, n - 1
    m = m + f(n - i) * dp(k, i) * merge(1, -1, mod(i, 2) == 0)
  end do
  write(*,'(i0)') m%get()
contains
  recursive subroutine dfs(this, v, p)
    type(unweighted_graph), intent(in) :: this
    integer, intent(in) :: v, p
    integer :: i, u
    r(v) = 1
    do i = 1, size(this%list(v)%at)
      u = this%list(v)%at(i)
      if (u == p) cycle
      call dfs(this, u, v)
      r(v) = r(v) + r(u)
    end do
  end
  pure elemental type(modint) function coef(n, k) result(res)
    integer, intent(in) :: n, k
    res = perm(n, k) * comb(n, k)
  end
end program squirrel_migration

Measured memory_bytes: 200168