Problem:


Score : 700 points

Problem Statement

There are N towns numbered 1, 2, \cdots, N.

Some roads are planned to be built so that each of them connects two distinct towns bidirectionally. Currently, there are no roads connecting towns.

In the planning of construction, each town chooses one town different from itself and requests the following: roads are built so that the chosen town is reachable from itself using one or more roads.

These requests from the towns are represented by an array P_1, P_2, \cdots, P_N. If P_i = -1, it means that Town i has not chosen the request; if 1 \leq P_i \leq N, it means that Town i has chosen Town P_i.

Let K be the number of towns i such that P_i = -1. There are (N-1)^K ways in which the towns can make the requests. For each way to make requests, find the minimum number of roads needed to meet all the requests, and print the sum of those (N-1)^K numbers, modulo (10^9+7).

Constraints

2 \leq N \leq 5000

P_i = -1 or 1 \leq P_i \leq N.

P_i \neq i

All values in input are integers.

Input

Input is given from Standard Input in the following format:

N
P_1 P_2 \cdots P_N

Output

For each way to make requests, find the minimum number of roads needed to meet all the requests, and print the sum of those (N-1)^K numbers, modulo (10^9+7).

Sample Input 1

4
2 1 -1 3

Sample Output 1

8

There are three ways to make requests, as follows:

Choose P_1 = 2, P_2 = 1, P_3 = 1, P_4 = 3. In this case, at least three roads - for example, (1,2),(1,3),(3,4) - are needed to meet the requests.

Choose P_1 = 2, P_2 = 1, P_3 = 2, P_4 = 3. In this case, at least three roads - for example, (1,2),(1,3),(3,4) - are needed to meet the requests.

Choose P_1 = 2, P_2 = 1, P_3 = 4, P_4 = 3. In this case, at least two roads - for example, (1,2),(3,4) - are needed to meet the requests.

Note that it is not mandatory to connect Town i and Town P_i directly.

The sum of the above numbers is 8.

Sample Input 2

2
2 1

Sample Output 2

1

There may be just one fixed way to make requests.

Sample Input 3

10
2 6 9 -1 6 9 -1 -1 -1 -1

Sample Output 3

527841



Sample input:
4
2 1 -1 3


Sample output:
8

Submission status: Wrong Answer

Solution code:
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor, cos, radians, pi, sin
from operator import mul
from functools import reduce
from operator import mul

mod = 10 ** 9 + 7
sys.setrecursionlimit(2147483647)
INF = 10 ** 13
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]


class UnionFind:
    def __init__(self, n):
        # 負  : 根であることを示す。絶対値はランクを示す
        # 非負: 根でないことを示す。値は親を示す
        self.table = [-1] * n
        self.size = [1] * n
        self.group_num = n

    def root(self, x):
        if self.table[x] < 0:
            return x
        else:
            self.table[x] = self.root(self.table[x])
            return self.table[x]

    def get_size(self, x):
        r = self.root(x)
        return self.size[r]

    def is_same(self, x, y):
        return self.root(x) == self.root(y)

    def union(self, x, y):
        r1 = self.root(x)
        r2 = self.root(y)
        if r1 == r2:
            return
        # ランクの取得
        d1 = self.table[r1]
        d2 = self.table[r2]
        if d1 <= d2:
            self.table[r2] = r1
            self.size[r1] += self.size[r2]
            if d1 == d2:
                self.table[r1] -= 1
        else:
            self.table[r1] = r2
            self.size[r2] += self.size[r1]
        self.group_num -= 1

n = I()
P = LI()

fac = [1] * (n + 1)
inv = [1] * (n + 1)
for j in range(1, n + 1):
    fac[j] = fac[j-1] * j % mod


inv[n] = pow(fac[n], mod-2, mod)
for j in range(n-1, -1, -1):
    inv[j] = inv[j+1] * (j+1) % mod


def comb(n, r):
    if r > n or n < 0 or r < 0:
        return 0
    return fac[n] * inv[n - r] * inv[r] % mod


U = UnionFind(n)
for i in range(n):
    if P[i] != -1:
        U.union(P[i] - 1, i)

