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

Solution (Python):
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
from functools import lru_cache

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)]


@lru_cache(maxsize=None)
def factorial(n):
    if n == 0:
        return 1
    else:
        return (n*factorial(n-1)) % mod


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.get_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):
    
    buf = dp[i] * factorial(i - 1)
    buf %= mod
    buf *= pow(n - 1, k - i, mod)
    buf %= mod
    cycle_count += buf



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

Solution (Python):
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]
        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)

Solution (Python):
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
        
        s1 = self.table[r1]
        s2 = self.table[r2]
        if s2 <= s1:
            self.table[r2] = r1
            self.size[r1] += self.size[r2]
        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)

Solution (Python):
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.n = n
        self.parents = [-1] * n

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

    def union(self, x, y):
        x = self.find(x)
        y = self.find(y)

        if x == y:
            return

        if self.parents[x] > self.parents[y]:
            x, y = y, x

        self.parents[x] += self.parents[y]
        self.parents[y] = x

    def size(self, x):
        return -self.parents[self.find(x)]

    def same(self, x, y):
        return self.find(x) == self.find(y)

    def members(self, x):
        root = self.find(x)
        return [i for i in range(self.n) if self.find(i) == root]

    def roots(self):
        return [i for i, x in enumerate(self.parents) if x < 0]

    def group_count(self):
        return len(self.roots())

    def all_group_members(self):
        return {r: self.members(r) for r in self.roots()}


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_count() - ret)) * pow(n - 1, ret, mod) % mod
print((ans - cycle_cnt) % mod)

Solution (Python):
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 ceil, floor, cos, radians, pi, sin
from operator import mul
from functools import reduce
from operator import mul
from functools import lru_cache

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)]


@lru_cache(maxsize=None)
def factorial(n):
    if n == 0:
        return 1
    else:
        return (n*factorial(n-1)) % mod


class UnionFind():
    def __init__(self, n):
        self.n = n
        self.parents = [-1] * n

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

    def union(self, x, y):
        x = self.find(x)
        y = self.find(y)

        if x == y:
            return

        if self.parents[x] > self.parents[y]:
            x, y = y, x

        self.parents[x] += self.parents[y]
        self.parents[y] = x

    def size(self, x):
        return -self.parents[self.find(x)]

    def same(self, x, y):
        return self.find(x) == self.find(y)

    def members(self, x):
        root = self.find(x)
        return [i for i in range(self.n) if self.find(i) == root]

    def roots(self):
        return [i for i, x in enumerate(self.parents) if x < 0]

    def group_count(self):
        return len(self.roots())

    def all_group_members(self):
        return {r: self.members(r) for r in self.roots()}


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):
    
    buf = dp[i] * factorial(i - 1)
    buf %= mod
    buf *= pow(n - 1, k - i, mod)
    buf %= mod
    cycle_count += buf



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

Solution (Python):
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

    
    def find(self,x):
        if self.par[x] == x:
            return x
        else:
            self.par[x] = y = self.find(self.par[x])
            return 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

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

    
    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:
    ans = 0
    for i in range(N):
        if i != uf.find(i):
            ans += 1
    print(ans)
    exit()

K = len(n)


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


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



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 = [[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 = [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)

Solution (Python):
import sys
INF = 1 << 60
MOD = 10**9 + 7 
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)

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

    
    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)

    
    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

    
    mf = modfact(5001)
    for i in range(1, l + 1):
        if i == 1:
            ans -= sum(s - 1 for s in S) % MOD * pow(n - 1, k - 1, MOD)
        else:
            ans -= dp[i] * mf.fact(i - 1) % MOD * pow(n - 1, k - i, MOD)
        ans %= MOD

    print(ans)
resolve()

Solution (Python):
import sys
import numpy as np

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

MOD = 10**9 + 7

N, *P = map(int, read().split())

class UnionFind:
    def __init__(self, N):
        self.root = list(range(N))
        self.size = [1] * (N)
        self.has_loop = [False] * N

    def __getitem__(self, x):
        root = self.root
        while root[x] != x:
            root[x] = root[root[x]]
            x = root[x]
        return x

    def merge(self, x, y):
        x = self[x]
        y = self[y]
        if x == y:
            self.has_loop[x] = True
            return False
        sx, sy = self.size[x], self.size[y]
        if sx < sy:
            x, y = y, x
            sx, sy = sy, sx
        self.root[y] = x
        self.size[x] += sy
        self.has_loop[x] |= self.has_loop[y]
        return True

uf = UnionFind(N + 1)
for i, x in enumerate(P, 1):
    if x == -1:
        continue
    uf.merge(i, x)

loop = 0
sizes = []
for v in range(1, N + 1):
    if v != uf[v]:
        continue
    if uf.has_loop[v]:
        loop += 1
    else:
        sizes.append(uf.size[v])

S = len(sizes)

f = np.zeros(S + 10, np.int64)
f[0] = 1
for x in sizes:
    f[1:] += x * f[:-1]
    f %= MOD
f[1] -= S

