Alex has a board game consisting of:

A chip for marking his current location on the board.
$n$ fields numbered from $1$ to $n$. Each position $\boldsymbol{i}$ has a value, $f_i$, denoting the next position for the chip to jump to from that field.
A die with $m$ faces numbered from $\mbox{o}$ to $m-1$. Each face $j$ has a probability, $p_j$, of being rolled.

Alex then performs the following actions:

Begins the game by placing the chip at a position in a field randomly and with equiprobability. 
Takes $\boldsymbol{\mbox{k}}$ turns; during each turn he:
Rolls the die. We'll denote the number rolled during a turn as $\boldsymbol{d}$.
Jumps the chip $\boldsymbol{d}$ times. Recall that each field contains a value denoting the next field number to jump to.
After completing $\boldsymbol{\mbox{k}}$ turns, the game ends and he must calculate the respective probabilities for each field as to whether the game ended with the chip in that field.

Given $n$, $m$, $\boldsymbol{\mbox{k}}$, the game board, and the probabilities for each die face, print $n$ lines where each line $\boldsymbol{i}$ contains the probability that the chip is on field $\boldsymbol{i}$ at the end of the game.

Note: All the probabilities in this task are rational numbers modulo $M=998244353$. That is, if the probability can be expressed as the irreducible fraction $\frac{p}{q}$ where $q\ \text{mod}\ M\neq0$, then it corresponds to the number $(p\times q^{-1})\:\text{mod}\:M$ (or, alternatively, $p\times q^{-1}\equiv x(\:\text{mod}\:M)$). Click here to learn about Modular Multiplicative Inverse.

Input Format

The first line contains three space-separated integers describing the respective values of $n$ (the number of positions), $m$ (the number of die faces), and $\boldsymbol{\mbox{k}}$ (the number of turns). 

The second line contains $n$ space-separated integers describing the respective values of each $f_i$ (i.e., the index of the field that field $\boldsymbol{i}$ can transition to). 

The third line contains $m$ space-separated integers describing the respective values of each $p_j$ (where $0\leq p_j<M$) describing the probabilities of the faces of the $m$-sided die.  

Constraints

$1\leq n\leq6\times10^4$  
$4\leq m\leq10^5$  
$1\leq k\leq1000$  
$1\leq i,f_i\leq n$
$0\leq p_j<M$
The sum of $p_j\:\text{mod}\:M$ is $1$  

Note: The time limit for this challenge is doubled for all languages. Read more about standard time limits at our environment page.

Output Format

Print $n$ lines of output in which each line $\boldsymbol{i}$ contains a single integer, $x_i$ (where $0\leq x_i<M$), denoting the probability that the chip will be on field $\boldsymbol{i}$ after $\boldsymbol{\mbox{k}}$ turns.

Sample Input 0
4 5 1
2 3 2 4
332748118 332748118 332748118 0 0

Sample Output 0
582309206
332748118
332748118
748683265

Explanation 0

The diagram below depicts the respective probabilities of each die face being rolled:

The diagram below depicts each field with an arrow pointing to the next field:

There are four equiprobable initial fields, so each field has a $\frac{1}{4}$ probability of being the chip's initial location. Next, we calculate the probability that the chip will end up in each field after $k=1$ turn:

The only way the chip ends up in this field is if it never jumps from the field, which only happens if Alex rolls a $0$. So, this field's probability is $\frac{1}{4}\cdot\frac{1}{3}=\frac{1}{12}$. We then calculate and print the result of $\frac{1}{12}\ \text{mod}\ 998244353=582309206$ on a new line.

The chip can end up in field $2$ after one turn in the following scenarios:

Start in field $1$ and roll a $1$, the probability for which is $\frac{1}{4}\cdot\frac{1}{3}=\frac{1}{12}$.
Start in field $2$ and roll a $0$ or a $2$, the probability for which is $\frac{1}{4}\cdot\frac{2}{3}=\frac{2}{12}$.
Start in field $3$ and roll a $1$, the probability for which is $\frac{1}{4}\cdot\frac{1}{3}=\frac{1}{12}$.