k = 0
L = []
for j in range(n):
    if P[j] == -1:
        L += [U.size[j]]
        k += 1

dp = [0]*(k+1)
dp[0] = 1
dp_new = dp[:]
for i in range(1, k+1):
    for j in range(i):
        dp_new[j+1] += dp[j]*L[i-1]
    dp = dp_new[:]


if k:
    dp[1] = sum(L) - len(L)

cycle_count = 0
for i in range(1, k + 1):
    # i個のグループからなるサイクルがひとつできる
    buf = dp[i] * factorial(i - 1)
    buf %= mod
    buf *= pow(n - 1, k - i, mod)
    buf %= mod
    cycle_count += buf
    # print("数え上げ", dp[i], factorial(i-1), pow(n-1, k-i, mod))
# print(f"cycle_count, {cycle_count}")

# すでに固定サイクルがあれば数える
cycle_count += (U.group_num - k) * pow(n - 1, k, mod)

# \sigma_{つなぎ方} n
ans = n * pow(n - 1, k, mod)
# print("sigma_{つなぎ方} n", ans)
ans %= mod
ans = ans - cycle_count
ans %= mod
print(ans)

Submission status: Wrong Answer

Solution code:
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor, cos, radians, pi, sin
from operator import mul
from functools import reduce
from operator import mul

mod = 10 ** 9 + 7
sys.setrecursionlimit(2147483647)
INF = 10 ** 13
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]


class UnionFind:
    def __init__(self, n):
        # 負  : 根であることを示す。絶対値はランクを示す
        # 非負: 根でないことを示す。値は親を示す
        self.table = [-1] * n
        self.size = [1] * n
        self.group_num = n

    def root(self, x):
        if self.table[x] < 0:
            return x
        else:
            self.table[x] = self.root(self.table[x])
            return self.table[x]

    def get_size(self, x):
        r = self.root(x)
        return self.size[r]

    def is_same(self, x, y):
        return self.root(x) == self.root(y)

    def union(self, x, y):
        r1 = self.root(x)
        r2 = self.root(y)
        if r1 == r2:
            return
        # ランクの取得
        d1 = self.table[r1]
        d2 = self.table[r2]
        if d1 <= d2:
            self.table[r2] = r1
            self.size[r1] += self.size[r2]
            if d1 == d2:
                self.table[r1] -= 1
        else:
            self.table[r1] = r2
            self.size[r2] += self.size[r1]
        self.group_num -= 1

n = I()
P = LI()

fac = [1] * (n + 1)
inv = [1] * (n + 1)
for j in range(1, n + 1):
    fac[j] = fac[j-1] * j % mod


inv[n] = pow(fac[n], mod-2, mod)
for j in range(n-1, -1, -1):
    inv[j] = inv[j+1] * (j+1) % mod


def comb(n, r):
    if r > n or n < 0 or r < 0:
        return 0
    return fac[n] * inv[n - r] * inv[r] % mod


U = UnionFind(n)
for i in range(n):
    if P[i] != -1:
        U.union(P[i] - 1, i)

k = 0
L = []
for j in range(n):
    if P[j] == -1:
        L += [U.size[j]]
        k += 1

dp = [0]*(k+1)
dp[0] = 1
dp_new = dp[:]
for i in range(1, k+1):
    for j in range(i):
        dp_new[j+1] += dp[j]*L[i-1]
    dp = dp_new[:]


if k:
    tot = 0
    for s in L:
        tot += s-1
    dp[1] = tot

cycle_cnt = 0
for i in range(1, len(L) + 1):
    cycle_cnt = (cycle_cnt + fac[i - 1] * dp[i] * pow(n - 1, len(L) - i, mod)) % mod


cycle_cnt = cycle_cnt + (U.group_num - k) * pow(n - 1, k, mod) % mod
ans = n * pow(n - 1, k, mod) % mod
print((ans - cycle_cnt) % mod)

Submission status: Wrong Answer

Solution code:
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor, cos, radians, pi, sin
from operator import mul
from functools import reduce
from operator import mul

mod = 10 ** 9 + 7
sys.setrecursionlimit(2147483647)
INF = 10 ** 13
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]