fact = np.empty(N + 1, np.int64)
fact[0] = 1
for n in range(1, N + 1):
    fact[n] = n * fact[n - 1] % MOD
power = np.empty(N + 1, np.int64)
power[0] = 1
for n in range(1, N + 1):
    power[n] = (N - 1) * power[n - 1] % MOD

coef = fact[:S] * power[:S][::-1] % MOD

x = (f[1:S+1] * coef % MOD).sum() % MOD
x += loop * pow(N-1, S, MOD) % MOD

x = N * pow(N-1, S, MOD) - x
x %= MOD
print(x)

Solution (Python):
import sys
input = lambda : sys.stdin.readline().rstrip()
sys.setrecursionlimit(max(1000, 10**9))
write = lambda x: sys.stdout.write(x+"\n")



class UnionFindTree:
    def __init__(self, n):
        self.n = n
        self.parent = list(range(n))
        self.size = [1] * n

    def root(self, i):
        inter = set()
        while self.parent[i]!=i:
            inter.add(i)
            i = self.parent[i]
        r = i
        for i in inter:
            self.parent[i] = r
        return r

    def connect(self, i, j):
        if i==j:
            return
        ri = self.root(i)
        rj = self.root(j)
        if ri==rj:
            return
        if self.size[ri]<self.size[rj]:
            self.parent[ri] = rj
            self.size[rj] += self.size[ri]
        else:
            self.parent[rj] = ri
            self.size[ri] += self.size[rj]


n = int(input())
uf = UnionFindTree(n)
ps = list(map(int, input().split()))
k = 0
for i,p in enumerate(ps):
    if p<0:
        k += 1
        continue
    uf.connect(i,p-1)
for i in range(n):
    uf.root(i)
from collections import Counter, defaultdict
M = 10**9+7
ncount = Counter(uf.parent)
mcount = defaultdict(int)
for i,g in enumerate(uf.parent):
    if ps[i]==-1:
        mcount[g] += 1

ans = sum(v-1 for v in ncount.values())
ans *= pow(n-1, k, M)

for kk in ncount.keys():
    ncount[kk] -= mcount[kk]

for kk,mm in mcount.items():
    if mm==0:
        continue
    nn = ncount[kk]
    ans += (n - mm - nn)*pow(n-1, k-1, M)*mm
    ans %= M

nums = [nn+mcount[k] for k,nn in ncount.items() if mcount[k]]
l = len(nums)
dp = [[0] * l for _ in range(l)]

if nums:
    dp[0][0] = nums[0]
    for i in range(1,l):
        dp[i][0] = dp[i-1][0] + nums[i]

    for i in range(1,l):
        for j in range(1,l):
            dp[i][j] += dp[i-1][j-1] * nums[i] * j + dp[i-1][j]
            dp[i][j] %= M

    for j in range(1, l):
        ans -= dp[-1][j]*pow(n-1, l-j-1, M)

print(ans%M)

Solution (Python):

import sys
sys.setrecursionlimit(10**8)
INF = float("inf")

MOD = 1000000007  


class UnionFind():
    def __init__(self, n):
        self.n = n
        self.parents = [-1] * n

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

    def union(self, x, y):
        x = self.find(x)
        y = self.find(y)

        if x == y:
            return

        if self.parents[x] > self.parents[y]:
            x, y = y, x

        self.parents[x] += self.parents[y]
        self.parents[y] = x

    def size(self, x):
        return -self.parents[self.find(x)]

    def same(self, x, y):
        return self.find(x) == self.find(y)

    def members(self, x):
        root = self.find(x)
        return [i for i in range(self.n) if self.find(i) == root]

    def roots(self):
        return [i for i, x in enumerate(self.parents) if x < 0]

    def group_count(self):
        return len(self.roots())

    def all_group_members(self):
        return {r: self.members(r) for r in self.roots()}

    def __str__(self):
        return '\n'.join('{}: {}'.format(r, self.members(r)) for r in self.roots())


def main():
    N = int(input())
    P = list(map(int, input().split()))

    
    uf = UnionFind(N)
    K = 0
    for i, p in enumerate(P):
        if p != -1:
            uf.union(i, p-1)
        else:
            K += 1

    
    factorial = [1]*K
    for i in range(K-1):
        factorial[i+1] = ((i+1)*factorial[i]) % MOD

    
    size_undef = [0]*K
    curr = 0
    for i, p in enumerate(P):
        if p == -1:
            size_undef[curr] = uf.size(i)
            curr += 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]*size_undef[i-1]
        dp = dp_new[:]

    
    if K >= 1:
        tot = 0
        for s in size_undef:
            tot += s-1
        dp[1] = tot

    
    cycle_count = 0
    for i in range(1, K+1):
        
        cycle_count += dp[i]*factorial[i-1]*pow(N-1, K-i, MOD)
        cycle_count %= MOD

    
    cycle_count += (uf.group_count()-K)*pow(N-1, K, MOD)

    
    ans = N*pow(N-1, K, MOD) - cycle_count
    ans %= MOD

    print(ans)

    return


if __name__ == '__main__':
    main()