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

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]
        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: Accepted

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
        # ランクの取得
        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)

Submission status: Accepted

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

Submission status: Accepted

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 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):
    # i個のグループからなるサイクルがひとつできる
    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)

Submission status: Accepted

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

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)

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

Submission status: Accepted

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

Submission status: Accepted

Solution code:
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)]
# dp[i][j] : i番目まで使うときの周期j+1の個数
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)

Submission status: Accepted

Solution code:
#!/usr/bin/env python3
import sys
sys.setrecursionlimit(10**8)
INF = float("inf")

MOD = 1000000007  # type: int


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[i, j]: i個の根付き木を見てj個を選んだ時の\sum{\prod_v S_v}
    # 「メモリ削減したdpとN^2確保するDP、どっちが早いんやろなあ」
    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[:]

    # ただしj=1個を選んだ時のものは特別にsize[i]-1を考える必要がある
    if K >= 1:
        tot = 0
        for s in size_undef:
            tot += s-1
        dp[1] = tot

    # cycleをカウントする
    cycle_count = 0
    for i in range(1, K+1):
        # i個のグループからなるサイクルがひとつできる
        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()

Submission status: Accepted

Solution code:
#!/usr/bin/env python3
import sys
sys.setrecursionlimit(10**8)
INF = float("inf")

MOD = 1000000007  # type: int


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[i, j]: i個の根付き木を見てj個を選んだ時の\sum{\prod_v S_v}
    # 「メモリ削減したdpとN^2確保するDP、どっちが早いんやろなあ」
    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[:]

    # ただしj=1個を選んだ時のものは特別にsize[i]-1を考える必要がある
    if K >= 1:
        tot = 0
        for s in size_undef:
            tot += s-1
        dp[1] = tot

    # cycleをカウントする
    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

    # すでに固定サイクルがあれば数える
    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()