class UnionFind:
    def __init__(self, n):
        # 負  : 根であることを示す。絶対値はランクを示す
        # 非負: 根でないことを示す。値は親を示す
        self.table = [-1] * n
        self.size = [1] * n
        self.group_num = n

    def root(self, x):
        if self.table[x] < 0:
            return x
        else:
            self.table[x] = self.root(self.table[x])
            return self.table[x]

    def get_size(self, x):
        r = self.root(x)
        return self.size[r]

    def is_same(self, x, y):
        return self.root(x) == self.root(y)

    def union(self, x, y):
        r1 = self.root(x)
        r2 = self.root(y)
        if r1 == r2:
            return
        # ランクの取得
        d1 = self.table[r1]
        d2 = self.table[r2]
        if d1 <= d2:
            self.table[r2] = r1
            self.size[r1] += self.size[r2]
            if d1 == d2:
                self.table[r1] -= 1
        else:
            self.table[r1] = r2
            self.size[r2] += self.size[r1]
        self.group_num -= 1

n = I()
P = LI()

fac = [1] * (n + 1)
inv = [1] * (n + 1)
for j in range(1, n + 1):
    fac[j] = fac[j-1] * j % mod


inv[n] = pow(fac[n], mod-2, mod)
for j in range(n-1, -1, -1):
    inv[j] = inv[j+1] * (j+1) % mod


def comb(n, r):
    if r > n or n < 0 or r < 0:
        return 0
    return fac[n] * inv[n - r] * inv[r] % mod


U = UnionFind(n)
for i in range(n):
    if P[i] != -1:
        U.union(P[i] - 1, i)

k = 0
L = []
for j in range(n):
    if P[j] == -1:
        L += [U.size[j]]
        k += 1

dp = [0]*(k+1)
dp[0] = 1
dp_new = dp[:]
for i in range(1, k+1):
    for j in range(i):
        dp_new[j+1] += dp[j]*L[i-1]
    dp = dp_new[:]

if L:
    dp[1] = sum(L) - len(L)
cycle_cnt = 0
for i in range(1, len(L) + 1):
    cycle_cnt = (cycle_cnt + fac[i - 1] * dp[i] * pow(n - 1, len(L) - i, mod)) % mod


cycle_cnt = cycle_cnt + (U.group_num - k) * pow(n - 1, k, mod) % mod
ans = n * pow(n - 1, k, mod) % mod
print((ans - cycle_cnt) % mod)

Submission status: Wrong Answer

Solution code:
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor, cos, radians, pi, sin
from operator import mul
from functools import reduce
from operator import mul

mod = 10 ** 9 + 7
sys.setrecursionlimit(2147483647)
INF = 10 ** 13
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]


class UnionFind:
    def __init__(self, n):
        # 負  : 根であることを示す。絶対値はランクを示す
        # 非負: 根でないことを示す。値は親を示す
        self.table = [-1] * n
        self.size = [1] * n
        self.group_num = n

    def root(self, x):
        if self.table[x] < 0:
            return x
        else:
            self.table[x] = self.root(self.table[x])
            return self.table[x]

    def get_size(self, x):
        r = self.root(x)
        return self.size[r]

    def is_same(self, x, y):
        return self.root(x) == self.root(y)

    def union(self, x, y):
        r1 = self.root(x)
        r2 = self.root(y)
        if r1 == r2:
            return
        # ランクの取得
        d1 = self.table[r1]
        d2 = self.table[r2]
        if d1 <= d2:
            self.table[r2] = r1
            self.size[r1] += self.size[r2]
            if d1 == d2:
                self.table[r1] -= 1
        else:
            self.table[r1] = r2
            self.size[r2] += self.size[r1]
        self.group_num -= 1

n = I()
P = LI()

fac = [1] * (n + 1)
inv = [1] * (n + 1)
for j in range(1, n + 1):
    fac[j] = fac[j-1] * j % mod


inv[n] = pow(fac[n], mod-2, mod)
for j in range(n-1, -1, -1):
    inv[j] = inv[j+1] * (j+1) % mod


def comb(n, r):
    if r > n or n < 0 or r < 0:
        return 0
    return fac[n] * inv[n - r] * inv[r] % mod


