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

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
    undef = []                  
    for i, p in enumerate(P):
        if p != -1:
            uf.union(i, p-1)
        else:
            K += 1
            undef.append(i)

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

    
    size_undef = []
    for u in undef:
        size_undef.append(uf.size(u))

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

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

Solution (Python):

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

MOD = 1000000007  


@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()}

    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
    undef = []
    for i, p in enumerate(P):
        if p != -1:
            uf.union(i, p-1)
        else:
            K += 1
            undef.append(i)

    groups = {}
    mapping = {}
    for i, r in enumerate(uf.roots()):
        groups[r] = [uf.size(r), True]
        mapping[i] = r
    for u in undef:
        groups[uf.find(u)][1] = False

    
    size_all = []
    looped = []
    size_undef = []
    for k, v in groups.items():
        size_all.append(v[0])
        looped.append(v[1])
        if v[1] == False:
            size_undef.append(v[0])
    
    

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

    
    cycle_count += (len(groups)-K)*pow(N-1, K, MOD)

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

    print(ans)

    return


if __name__ == '__main__':
    main()

Solution (Python):
def productall(a):
    N = len(a)
    if N == 0: return [1]
    A = [1] * N + a[:]
    P = 10 ** 9 + 7
    k = 12
    K = k * 8
    pa1 = (1 << k * 4 + 16) - ((1 << k * 4 + 16) % P)
    pa2 = (1 << k * 2 + 24) - ((1 << k * 2 + 24) % P)
    pa3 = (1 << k + 28) - ((1 << k + 28) % P)
    def modP(x):
        x -= ((x & m1) >> k * 4 + 16) * pa1
        x -= ((x & m2) >> k * 2 + 24) * pa2
        x -= ((x & m3) >> k + 28) * pa3
        return x
    
    c = 3
    m1 = int(("1" * (k * 4 - 16) + "0" * (k * 4 + 16)) * c, 2)
    m2 = int(("1" * (k * 6 - 24) + "0" * (k * 2 + 24)) * c, 2)
    m3 = int(("1" * (k * 7 - 28) + "0" * (k + 28)) * c, 2)
    for i in range(1, N)[::-1]:
        if i == 1 << i.bit_length() - 1:
            c *= 2
            m1 = int(("1" * (k * 4 - 16) + "0" * (k * 4 + 16)) * c, 2)
            m2 = int(("1" * (k * 6 - 24) + "0" * (k * 2 + 24)) * c, 2)
            m3 = int(("1" * (k * 7 - 28) + "0" * (k + 28)) * c, 2)
        A[i] = modP(A[2*i] * A[2*i+1])

    t = bin(A[1])[2:] + "_"
    return [int(t[-(i+1) * K - 1:-i * K - 1], 2) % P for i in range((len(t)+K-2) // K)]
def par(a):
    L = []
    while P[a] != a:
        L.append(a)
        a = P[a]
    for l in L:
        P[l] = a
    return a
def unite(a, b):
    pa = par(a)
    pb = par(b)
    if pa == pb: return
    if LEN[pa] < LEN[pb]:
        a, b, pa, pb = b, a, pb, pa
    P[pb] = pa
    if LEN[pa] == LEN[pb]: LEN[pa] += 1
    CNT[pa] += CNT[pb]
def cnt(a):
    return CNT[par(a)]

N = int(input())
P = [i for i in range(N)]
LEN = [1] * N
CNT = [1] * N
FULL = [0] * N
A = [int(a) - 1 for a in input().split()]
for i, a in enumerate(A):
    if a < 0: continue
    if par(i) != par(a):
        unite(i, a)
    else:
        FULL[i] = 1
for i in range(N):
    if FULL[i]:
        FULL[par(i)] = 1
X = []
Y = []
for i in range(N):
    if par(i) == i:
        if FULL[i] == 0:
            X.append(CNT[i])
        else:
            Y.append(CNT[i])
M = len(X)
mod = 10 ** 9 + 7
ans = (sum(X) + sum(Y) - len(Y)) * pow(N - 1, M, mod) % mod
L = productall([(a << 96) + 1 for a in X])
fa = 1
if M: ans = (ans + M * pow(N - 1, M - 1, mod)) % mod
for i, l in enumerate(L):
    if i == 0: continue
    ans = (ans - l * fa * pow(N - 1, M - i + mod - 1, mod)) % mod
    fa = fa * i % mod
print(ans)

Solution (Python):
def productall(a):
    N = len(a)
    if N == 0: return [1]
    A = [1] * N + a[:]
    P = 10 ** 9 + 7
    k = 12
    K = k * 8
    pa1 = (1 << k * 4 + 16) - ((1 << k * 4 + 16) % P)
    pa2 = (1 << k * 2 + 24) - ((1 << k * 2 + 24) % P)
    pa3 = (1 << k + 28) - ((1 << k + 28) % P)
    def modP(x):
        x -= ((x & m1) >> k * 4 + 16) * pa1
        x -= ((x & m2) >> k * 2 + 24) * pa2
        x -= ((x & m3) >> k + 28) * pa3
        return x
    
    c = 4
    m1 = int(("1" * (k * 4 - 16) + "0" * (k * 4 + 16)) * c, 2)
    m2 = int(("1" * (k * 6 - 24) + "0" * (k * 2 + 24)) * c, 2)
    m3 = int(("1" * (k * 7 - 28) + "0" * (k + 28)) * c, 2)
    for i in range(1, N)[::-1]:
        if i == 1 << i.bit_length() - 1:
            c *= 2
            m1 = int(("1" * (k * 4 - 16) + "0" * (k * 4 + 16)) * c, 2)
            m2 = int(("1" * (k * 6 - 24) + "0" * (k * 2 + 24)) * c, 2)
            m3 = int(("1" * (k * 7 - 28) + "0" * (k + 28)) * c, 2)
        A[i] = modP(A[2*i] * A[2*i+1])

    t = bin(A[1])[2:] + "_"
    return [int(t[-(i+1) * K - 1:-i * K - 1], 2) % P for i in range((len(t)+K-2) // K)]
def par(a):
    L = []
    while P[a] != a:
        L.append(a)
        a = P[a]
    for l in L:
        P[l] = a
    return a
def unite(a, b):
    pa = par(a)
    pb = par(b)
    if pa == pb: return
    if LEN[pa] < LEN[pb]:
        a, b, pa, pb = b, a, pb, pa
    P[pb] = pa
    if LEN[pa] == LEN[pb]: LEN[pa] += 1
    CNT[pa] += CNT[pb]
def cnt(a):
    return CNT[par(a)]

N = int(input())
P = [i for i in range(N)]
LEN = [1] * N
CNT = [1] * N
FULL = [0] * N
A = [int(a) - 1 for a in input().split()]
for i, a in enumerate(A):
    if a < 0: continue
    if par(i) != par(a):
        unite(i, a)
    else:
        FULL[i] = 1
for i in range(N):
    if FULL[i]:
        FULL[par(i)] = 1
X = []
Y = []
for i in range(N):
    if par(i) == i:
        if FULL[i] == 0:
            X.append(CNT[i])
        else:
            Y.append(CNT[i])
M = len(X)
mod = 10 ** 9 + 7
ans = (sum(X) + sum(Y) - len(Y)) * pow(N - 1, M, mod) % mod
L = productall([(a << 96) + 1 for a in X])
fa = 1
if M: ans = (ans + M * pow(N - 1, M - 1, mod)) % mod
for i, l in enumerate(L):
    if i == 0: continue
    ans = (ans - l * fa * pow(N - 1, M - i + mod - 1, mod)) % mod
    fa = fa * i % mod
print(ans)

Solution (Python):
def productall(a):
    N = len(a)
    if N == 0: return [1]
    A = [1] * N + a[:]
    P = 10 ** 9 + 7
    k = 12
    K = k * 8
    pa1 = (1 << k * 4 + 16) - ((1 << k * 4 + 16) % P)
    pa2 = (1 << k * 2 + 24) - ((1 << k * 2 + 24) % P)
    pa3 = (1 << k + 28) - ((1 << k + 28) % P)
    m1 = int(("1" * (k * 4 - 16) + "0" * (k * 4 + 16)) * 5050, 2)
    m2 = int(("1" * (k * 6 - 24) + "0" * (k * 2 + 24)) * 5050, 2)
    m3 = int(("1" * (k * 7 - 28) + "0" * (k + 28)) * 5050, 2)
    def modP(x):
        x -= ((x & m1) >> k * 4 + 16) * pa1
        x -= ((x & m2) >> k * 2 + 24) * pa2
        x -= ((x & m3) >> k + 28) * pa3
        return x

    for i in range(N)[::-1]:
        A[i] = modP(A[2*i] * A[2*i+1])

    t = bin(A[1])[2:] + "_"
    return [int(t[-(i+1) * K - 1:-i * K - 1], 2) % P for i in range((len(t)+K-2) // K)]
def par(a):
    L = []
    while P[a] != a:
        L.append(a)
        a = P[a]
    for l in L:
        P[l] = a
    return a
def unite(a, b):
    pa = par(a)
    pb = par(b)
    if pa == pb: return
    if LEN[pa] < LEN[pb]:
        a, b, pa, pb = b, a, pb, pa
    P[pb] = pa
    if LEN[pa] == LEN[pb]: LEN[pa] += 1
    CNT[pa] += CNT[pb]
def cnt(a):
    return CNT[par(a)]

N = int(input())
P = [i for i in range(N)]
LEN = [1] * N
CNT = [1] * N
FULL = [0] * N
A = [int(a) - 1 for a in input().split()]
for i, a in enumerate(A):
    if a < 0: continue
    if par(i) != par(a):
        unite(i, a)
    else:
        FULL[i] = 1
for i in range(N):
    if FULL[i]:
        FULL[par(i)] = 1
X = []
Y = []
for i in range(N):
    if par(i) == i:
        if FULL[i] == 0:
            X.append(CNT[i])
        else:
            Y.append(CNT[i])
M = len(X)
mod = 10 ** 9 + 7
ans = (sum(X) + sum(Y) - len(Y)) * pow(N - 1, M, mod) % mod
L = productall([(a << 96) + 1 for a in X])
fa = 1
if M: ans = (ans + M * pow(N - 1, M - 1, mod)) % mod
for i, l in enumerate(L):
    if i == 0: continue
    ans = (ans - l * fa * pow(N - 1, M - i + mod - 1, mod)) % mod
    fa = fa * i % mod
print(ans)

Solution (Python):
def productall(a):
    N = len(a)
    if N == 0: return [1]
    A = [1] * N + a[:]
    P = 10 ** 9 + 7
    k = 12
    K = k * 8
    pa1 = (1 << k * 4 + 16) - ((1 << k * 4 + 16) % P)
    pa2 = (1 << k * 2 + 24) - ((1 << k * 2 + 24) % P)
    pa3 = (1 << k + 28) - ((1 << k + 28) % P)
    m1 = int(("1" * (k * 4 - 16) + "0" * (k * 4 + 16)) * 5050, 2)
    m2 = int(("1" * (k * 6 - 24) + "0" * (k * 2 + 24)) * 5050, 2)
    m3 = int(("1" * (k * 7 - 28) + "0" * (k + 28)) * 5050, 2)
    def modP(x):
        x -= ((x & m1) >> k * 4 + 16) * pa1
        x -= ((x & m2) >> k * 2 + 24) * pa2
        x -= ((x & m3) >> k + 28) * pa3
        return x

    for i in range(N)[::-1]:
        A[i] = modP(A[2*i] * A[2*i+1])

    t = bin(A[1])[2:] + "_"
    return [int(t[-(i+1) * K - 1:-i * K - 1], 2) % P for i in range((len(t)+K-2) // K)]
def par(a):
    L = []
    while P[a] != a:
        L.append(a)
        a = P[a]
    for l in L:
        P[l] = a
    return a
def unite(a, b):
    pa = par(a)
    pb = par(b)
    if pa == pb: return
    if LEN[pa] < LEN[pb]:
        a, b, pa, pb = b, a, pb, pa
    P[pb] = pa
    if LEN[pa] == LEN[pb]: LEN[pa] += 1
    CNT[pa] += CNT[pb]
def cnt(a):
    return CNT[par(a)]

N = int(input())
P = [i for i in range(N)]
LEN = [1] * N
CNT = [1] * N
FULL = [0] * N
A = [int(a) - 1 for a in input().split()]
for i, a in enumerate(A):
    if a < 0: continue
    if par(i) != par(a):
        unite(i, a)
    else:
        FULL[i] = 1
for i in range(N):
    if FULL[i]:
        FULL[par(i)] = 1
X = []
Y = []
for i in range(N):
    if par(i) == i:
        if FULL[i] == 0:
            X.append(CNT[i])
        else:
            Y.append(CNT[i])
M = len(X)
mod = 10 ** 9 + 7
ans = (sum(X) + sum(Y) - len(Y)) * pow(N - 1, M, mod) % mod
L = productall([(a << 96) + 1 for a in X])
if M:
    fa = 1
    ans = (ans + M * pow(N - 1, M - 1, mod)) % mod
    for i, l in enumerate(L):
        if i == 0: continue
        ans = (ans - l * fa * pow(N - 1, M - i, mod)) % mod
        fa = fa * i % mod
print(ans)

Solution (Python):
import sys;input=sys.stdin.readline
class UnionFindNode:
    def __init__(self, group_id, parent=None, value=None):
        self.group_id_ = group_id
        self.parent_ = parent
        self.value = value
        self.rank_ = 1
        self.member_num_ = 1

    def is_root(self):
        return not self.parent_

    def root(self):
        parent = self
        while not parent.is_root():
            parent = parent.parent_
            self.parent_ = parent
        return parent

    def find(self):
        root = self.root()
        return root.group_id_

    def rank(self):
        root = self.root()
        return root.rank_

    def unite(self, unite_node):
        root = self.root()
        unite_root = unite_node.root()

        if root.group_id_ != unite_root.group_id_:
            if root.rank() > unite_root.rank():
                n_root, child = root, unite_root
            else:
                n_root, child = unite_root, root
            child.parent_ = n_root
            n_root.rank_ = max(n_root.rank_, child.rank_ + 1)
            n_root.member_num_ = n_root.member_num_ + child.member_num_

mod = 10**9+7
NNN = (10**6)
g1 = [1, 1]
for i in range( 2, NNN + 1 ):
    g1.append( ( g1[-1] * i ) % mod )
N, = map(int, input().split())
P = list(map(int, input().split()))
K = P.count(-1)

nodes = [UnionFindNode(i) for i in range(N)]
for i, p in enumerate(P):
    p-=1
    if p >= 0:
        nodes[p].unite(nodes[i])

rs = []
X = [0] * N
for i in range(N):
    if nodes[i].is_root():
        rs.append(i)
    if P[i] == -1:
        X[nodes[i].find()] = 1

M = len(rs)
R = N*pow(N-1, K, mod)%mod

Y = []
for i in range(M):
    if X[rs[i]]:
        Y.append(nodes[rs[i]].member_num_)
        R = (R-(nodes[rs[i]].member_num_-1)*pow(N-1, K-1, mod))%mod
    else:
        R = (R-pow(N-1, K, mod))%mod
nY = len(Y)
dp = [[0]*(nY+1) for _ in range(nY+1)]
dp[0][0] = 1
for i in range(nY):
    for j in range(nY):
        dp[i+1][j] = (dp[i][j] * Y[i] + (dp[i][j-1] if j > 0 else 0)) % mod
for i in range(nY):
    if nY-i == 1:
        continue

    R -= pow(N-1, K-(nY-i), mod) * g1[nY-i-1] * dp[-1][i]
    R %= mod
print(R)

Solution (Python):
import sys
sys.setrecursionlimit(10**9)
mod = 10**9+7
from collections import defaultdict
class UnionFind:
    def __init__(self, num):
        self.table = [-1 for _ in range(num)]

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

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

        if s1 != s2:
            if self.table[s1] <= self.table[s2]:
                self.table[s1] += self.table[s2]
                self.table[s2] = s1
            else:
                self.table[s2] += self.table[s1]
                self.table[s1] = s2
            return True
        return False

memo = [-1]*5001
memo[0] = 1
def fact(n):
    if memo[n] != -1:
        return memo[n]
    else:
        res = n*fact(n-1)%mod
        memo[n] = res
        return res

N = int(input())
P = list(map(int,input().split()))
K = 0
u = UnionFind(N)
loop = 0
for i in range(N):
    if P[i]==-1:
        K += 1
        continue
    if u.find(i)==u.find(P[i]-1):
        loop += 1
    u.union(i,P[i]-1)

d = defaultdict(lambda:[0,0])
for i in range(N):
    if P[i]==-1:
        d[u.find(i)][0] += 1
    d[u.find(i)][1] += 1

n = []
for k in d:
    n.append(d[k])

dp = [[0]*(len(n)+1) for i in range(len(n))]
for i in range(len(n)):
    dp[i][0] = 1
for i in range(len(n)):
    for j in range(1,i+2):
        if j != i+1:
            dp[i][j] += dp[i-1][j]
        dp[i][j] += dp[i-1][j-1]*n[i][0]*n[i][1]
        dp[i][j] %= mod

ans = 0
for i in range(2,len(n)+1):
    if K-i<0:
        break
    ans += dp[len(n)-1][i]*fact(i-1)*pow(N-1,K-i,mod)
    ans %= mod
for i in range(len(n)):
    if K==0:
        break
    ans += n[i][0]*(n[i][1]-1)*pow(N-1,K-1,mod)
    ans %= mod

ans = pow(N-1,K,mod)*(N-loop)-ans
ans %= mod
print(ans)

Solution (Python):
def productall(a):
    N = len(a)
    if N == 0: return [1]
    A = [1] * N + a[:]
    P = 10 ** 9 + 7
    k = 100
    K = k * 8
    pa1 = (1 << k * 4 + 16) - ((1 << k * 4 + 16) % P)
    pa2 = (1 << k * 2 + 24) - ((1 << k * 2 + 24) % P)
    pa3 = (1 << k + 28) - ((1 << k + 28) % P)
    def modP(x):
        x -= ((x & m1) >> k * 4 + 16) * pa1
        x -= ((x & m2) >> k * 2 + 24) * pa2
        x -= ((x & m3) >> k + 28) * pa3
        return x

    for i in range(N)[::-1]:
        if i == N - 1 or i and 1 << i.bit_length() - 1 == i:
            if i == N - 1:
                c = 2
            else:
                c *= 2
            m1 = int(("1" * (k * 4 - 16) + "0" * (k * 4 + 16)) * c, 2)
            m2 = int(("1" * (k * 6 - 24) + "0" * (k * 2 + 24)) * c, 2)
            m3 = int(("1" * (k * 7 - 28) + "0" * (k + 28)) * c, 2)
        A[i] = modP(A[2*i] * A[2*i+1])

    t = bin(A[1])[2:] + "_"
    return [int(t[-(i+1) * K - 1:-i * K - 1], 2) % P for i in range((len(t)+K-2) // K)]
def par(a):
    L = []
    while P[a] != a:
        L.append(a)
        a = P[a]
    for l in L:
        P[l] = a
    return a
def unite(a, b):
    pa = par(a)
    pb = par(b)
    if pa == pb: return
    if LEN[pa] < LEN[pb]:
        a, b, pa, pb = b, a, pb, pa
    P[pb] = pa
    if LEN[pa] == LEN[pb]: LEN[pa] += 1
    CNT[pa] += CNT[pb]
def cnt(a):
    return CNT[par(a)]

N = int(input())
P = [i for i in range(N)]
LEN = [1] * N
CNT = [1] * N
FULL = [0] * N
A = [int(a) - 1 for a in input().split()]
for i, a in enumerate(A):
    if a < 0: continue
    if par(i) != par(a):
        unite(i, a)
    else:
        FULL[i] = 1
for i in range(N):
    if FULL[i]:
        FULL[par(i)] = 1
X = []
Y = []
for i in range(N):
    if par(i) == i:
        if FULL[i] == 0:
            X.append(CNT[i])
        else:
            Y.append(CNT[i])
M = len(X)
mod = 10 ** 9 + 7
ans = (sum(X) + sum(Y) - len(Y)) * pow(N - 1, M, mod) % mod
L = productall([(a << 800) + 1 for a in X])
if M:
    fa = 1
    ans = (ans + M * pow(N - 1, M - 1, mod)) % mod
    for i, l in enumerate(L):
        if i == 0: continue
        ans = (ans - l * fa * pow(N - 1, M - i, mod)) % mod
        fa = fa * i % mod
print(ans)

Solution (Python):
def par(a):
    L = []
    while P[a] != a:
        L.append(a)
        a = P[a]
    for l in L:
        P[l] = a
    return a
def unite(a, b):
    pa = par(a)
    pb = par(b)
    if pa == pb: return
    if LEN[pa] < LEN[pb]:
        a, b, pa, pb = b, a, pb, pa
    P[pb] = pa
    if LEN[pa] == LEN[pb]: LEN[pa] += 1
    CNT[pa] += CNT[pb]
def cnt(a):
    return CNT[par(a)]

N = int(input())
P = [i for i in range(N)]
LEN = [1] * N
CNT = [1] * N
FULL = [0] * N
A = [int(a) - 1 for a in input().split()]
for i, a in enumerate(A):
    if a < 0: continue
    if par(i) != par(a):
        unite(i, a)
    else:
        FULL[i] = 1
for i in range(N):
    if FULL[i]:
        FULL[par(i)] = 1
X = []
Y = []
for i in range(N):
    if par(i) == i:
        if FULL[i] == 0:
            X.append(CNT[i])
        else:
            Y.append(CNT[i])
M = len(X)
mod = 10 ** 9 + 7
K = 96
m = int(("1" * 32 + "0" * 64) * 5050, 2)
pa = (1 << 64) - ((1 << 64) % mod)
modP = lambda x: x - ((x & m) >> 64) * pa
ans = (sum(X) + sum(Y) - len(Y)) * pow(N - 1, M, mod) % mod
x = 1
for i, a in enumerate(X):
    x *= (a << K) + 1
    x = modP(x)

sx = bin(x)[2:] + "_"
L = [int(sx[-(i+1) * K - 1:-i * K - 1], 2) % mod for i in range((len(sx)+K-2) // K)]
if M:
    fa = 1
    ans = (ans + M * pow(N - 1, M - 1, mod)) % mod
    for i, l in enumerate(L):
        if i == 0: continue
        ans = (ans - l * fa * pow(N - 1, M - i, mod)) % mod
        fa = fa * i % mod
print(ans)

Solution (Python):
def par(a):
    L = []
    while P[a] != a:
        L.append(a)
        a = P[a]
    for l in L:
        P[l] = a
    return a
def unite(a, b):
    pa = par(a)
    pb = par(b)
    if pa == pb: return
    if LEN[pa] < LEN[pb]:
        a, b, pa, pb = b, a, pb, pa
    P[pb] = pa
    if LEN[pa] == LEN[pb]: LEN[pa] += 1
    CNT[pa] += CNT[pb]
def cnt(a):
    return CNT[par(a)]

N = int(input())
P = [i for i in range(N)]
LEN = [1] * N
CNT = [1] * N
FULL = [0] * N
A = [int(a) - 1 for a in input().split()]
for i, a in enumerate(A):
    if a < 0: continue
    if par(i) != par(a):
        unite(i, a)
    else:
        FULL[i] = 1
for i in range(N):
    if FULL[i]:
        FULL[par(i)] = 1
X = []
Y = []
for i in range(N):
    if par(i) == i:
        if FULL[i] == 0:
            X.append(CNT[i])
        else:
            Y.append(CNT[i])
M = len(X)
mod = 10 ** 9 + 7
K = 96
m = int(("1" * 32 + "0" * 64) * 5050, 2)
pa = (1 << 64) - ((1 << 64) % mod)
modP = lambda x: x - ((x & m) >> 64) * pa
ans = (sum(X) + sum(Y) - len(Y)) * pow(N - 1, M, mod) % mod
x = 1
for i, a in enumerate(X):
    x *= (a << K) + 1
    x = modP(x)

sx = bin(x)[2:] + "_"
L = [int(sx[-(i+1) * K - 1:-i * K - 1], 2) % mod for i in range((len(sx)+K-2) // K)]
fa = 1
ans = (ans + M * pow(N - 1, M - 1, mod)) % mod
for i, l in enumerate(L):
    if i == 0: continue
    ans = (ans - l * fa * pow(N - 1, M - i, mod)) % mod
    fa = fa * i % mod
print(ans)

Solution (Python):
import numpy as np


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

    def _root(self, x):
        stack = []
        tbl = self.table
        while tbl[x] >= 0:
            stack.append(x)
            x = tbl[x]
        for y in stack:
            tbl[y] = x
        return x

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

    def unite(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.table[r1] += d2
        else:
            self.table[r1] = r2
            self.table[r2] += d1

    def get_size(self, x):
        return -self.table[self._root(x)]


n = int(input())
ppp = list(map(int, input().split()))
MOD = 10 ** 9 + 7
uft = UnionFind(n)
base = 0
undefined = []
for i, p in enumerate(ppp):
    if p == -1:
        undefined.append(i)
        continue
    p -= 1
    if not uft.find(i, p):
        base += 1
        uft.unite(i, p)

if len(undefined) == 0:
    print(base % MOD)
    exit()
elif len(undefined) == 1:
    c = uft.get_size(undefined[0])
    others = n - c
    print((base * (n - 1) + others) % MOD)
    exit()

m = len(undefined)
dp = np.zeros(m + 1, dtype=np.int64)  
dp[0] = 1
additional = 0  

for i in undefined:
    c = uft.get_size(i)
    dp[1:] += dp[:-1] * c
    dp %= MOD
    additional += n - c
dp = dp.tolist()

duplicated = 0
pat = pow(n - 1, m, MOD)
inv = pow(n - 1, MOD - 2, MOD)
loop_permutation = 1
loop_other_pattern = pat * inv * inv % MOD
for loop_size in range(2, m + 1):
    duplicated = (duplicated + dp[loop_size] * loop_permutation * loop_other_pattern) % MOD
    loop_other_pattern = loop_other_pattern * inv % MOD
    loop_permutation = loop_permutation * loop_size % MOD

ans = (base * pat + additional * pat * inv - duplicated) % MOD

print(ans)

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

    def _root(self, x):
        stack = []
        tbl = self.table
        while tbl[x] >= 0:
            stack.append(x)
            x = tbl[x]
        for y in stack:
            tbl[y] = x
        return x

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

    def unite(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.table[r1] += d2
        else:
            self.table[r1] = r2
            self.table[r2] += d1

    def get_size(self, x):
        return -self.table[self._root(x)]


n = int(input())
ppp = list(map(int, input().split()))
MOD = 10 ** 9 + 7
uft = UnionFind(n)
base = 0
undefined = []
for i, p in enumerate(ppp):
    if p == -1:
        undefined.append(i)
        continue
    p -= 1
    if not uft.find(i, p):
        base += 1
        uft.unite(i, p)

if len(undefined) == 0:
    print(base % MOD)
    exit()
elif len(undefined) == 1:
    c = uft.get_size(undefined[0])
    others = n - c
    print((base * (n - 1) + others) % MOD)
    exit()

m = len(undefined)
dp = [0] * (m + 1)  
dp[0] = 1
additional = 0  

for i in undefined:
    c = uft.get_size(i)
    for j in range(m - 1, -1, -1):
        dp[j + 1] += dp[j] * c
    additional += n - c

duplicated = 0
pat = pow(n - 1, m, MOD)
inv = pow(n - 1, MOD - 2, MOD)
loop_permutation = 1
loop_other_pattern = pat * inv * inv % MOD
for loop_size in range(2, m + 1):
    duplicated = (duplicated + dp[loop_size] * loop_permutation * loop_other_pattern) % MOD
    loop_other_pattern = loop_other_pattern * inv % MOD
    loop_permutation = loop_permutation * loop_size % MOD

ans = (base * pat + additional * pat * inv - duplicated) % MOD

print(ans)

Solution (Python):
def modfac(n, MOD):
 
    f = 1
    factorials = [1]
    for m in range(1, n + 1):
        f *= m
        f %= MOD
        factorials.append(f)
    inv = pow(f, MOD - 2, MOD)
    invs = [1] * (n + 1)
    invs[n] = inv
    for m in range(n, 1, -1):
        inv *= m
        inv %= MOD
        invs[m - 1] = inv
    return factorials, invs


def modnCr(n,r,mod,fac,inv): 

    return fac[n] * inv[n-r] * inv[r] % mod

def uf_find(n,p):

    ufl = []

    while p[n] != n:
        ufl.append(n)
        n = p[n]

    for i in ufl:
        p[i] = n

    return n


def uf_union(a,b,p,rank,tsize):

    ap = uf_find(a,p)
    bp = uf_find(b,p)

    if ap == bp:
        return True
    else:

        if rank[ap] > rank[bp]:
            p[bp] = ap
            tsize[ap] += tsize[bp]
        elif rank[ap] < rank[bp]:
            p[ap] = bp
            tsize[bp] += tsize[ap]
        else:
            p[bp] = ap
            rank[ap] += 1
            tsize[ap] += tsize[bp]

        return False

mod = 10**9+7
N = int(input())
P = list(map(int,input().split()))
fac,inv = modfac(N+10,mod)

ans = 0

p = [i for i in range(N)]
rank = [1] * N
tsize = [1] * N

for i in range(N):

    if P[i] == -1:
        continue
    nex = P[i] - 1
    uf_union(i,nex,p,rank,tsize)

able = [True] * N
c = []
d = []

for i in range(N):

    if P[i] == -1:
        nowp = uf_find(i,p)
        c.append(tsize[nowp])
        able[nowp] = False

for i in range(N):

    if uf_find(i,p) == i and able[i]:
        d.append(tsize[i])


ans = 0
for i in c:
    ans += i
for i in d:
    ans += i-1
ans *= pow(N-1,len(c),mod)
ans %= mod




for i in c:
    ans -= (i-1) * pow(N-1,len(c)-1,mod)
    ans %= mod



dp = [0] * (len(c)+1)
dp[0] = 1

for i in c:
    for j in range(len(c)-1,-1,-1):
        dp[j+1] += dp[j] * i
        


for i in range(2,len(c)+1):

    ans -= dp[i] * fac[i-1] * pow(N-1,len(c)-i,mod)
    ans %= mod

print (ans % mod)


"""
連結成分で考える
全ての街が -1 以外を選んだ時、
最小数 + -1の数
-1同士が選ぶことを考える
-1同士が選んだ時のみ答えが1減る
そうでなければ変化なし？

連結成分に-1は1つ以下しかない？
→正しい

-1を含まない連結成分は分けておく
-1の連結成分同士が求め合えば、答えは-1される

a個含む & bこ含む
場合

a*b * (N-1)^(K-2)
引けばいい

自分の連結成分を求めるときも1減る
ループを作ると1減るんだ！



"""