After summing these probabilities, we get a total probability of $\frac1{12}+\frac2{12}+\frac1{12}=\frac13$ for the field. We then calculate and print the result of $\frac{1}{3}\ \text{mod}\ 998244353=332748118$ on a new line.

The chip can end up in field $3$ after one turn in the following scenarios:

Start in field $1$ and roll a $2$, the probability for which is $\frac{1}{4}\cdot\frac{1}{3}=\frac{1}{12}$.
Start in field $2$ and roll a $1$, the probability for which is $\frac{1}{4}\cdot\frac{1}{3}=\frac{1}{12}$.
Start in field $3$ and roll a $0$ or a $2$, the probability for which is $\frac{1}{4}\cdot\frac{2}{3}=\frac{2}{12}$.

After summing these probabilities, we get a total probability of $\frac1{12}+\frac1{12}+\frac2{12}=\frac13$ for the field. We then calculate and print the result of $\frac{1}{3}\ \text{mod}\ 998244353=332748118$ on a new line.

If the chip is initially placed in field $4$, it will always end up in field $4$ regardless of how many turns are taken (because this field loops back onto itself). Thus, this field's probability is $\frac{1}{4}$. We then calculate and print the result of $\frac{1}{4}\ \text{mod}\ 998244353=748683265$ on a new line.

Difficulty: HARD
Tags: Graph algorithms

Starter code:

Solutions:
MOD = 998244353

def inv(v):
	return pow(v % MOD, MOD - 2, MOD) % MOD

class Node(object):

	def __init__(self, i, p0, rolls):
		self.child = None
		self.parents = []
		self.lastProb = p0
		self.currentProb = 0
		self.index = i

	def clear(self):
		(self.lastProb, self.currentProb) = (self.currentProb, 0)

	def update(self, diceProbs, depth=0, baseProb=None):
		if baseProb == None:
			baseProb = self.lastProb
		if baseProb == 0 or depth >= len(diceProbs):
			return
		self.currentProb += diceProbs[depth] * baseProb
		self.currentProb %= MOD
		self.child.update(diceProbs, depth + 1, baseProb)

def DRW(nodeChildren, diceProbs, rolls):
	p0 = inv(len(nodeChildren))
	nodes = [Node(i, p0, rolls) for i in range(len(nodeChildren))]
	for (i, c) in enumerate(nodeChildren):
		nodes[i].child = nodes[c]
		nodes[c].parents.append(nodes[i])
	for r in range(1, rolls + 1):
		for n in nodes:
			n.update(diceProbs)
		for n in nodes:
			n.clear()
	return [n.lastProb for n in nodes]
import sys
data = sys.stdin.read().split('\n')
(n, m, k) = map(int, data[0].strip().split(' '))
children = list(map(lambda x: int(x) - 1, data[1].strip().split(' ')))
diceProbs = list(map(int, data[2].strip().split(' ')))
for p in DRW(children, diceProbs, k):
	print(p)


import sys