U = UnionFind(n)
for i in range(n):
    if P[i] != -1:
        U.union(P[i] - 1, i)

ret = 0
L = []
for j in range(n):
    if P[j] == -1:
        L += [U.size[j]]
        ret += 1

dp = [0] * (len(L) + 1)
dp[0] = 1
for l in range(1, len(L) + 1):
    for m in range(l, 0, -1):
        dp[m] += dp[m - 1] * L[l - 1]

if L:
    dp[1] = sum(L) - len(L)
cycle_cnt = 0
for i in range(1, len(L) + 1):
    cycle_cnt = (cycle_cnt + fac[i - 1] * dp[i] * pow(n - 1, len(L) - i, mod)) % mod


cycle_cnt = cycle_cnt + (U.group_num - ret) * pow(n - 1, ret, mod) % mod
ans = n * pow(n - 1, ret, mod) % mod
print((ans - cycle_cnt) % mod)

Submission status: Wrong Answer

Solution code:
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor, cos, radians, pi, sin
from operator import mul
from functools import reduce
from operator import mul

mod = 10 ** 9 + 7
sys.setrecursionlimit(2147483647)
INF = 10 ** 13
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]


class UnionFind:
    def __init__(self, n):
        # 負  : 根であることを示す。絶対値はランクを示す
        # 非負: 根でないことを示す。値は親を示す
        self.table = [-1] * n
        self.size = [1] * n
        self.group_num = n

    def root(self, x):
        if self.table[x] < 0:
            return x
        else:
            self.table[x] = self.root(self.table[x])
            return self.table[x]

    def get_size(self, x):
        r = self.root(x)
        return self.size[r]

    def is_same(self, x, y):
        return self.root(x) == self.root(y)

    def union(self, x, y):
        r1 = self.root(x)
        r2 = self.root(y)
        if r1 == r2:
            return
        # ランクの取得
        d1 = self.table[r1]
        d2 = self.table[r2]
        if d1 <= d2:
            self.table[r2] = r1
            self.size[r1] += self.size[r2]
            if d1 == d2:
                self.table[r1] -= 1
        else:
            self.table[r1] = r2
            self.size[r2] += self.size[r1]
        self.group_num -= 1

n = I()
P = LI()

fac = [1] * (n + 1)
inv = [1] * (n + 1)
for j in range(1, n + 1):
    fac[j] = fac[j-1] * j % mod


inv[n] = pow(fac[n], mod-2, mod)
for j in range(n-1, -1, -1):
    inv[j] = inv[j+1] * (j+1) % mod


def comb(n, r):
    if r > n or n < 0 or r < 0:
        return 0
    return fac[n] * inv[n - r] * inv[r] % mod


U = UnionFind(n)
for i in range(n):
    if P[i] != -1:
        U.union(P[i] - 1, i)

ret = 0
L = []
for j in range(n):
    if P[j] == -1:
        L += [U.size[j]]
        ret += 1

dp = [0] * (len(L) + 1)
dp[0] = 1
for l in range(1, len(L) + 1):
    for m in range(l, 0, -1):
        dp[m] += dp[m - 1] * L[l - 1]

if L:
    dp[1] = sum(L) - len(L)
cycle_cnt = 0
for i in range(1, len(L) + 1):
    cycle_cnt = (cycle_cnt + fac[i - 1] * dp[i] * pow(n - 1, len(L) - i, mod)) % mod

ans = (n + ret - U.group_num) * pow(n - 1, ret, mod) % mod
print((ans - cycle_cnt) % mod)

Submission status: Wrong Answer

Solution code:
import sys
import math
from collections import defaultdict

sys.setrecursionlimit(10**7)
def input():
    return sys.stdin.readline()[:-1]

mod = 10**9 + 7

def I(): return int(input())
def LI(): return list(map(int, input().split()))
def LIR(row,col):
    if row <= 0:
        return [[] for _ in range(col)]
    elif col == 1:
        return [I() for _ in range(row)]
    else:
        read_all = [LI() for _ in range(row)]
        return map(list, zip(*read_all))

#################

class UF():
    def __init__(self,n):
        self.par = [i for i in range(n)]
        self.rank = [0]*n
        self.size = [1]*n

    # x が属する集合の代表値を返す
    def find(self,x):
        if self.par[x] == x:
            return x
        else:
            self.par[x] = y = self.find(self.par[x])
            return y

    # x と y の集合を結合する
    def union(self,x,y):
        px = self.find(x)
        py = self.find(y)
        if px == py:
            return
        else:
            if self.rank[px] < self.rank[py]:
                px, py = py, px
            if self.rank[px] == self.rank[py]:
                self.rank[px] += 1
            self.par[py] = px
            self.size[px] += self.size[py]
            return

    # x と y が同じ集合に属するか判定する
    def is_same(self,x,y):
        return self.find(x) == self.find(y)

    # x が属する集合の要素数を返す
    def get_size(self,x):
        return self.size[self.find(x)]

N = I()
P = LI()

uf = UF(N)
for i in range(N):
    if P[i] != -1:
        uf.union(i,P[i]-1)

n = []
for i in range(N):
    if P[i] == -1:
        n.append(uf.get_size(i))

if not n:
    print(1)
    exit()

K = len(n)

# beki[i] = i!
beki = [1]*K
for i in range(1,K):
    beki[i] = beki[i-1]*i
    beki[i] %= mod

# fact[i] = (N-1)^i
fact = [1]*(K+1)
for i in range(1,K+1):
    fact[i] = fact[i-1]*(N-1)
    fact[i] %= mod

# -1を含むグループのみを取り出して考える
# dp[i][j]: 0~i番目のグループのうち，計jグループから1個ずつ取り出す場合の数
dp = [[0]*(K+1) for _ in range(K)]
for i in range(K):
    for j in range(K+1):
        if j == 0:
            dp[i][0] = 1
        elif i == 0:
            if j == 1:
                dp[0][1] = n[0]
        elif i+1 < j:
            continue
        else:
            dp[i][j] = dp[i-1][j] + n[i]*dp[i-1][j-1]
            dp[i][j] %= mod

# roop[i][j]: 0~i番目のグループのうち，jグループを使ったループの場合の数
roop = [[0]*(K+1) for _ in range(K)]
for i in range(K):
    for j in range(K+1):
        if j == 0:
            roop[i][0] = 1
        elif i == 0:
            if j == 1:
                roop[0][1] = n[0]-1
        elif i+1 < j:
            continue
        else:
            if j == 1:
                roop[i][j] = roop[i-1][j] + beki[j-1]*(n[i]-1)*dp[i-1][j-1]
            else:
                roop[i][j] = roop[i-1][j] + beki[j-1]*n[i]*dp[i-1][j-1]
            roop[i][j] %= mod

# cases[j]: Kグループのうち，jグループを使ったループの場合の数
# ループに使わない要素の行き先（N-1通り）も考慮する
cases = [0]*(K+1)
for j in range(K+1):
    cases[j] = roop[K-1][j]*fact[K-j]
    cases[j] %= mod

num = pow(N-1,K,mod)
default = 0
for i in range(N):
    if i == uf.find(i):
        default += uf.get_size(i)-1
ans = num*(default+K)
ans %= mod

for j in range(1,K+1):
    ans -= cases[j]
    ans %= mod

print(ans)

Submission status: Wrong Answer

Solution code:
import sys
INF = 1 << 60
MOD = 10**9 + 7 # 998244353
sys.setrecursionlimit(2147483647)
input = lambda:sys.stdin.readline().rstrip()

class modfact(object):
    def __init__(self, n):
        fact, invfact = [1] * (n + 1), [1] * (n + 1)
        for i in range(1, n + 1): fact[i] = i * fact[i - 1] % MOD
        invfact[n] = pow(fact[n], MOD - 2, MOD)
        for i in range(n - 1, -1, -1): invfact[i] = invfact[i + 1] * (i + 1) % MOD
        self._fact, self._invfact = fact, invfact

    def inv(self, n):
        return self._fact[n - 1] * self._invfact[n] % MOD

    def fact(self, n):
        return self._fact[n]

    def invfact(self, n):
        return self._invfact[n]

    def comb(self, n, k):
        if k < 0 or n < k: return 0
        return self._fact[n] * self._invfact[k] % MOD * self._invfact[n - k] % MOD

    def perm(self, n, k):
        if k < 0 or n < k: return 0
        return self._fact[n] * self._invfact[n - k] % MOD