def egcd(a, b):
	if a == 0:
		return (b, 0, 1)
	(g, x, y) = egcd(b % a, a)
	return (g, y - b // a * x, x)

def mulinv(b, m):
	(g, x, _) = egcd(b, m)
	if g == 1:
		return x % m
(n, m, k) = (int(x) for x in next(sys.stdin).split())
fs = [int(x) - 1 for x in next(sys.stdin).split()]
ps = [int(x) for x in next(sys.stdin).split()]
M = 998244353

def jump(cells, fs):
	ncells = [0] * n
	for (f, c) in zip(fs, cells):
		ncells[f] += c
	return ncells
cells = [mulinv(n, M)] * n
for turn in range(k):
	jumps = [cells[:]]
	for _ in range(m - 1):
		jumps.append(jump(jumps[-1], fs))
		new_cells = [0] * n
	for (p, jcells) in zip(ps, jumps):
		for (i, c) in enumerate(jcells):
			new_cells[i] += p * c
	cells = [c % M for c in new_cells]
for c in cells:
	print(c)


import sys

def probs(next_pos, die_probs, num_positions, die_size, num_turns, mod):
	next_pos = [x - 1 for x in next_pos]
	step = [pow(num_positions, mod - 2, mod)] * num_positions
	for _ in range(num_turns):
		totals = [0] * num_positions
		for p in die_probs:
			for i in range(num_positions):
				totals[i] = (totals[i] + p * step[i]) % mod
			step2 = [0] * num_positions
			for i in range(num_positions):
				i2 = next_pos[i]
				step2[i2] = (step2[i2] + step[i]) % mod
			step = step2
		step = totals
	return totals
(n, m, k) = input().strip().split(' ')
(n, m, k) = [int(n), int(m), int(k)]
f = list(map(int, input().strip().split(' ')))
p = list(map(int, input().strip().split(' ')))
mod = 998244353
for x in probs(f, p, n, m, k, mod):
	print(x)


import sys
mod_base = 998244353

def gcd(x, y):
	while y != 0:
		(x, y) = (y, x % y)
	return x

def xgcd(b, n):
	(x0, x1, y0, y1) = (1, 0, 0, 1)
	while n != 0:
		(q, b, n) = (b // n, n, b % n)
		(x0, x1) = (x1, x0 - q * x1)
		(y0, y1) = (y1, y0 - q * y1)
	return (b, x0, y0)

def mulinv(b, n):
	(g, x, _) = xgcd(b, n)
	if g == 1:
		return x % n

def reduce(num, denom):
	d = gcd(num, denom)
	return (num // d, denom // d)

def prob_to_int(num, denom):
	(num, denom) = reduce(num, denom)
	return num * mulinv(denom, mod_base)
(n, m, k) = [int(num) for num in input().split()]
trans_to = [None] + [int(num) for num in input().split()]
die_prob = [int(num) for num in input().split()]

def find_paths(trans_to, n):
	path = {i: [i] for i in range(1, n + 1)}
	path_complete = [False] * (n + 1)
	origins = set([i for i in range(1, n + 1)])
	for i in range(1, n + 1):
		if trans_to[i] in origins:
			origins.remove(trans_to[i])
	for node in origins:
		(path, path_complete) = find_path_from_node(node, trans_to, n, path, path_complete)
	for node in range(1, n + 1):
		if not path_complete[node]:
			(path, path_complete) = find_path_from_node(node, trans_to, n, path, path_complete)
	return path

def find_path_from_node(node, trans_to, n, path, path_complete):
	cur_node = node
	cur_path = [cur_node]
	nodes_in_path = set([cur_node])
	while True:
		cur_node = trans_to[cur_node]
		if path_complete[cur_node]:
			for p in cur_path:
				path[p] = [path[p] + path[cur_node][0], path[cur_node][1]]
				path_complete[p] = True
			break
		elif cur_node in nodes_in_path:
			loop_idx = cur_path.index(cur_node)
			for i in range(loop_idx + 1, len(cur_path)):
				path[cur_path[i]] += cur_path[loop_idx:i]
			for (i, p) in enumerate(cur_path):
				if i < loop_idx:
					path[p] = [path[p][:loop_idx - i], path[p][loop_idx - i:]]
				else:
					path[p] = [[], path[p]]
				path_complete[p] = True
			break
		else:
			for n in cur_path:
				path[n].append(cur_node)
			cur_path.append(cur_node)
			nodes_in_path.add(cur_node)
	return (path, path_complete)
paths = find_paths(trans_to, n)

def destination(start, num_moves):
	(transient, cycle) = paths[start]
	if num_moves < len(transient):
		return transient[num_moves]
	else:
		return cycle[(num_moves - len(transient)) % len(cycle)]
trans_probs = {i: {} for i in range(1, n + 1)}
for node in range(1, n + 1):
	for num_moves in range(m):
		dest = destination(node, num_moves)
		trans_probs[node][dest] = (trans_probs[node].get(dest, 0) + die_prob[num_moves]) % mod_base
cell_probs = [None] + [prob_to_int(1, n) for _ in range(n)]

def update_probs(cell_probs, trans_probs):
	new_probs = [None] + [0 for _ in range(n)]
	for start_cell in range(1, n + 1):
		if cell_probs[start_cell] > 0:
			cell_prob = cell_probs[start_cell]
			for end_cell in trans_probs[start_cell]:
				trans_prob = trans_probs[start_cell][end_cell]
				new_prob = cell_prob * trans_prob % mod_base
				new_probs[end_cell] = (new_probs[end_cell] + new_prob) % mod_base
	return new_probs
for _ in range(k):
	cell_probs = update_probs(cell_probs, trans_probs)
for i in range(1, n + 1):
	print(cell_probs[i])


import sys
(n, m, k) = input().strip().split(' ')
(n, m, k) = [int(n), int(m), int(k)]
f = list(map(lambda x: int(x) - 1, input().strip().split(' ')))
p = list(map(int, input().strip().split(' ')))

def ExtendedEuclid(a, b):
	s = 0
	s_old = 1
	t = 1
	t_old = 0
	r = b
	r_old = a
	while r != 0:
		quot = r_old // r
		(r_old, r) = (r, r_old - quot * r)
		(s_old, s) = (s, s_old - quot * s)
		(t_old, t) = (t, t_old - quot * t)
	return t_old

def BruteForce(pField, f, p, n, M, TargetList):
	pFNew = [0 for i in range(n)]
	for (I_F, i_F) in enumerate(pField):
		for (I_D, i_D) in enumerate(p):
			pFNew[TargetList[I_F][I_D]] += i_F * i_D % M
	Out = [i % M for i in pFNew]
	return Out
M = 998244353
pStart = ExtendedEuclid(M, n)
pField = [pStart for m in range(n)]
TargetList = {i: [] for i in range(n)}
for StartPos in range(n):
	Pos = StartPos
	TargetList[StartPos].append(Pos)
	for MoveNum in range(1, m):
		Pos = f[Pos]
		TargetList[StartPos].append(Pos)
for kLoop in range(k):
	pField = BruteForce(pField, f, p, n, M, TargetList)
for i in pField:
	print(i)


import sys
(n, m, k) = input().strip().split(' ')
(n, m, k) = [int(n), int(m), int(k)]
f = list(map(int, input().strip().split(' ')))
p = list(map(int, input().strip().split(' ')))
M = 998244353
prob_init = [pow(n, M - 2, M)] * n
prob = [0] * n
for turn in range(k):
	for field in range(n):
		place = field + 1
		for face in range(m):
			if p[face] == 0:
				place = f[place - 1]
				continue
			prob[place - 1] = (prob[place - 1] + prob_init[field] * p[face] % M) % M
			place = f[place - 1]
	prob_init = prob
	prob = [0] * n
for pr in prob_init:
	print(pr)


MOD = 998244353

def inv(v):
	return pow(v % MOD, MOD - 2, MOD) % MOD

class Node(object):

	def __init__(self, i, p0, rolls):
		self.child = None
		self.parents = []
		self.lastProb = p0
		self.currentProb = 0
		self.index = i

	def clear(self):
		(self.lastProb, self.currentProb) = (self.currentProb, 0)

	def update(self, diceProbs, depth=0, baseProb=None):
		if baseProb == None:
			baseProb = self.lastProb
		if baseProb == 0 or depth >= len(diceProbs):
			return
		self.currentProb += diceProbs[depth] * baseProb
		self.currentProb %= MOD
		self.child.update(diceProbs, depth + 1, baseProb)

def DRW(nodeChildren, diceProbs, rolls):
	p0 = inv(len(nodeChildren))
	nodes = [Node(i, p0, rolls) for i in range(len(nodeChildren))]
	for (i, c) in enumerate(nodeChildren):
		nodes[i].child = nodes[c]
		nodes[c].parents.append(nodes[i])
	for r in range(1, rolls + 1):
		for n in nodes:
			n.update(diceProbs)
		for n in nodes:
			n.clear()
	return [n.lastProb for n in nodes]
import sys
data = sys.stdin.read().split('\n')
(n, m, k) = map(int, data[0].strip().split(' '))
children = list(map(lambda x: int(x) - 1, data[1].strip().split(' ')))
diceProbs = list(map(int, data[2].strip().split(' ')))
for p in DRW(children, diceProbs, k):
	print(p)