class UnionFind(object):
    def __init__(self, n, recursion = False):
        self._par = list(range(n))
        self._size = [1] * n
        self._recursion = recursion

    def root(self, k):
        if self._recursion:
            if k == self._par[k]:
                return k
            self._par[k] = self.root(self._par[k])
            return self._par[k]
        else:
            root = k
            while root != self._par[root]: root = self._par[root]
            while k != root: k, self._par[k] = self._par[k], root
            return root

    def unite(self, i, j):
        i, j = self.root(i), self.root(j)
        if i == j: return False
        if self._size[i] < self._size[j]: i, j = j, i
        self._par[j] = i
        self._size[i] += self._size[j]
        return True

    def is_connected(self, i, j):
        return self.root(i) == self.root(j)

    def size(self, k):
        return self._size[self.root(k)]

def resolve():
    n = int(input())
    P = [a - 1 if a != -1 else -1 for a in map(int, input().split())]
    k = P.count(-1)

    if k == 0:
        print(1)
        return

    E = [[] for _ in range(n)]
    for u, v in enumerate(P):
        if v != -1:
            E[u].append(v)
            E[v].append(u)

    # divide E into connected component
    color = [-1] * n
    cnt = 0
    def dfs(v):
        if color[v] != -1:
            return False
        color[v] = cnt
        stack = [v]
        while stack:
            v = stack.pop()
            for nv in E[v]:
                if color[nv] == -1:
                    color[nv] = cnt
                    stack.append(nv)
        return True
    for v in range(n):
        cnt += dfs(v)

    # detect loops for each color with union find
    has_loop = [False] * cnt
    uf = UnionFind(n)
    for v in range(n):
        nv = P[v]
        if nv == -1:
            continue
        if not uf.unite(v, nv):
            has_loop[color[v]] = True

    S = []
    for root in set(uf.root(v) for v in range(n)):
        if not has_loop[color[root]]:
            S.append(uf.size(root))

    ans = pow(n - 1, k, MOD) * (n - sum(has_loop))
    l = len(S)
    dp = [0] * (l + 1)
    dp[0] = 1
    for s in S:
        ndp = dp[:]
        for i in range(1, l + 1):
            ndp[i] += s * dp[i - 1]
            ndp[i] %= MOD
        dp = ndp

    # i = 1
    ans -= sum(s - 1 for s in S) % MOD * pow(n - 1, k - 1, MOD)
    ans %= MOD

    # i >= 2
    mf = modfact(5001)
    for i in range(2, l + 1):
        ans -= dp[i] * mf.fact(i - 1) % MOD * pow(n - 1, k - i, MOD)
        ans %= MOD

    print(ans)
resolve()

Submission status: Wrong Answer

Solution code:
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor, cos, radians, pi, sin
from operator import mul
from functools import reduce
from operator import mul

mod = 10 ** 9 + 7
sys.setrecursionlimit(2147483647)
INF = 10 ** 13
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]


class UnionFind:
    def __init__(self, n):
        # 負  : 根であることを示す。絶対値はランクを示す
        # 非負: 根でないことを示す。値は親を示す
        self.table = [-1] * n
        self.size = [1] * n
        self.group_num = n

    def root(self, x):
        if self.table[x] < 0:
            return x
        else:
            self.table[x] = self.root(self.table[x])
            return self.table[x]

    def get_size(self, x):
        r = self.root(x)
        return self.size[r]

    def is_same(self, x, y):
        return self.root(x) == self.root(y)

    def union(self, x, y):
        r1 = self.root(x)
        r2 = self.root(y)
        if r1 == r2:
            return
        # ランクの取得
        d1 = self.table[r1]
        d2 = self.table[r2]
        if d1 <= d2:
            self.table[r2] = r1
            self.size[r1] += self.size[r2]
            if d1 == d2:
                self.table[r1] -= 1
        else:
            self.table[r1] = r2
            self.size[r2] += self.size[r1]
        self.group_num -= 1

n = I()
P = LI()

fac = [1] * (n + 1)
inv = [1] * (n + 1)
for j in range(1, n + 1):
    fac[j] = fac[j-1] * j % mod


inv[n] = pow(fac[n], mod-2, mod)
for j in range(n-1, -1, -1):
    inv[j] = inv[j+1] * (j+1) % mod


def comb(n, r):
    if r > n or n < 0 or r < 0:
        return 0
    return fac[n] * inv[n - r] * inv[r] % mod


ret = 0
for p in P:
    if p == -1:
        ret += 1

U = UnionFind(n)
for i in range(n):
    if P[i] != -1:
        U.union(P[i] - 1, i)

L = []
for j in range(n):
    if P[j] == -1:
        L += [U.size[j]]

dp = [0] * (len(L) + 1)
dp[0] = 1
for l in range(1, len(L) + 1):
    for m in range(l, 0, -1):
        dp[m] += dp[m - 1] * L[l - 1]

if L:
    dp[1] = sum(L) - len(L)
cycle_cnt = 0
for i in range(1, len(L) + 1):
    cycle_cnt = (cycle_cnt + fac[i - 1] * dp[i] * pow(n - 1, len(L) - i, mod)) % mod

ans = (n - (U.group_num - ret)) * pow(n - 1, ret, mod) % mod
print((ans - cycle_cnt) % mod)

Submission status: Wrong Answer

Solution code:
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor, cos, radians, pi, sin
from operator import mul
from functools import reduce
from operator import mul

mod = 10 ** 9 + 7
sys.setrecursionlimit(2147483647)
INF = 10 ** 13
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]


class UnionFind:
    def __init__(self, n):
        # 負  : 根であることを示す。絶対値はランクを示す
        # 非負: 根でないことを示す。値は親を示す
        self.table = [-1] * n
        self.size = [1] * n
        self.group_num = n

    def root(self, x):
        if self.table[x] < 0:
            return x
        else:
            self.table[x] = self.root(self.table[x])
            return self.table[x]

    def get_size(self, x):
        r = self.root(x)
        return self.size[r]

    def is_same(self, x, y):
        return self.root(x) == self.root(y)

    def union(self, x, y):
        r1 = self.root(x)
        r2 = self.root(y)
        if r1 == r2:
            return
        # ランクの取得
        d1 = self.table[r1]
        d2 = self.table[r2]
        if d1 <= d2:
            self.table[r2] = r1
            self.size[r1] += self.size[r2]
            if d1 == d2:
                self.table[r1] -= 1
        else:
            self.table[r1] = r2
            self.size[r2] += self.size[r1]
        self.group_num -= 1

n = I()
P = LI()

fac = [1] * (n + 1)
inv = [1] * (n + 1)
for j in range(1, n + 1):
    fac[j] = fac[j-1] * j % mod


inv[n] = pow(fac[n], mod-2, mod)
for j in range(n-1, -1, -1):
    inv[j] = inv[j+1] * (j+1) % mod


def comb(n, r):
    if r > n or n < 0 or r < 0:
        return 0
    return fac[n] * inv[n - r] * inv[r] % mod


ret = 0
for p in P:
    if p == -1:
        ret += 1

U = UnionFind(n)
for i in range(n):
    if P[i] != -1:
        U.union(P[i] - 1, i)



L = []
for j in range(n):
    if P[j] == -1:
        L += [U.size[j]]

dp = [0] * (len(L) + 1)
dp[0] = 1
for l in range(1, len(L) + 1):
    for m in range(l, -1, -1):
        dp[m] += dp[m - 1] * L[l - 1]

if L:
    dp[1] = sum(L) - len(L)
cycle_cnt = 0
for i in range(1, len(L) + 1):
    cycle_cnt = (cycle_cnt + fac[i - 1] * dp[i] * pow(n - 1, len(L) - i, mod)) % mod

ans = (n - (U.group_num - ret)) * pow(n - 1, ret, mod) % mod
print((ans - cycle_cnt) % mod)