from copy import deepcopy from random import randint from math import pi, sin, cos from __pypy__.builders import StringBuilder import sys from os import read as os_read, write as os_write from atexit import register as atexist_register from typing import Generic, Iterator, List, Tuple, Dict, Iterable, Sequence, Callable, Union, Optional, TypeVar T = TypeVar('T') Graph = List[List[int]] Poly = List[int] Vector = List[int] Matrix = List[List[int]] Func10 = Callable[[int], None] Func20 = Callable[[int, int], None] Func11 = Callable[[int], int] Func21 = Callable[[int, int], int] Func31 = Callable[[int, int, int], int] class Fastio: ibuf = bytes() pil = pir = 0 sb = StringBuilder() def load(self): self.ibuf = self.ibuf[self.pil:] self.ibuf += os_read(0, 131072) self.pil = 0; self.pir = len(self.ibuf) def flush_atexit(self): os_write(1, self.sb.build().encode()) def flush(self): os_write(1, self.sb.build().encode()) self.sb = StringBuilder() def fastin(self): if self.pir - self.pil < 64: self.load() minus = x = 0 while self.ibuf[self.pil] < 45: self.pil += 1 if self.ibuf[self.pil] == 45: minus = 1; self.pil += 1 while self.ibuf[self.pil] >= 48: x = x * 10 + (self.ibuf[self.pil] & 15) self.pil += 1 if minus: return -x return x def fastin_string(self): if self.pir - self.pil < 64: self.load() while self.ibuf[self.pil] <= 32: self.pil += 1 res = bytearray() while self.ibuf[self.pil] > 32: if self.pir - self.pil < 64: self.load() res.append(self.ibuf[self.pil]) self.pil += 1 return res def fastout(self, x): self.sb.append(str(x)) def fastoutln(self, x): self.sb.append(str(x)); self.sb.append('\n') fastio = Fastio() rd = fastio.fastin; rds = fastio.fastin_string; wt = fastio.fastout; wtn = fastio.fastoutln; flush = fastio.flush atexist_register(fastio.flush_atexit) sys.stdin = None; sys.stdout = None def rdl(n): return [rd() for _ in range(n)] def wtnl(l): wtn(' '.join(map(str, l))) def wtn_yes(): wtn("Yes") def wtn_no(): wtn("No") # https://judge.yosupo.jp/submission/46633 # cooley_tukey ? # maybe https://nyaannyaan.github.io/library/ntt/complex-fft.hpp class CooleyTukey: isbuilt = 0 def __init__(self, k: int=20) -> None: self.wr = [0] * (1 << 20) self.wi = [0] * (1 << 20) self.baser = [0] * 20 self.basei = [0] * 20 self.setw(k) self.isbuilt = 1 @staticmethod def mul(xr: float, xi: float, yr: float, yi: float) -> tuple: return xr * yr - xi * yi, xr * yi + yr * xi def genw(self, i: int, b: int, zr: float, zi: float) -> None: if b == -1: self.wr[i] = zr self.wi[i] = zi else: self.genw(i, b - 1, zr, zi) wr, wi = self.baser[b], self.basei[b] self.genw(i | (1 << b), b - 1, zr * wr - zi * wi, zr * wi + zi * wr) def setw(self, k: int) -> None: k -= 1 arg = pi / (1 << k) i = 0 j = 1 << (k - 1) while j: self.baser[i] = cos(arg * j) self.basei[i] = sin(arg * j) i += 1 j >>= 1 self.genw(0, k - 1, 1, 0) def fft(self, ar: List[float], ai: List[float], k: int) -> None: if k == 0: return if k == 1: ar[0], ar[1] = ar[0] + ar[1], ar[0] - ar[1] ai[0], ai[1] = ai[0] + ai[1], ai[0] - ai[1] return if k & 1: v = 1 << (k - 1) for j in range(v): ar[j], ar[j + v] = ar[j] + ar[j + v], ar[j] - ar[j + v] ai[j], ai[j + v] = ai[j] + ai[j + v], ai[j] - ai[j + v] u = 1 << (k & 1) v = 1 << (k - 2 - (k & 1)) wr1, wi1 = self.wr[1], self.wi[1] while v: for j0 in range(v): t0r = ar[j0]; t0i = ai[j0] t1r = ar[j0 + v]; t1i = ai[j0 + v] t2r = ar[j0 + v * 2]; t2i = ai[j0 + v * 2] t3r = ar[j0 + v * 3]; t3i = ai[j0 + v * 3] t1m3r, t1m3i = self.mul(t1r - t3r, t1i - t3i, wr1, wi1) ar[j0] = (t0r + t2r) + (t1r + t3r); ai[j0] = (t0i + t2i) + (t1i + t3i) ar[j0 + v] = (t0r + t2r) - (t1r + t3r); ai[j0 + v] = (t0i + t2i) - (t1i + t3i) ar[j0 + v * 2] = (t0r - t2r) + t1m3r; ai[j0 + v * 2] = (t0i - t2i) + t1m3i ar[j0 + v * 3] = (t0r - t2r) - t1m3r; ai[j0 + v * 3] = (t0i - t2i) - t1m3i for jh in range(1, u): p = jh * v * 4 Wr = self.wr[jh]; Wi = self.wi[jh] Xr = self.wr[jh << 1]; Xi = self.wi[jh << 1] WXr, WXi = self.mul(Wr, Wi, Xr, Xi) for offset in range(v): t0r = ar[p + offset]; t0i = ai[p + offset] t1r, t1i = self.mul(ar[p + offset + v], ai[p + offset + v], Xr, Xi) t2r, t2i = self.mul(ar[p + offset + v * 2], ai[p + offset + v * 2], Wr, Wi) t3r, t3i = self.mul(ar[p + offset + v * 3], ai[p + offset + v * 3], WXr, WXi) t1m3r, t1m3i = self.mul(t1r - t3r, t1i - t3i, wr1, wi1) ar[p + offset] = (t0r + t2r) + (t1r + t3r); ai[p + offset] = (t0i + t2i) + (t1i + t3i) ar[p + offset + v] = (t0r + t2r) - (t1r + t3r); ai[p + offset + v] = (t0i + t2i) - (t1i + t3i) ar[p + offset + v * 2] = (t0r - t2r) + t1m3r; ai[p + offset + v * 2] = (t0i - t2i) + t1m3i ar[p + offset + v * 3] = (t0r - t2r) - t1m3r; ai[p + offset + v * 3] = (t0i - t2i) - t1m3i u <<= 2 v >>= 2 def ifft(self, ar: List[float], ai: List[float], k: int) -> None: if k == 0: return if k == 1: ar[0], ar[1] = ar[0] + ar[1], ar[0] - ar[1] ai[0], ai[1] = ai[0] + ai[1], ai[0] - ai[1] return u = 1 << (k - 2) v = 1 wr1, mwi1 = self.wr[1], -self.wi[1] while u: for j0 in range(v): t0r = ar[j0]; t0i = ai[j0] t1r = ar[j0 + v]; t1i = ai[j0 + v] t2r = ar[j0 + v * 2]; t2i = ai[j0 + v * 2] t3r = ar[j0 + v * 3]; t3i = ai[j0 + v * 3] t2m3r, t2m3i = self.mul(t2r - t3r, t2i - t3i, wr1, mwi1) ar[j0] = (t0r + t1r) + (t2r + t3r); ai[j0] = (t0i + t1i) + (t2i + t3i) ar[j0 + v * 2] = (t0r + t1r) - (t2r + t3r); ai[j0 + v * 2] = (t0i + t1i) - (t2i + t3i) ar[j0 + v] = (t0r - t1r) + t2m3r; ai[j0 + v] = (t0i - t1i) + t2m3i ar[j0 + v * 3] = (t0r - t1r) - t2m3r; ai[j0 + v * 3] = (t0i - t1i) - t2m3i for jh in range(1, u): p = jh * v * 4 Wr = self.wr[jh]; Wi = -self.wi[jh] Xr = self.wr[(jh << 1) + 0]; Xi = -self.wi[(jh << 1) + 0] Yr = self.wr[(jh << 1) + 1]; Yi = -self.wi[(jh << 1) + 1] for offset in range(v): t0r = ar[p + offset]; t0i = ai[p + offset] t1r = ar[p + offset + v]; t1i = ai[p + offset + v] t2r = ar[p + offset + v * 2]; t2i = ai[p + offset + v * 2] t3r = ar[p + offset + v * 3]; t3i = ai[p + offset + v * 3] t0m1r, t0m1i = self.mul(t0r - t1r, t0i - t1i, Xr, Xi) t2m3r, t2m3i = self.mul(t2r - t3r, t2i - t3i, Yr, Yi) ar[p + offset] = (t0r + t1r) + (t2r + t3r); ai[p + offset] = (t0i + t1i) + (t2i + t3i) ar[p + offset + v] = t0m1r + t2m3r; ai[p + offset + v] = t0m1i + t2m3i ar[p + offset + v * 2], ai[p + offset + v * 2] = self.mul((t0r + t1r) - (t2r + t3r), (t0i + t1i) - (t2i + t3i), Wr, Wi) ar[p + offset + v * 3], ai[p + offset + v * 3] = self.mul(t0m1r - t2m3r, t0m1i - t2m3i, Wr, Wi) u >>= 2 v <<= 2 if k & 1: u = 1 << (k - 1) for j in range(u): ar[j], ar[j + u] = ar[j] + ar[j + u], ar[j] - ar[j + u] ai[j], ai[j + u] = ai[j] + ai[j + u], ai[j] - ai[j + u] def fft_real(self, ALr: List[float], ALi: List[float], AHr: List[float], AHi: List[float], k: int) -> None: self.fft(ALr, ALi, k) AHr[0] = ALi[0] * 2; AHi[0] = 0 ALr[0] = ALr[0] * 2; ALi[0] = 0 AHr[1] = ALi[1] * 2; AHi[1] = 0 ALr[1] = ALr[1] * 2; ALi[1] = 0 i = 2; y = 2 while y < 1 << k: while i < y << 1: j = i ^ (y - 1) AHr[i] = ALi[j] + ALi[i]; AHi[i] = ALr[j] - ALr[i] ALr[i] = ALr[j] + ALr[i]; ALi[i] = -ALi[j] + ALi[i] AHr[j] = AHr[i]; AHi[j] = -AHi[i] ALr[j] = ALr[i]; ALi[j] = -ALi[i] i += 2 y <<= 1 def karatsuba(self, a: Vector, b: Vector, mod: int) -> Vector: B = 32000 bbmod = B * B % mod l = len(a) + len(b) - 1 k = max(2, l.bit_length()) M = 1 << k if not self.isbuilt: self.setw(k) alr = [float()] * M ali = [float()] * M ahr = [float()] * M ahi = [float()] * M for i, x in enumerate(a): quo, rem = divmod(x, B) alr[i], ali[i] = float(rem), float(quo) self.fft_real(alr, ali, ahr, ahi, k) blr = [float()] * M bli = [float()] * M bhi = [float()] * M bhr = [float()] * M for i, x in enumerate(b): quo, rem = divmod(x, B) blr[i], bli[i] = float(rem), float(quo) self.fft_real(blr, bli, bhr, bhi, k) for i in range(M): alri = alr[i]; alii = ali[i] mahii = -ahi[i]; ahri = ahr[i] tmp1r, tmp1i = self.mul(alri, alii, blr[i], bli[i]) tmp2r, tmp2i = self.mul(mahii, ahri, bhr[i], bhi[i]) tmp3r, tmp3i = self.mul(alri, alii, bhr[i], bhi[i]) tmp4r, tmp4i = self.mul(mahii, ahri, blr[i], bli[i]) blr[i] = tmp1r + tmp2r; bli[i] = tmp1i + tmp2i bhr[i] = tmp3r + tmp4r; bhi[i] = tmp3i + tmp4i self.ifft(blr, bli, k) self.ifft(bhr, bhi, k) u = [0] * l im = float(1 / (4 * M)) for i in range(l): x1 = round(blr[i] * im) % mod x2 = (round(bhr[i] * im) + round(bhi[i] * im)) % mod * B % mod x3 = round(bli[i] * im) % mod * bbmod % mod x = x1 + x2 + x3 if x >= mod: x -= mod if x >= mod: x -= mod u[i] = x return u def karatsuba_pow2(self, a: Vector, mod: int) -> Vector: B = 32000 l = len(a) * 2 - 1 k = 2; M = 4 while M < l: M <<= 1 k += 1 self.setw(k) alr = [float()] * M ali = [float()] * M ahr = [float()] * M ahi = [float()] * M for i, x in enumerate(a): quo, rem = divmod(x, B) alr[i], ali[i] = float(rem), float(quo) self.fft_real(alr, ali, ahr, ahi, k) for i in range(M): tmp1r = alr[i]; tmp1i = ali[i] tmp2r = -ahi[i]; tmp2i = ahr[i] tmp3r = tmp1r; tmp3i = tmp1i tmp4r = tmp2r; tmp4i = tmp2i tmp1r, tmp1i = self.mul(tmp1r, tmp1i, alr[i], ali[i]) tmp2r, tmp2i = self.mul(tmp2r, tmp2i, ahr[i], ahi[i]) tmp3r, tmp3i = self.mul(tmp3r, tmp3i, ahr[i], ahi[i]) tmp4r, tmp4i = self.mul(tmp4r, tmp4i, alr[i], ali[i]) alr[i] = tmp1r + tmp2r; ali[i] = tmp1i + tmp2i ahr[i] = tmp3r + tmp4r; ahi[i] = tmp3i + tmp4i self.ifft(alr, ali, k) self.ifft(ahr, ahi, k) u = [0] * l im = float(1 / (4 * M)) for i in range(l): alr[i] *= im; ali[i] *= im ahr[i] *= im; ahi[i] *= im x1 = round(alr[i]) % mod x2 = (round(ahr[i]) + round(ahi[i])) % mod * B % mod x3 = round(ali[i]) % mod * (B * B % mod) % mod x1 += x2 if x1 >= mod: x1 -= mod x1 += x3 if x1 >= mod: x1 -= mod u[i] = x1 return u def modinv(a: int, m: int) -> int: '''return x s.t. x == a^(-1) (mod m)''' b = m; u = 1; v = 0 while b: t = a // b a, b = b, a - t * b u, v = v, u - t * v u %= m return u # https://nyaannyaan.github.io/library/fps/berlekamp-massey.hpp def berlekamp_massey(s: Vector, mod: int) -> Vector: N = len(s) b = [1] c = [1] y = 1 for ed in range(1, N + 1): l = len(c) m = len(b) x = 0 for i, a in enumerate(c): x += a * s[ed - l + i] x %= mod b.append(0) m += 1 if x == 0: continue freq = x * modinv(y, mod) % mod if l < m: tmp = c[:] c[:0] = [0] * (m - l) for i in range(m): c[m - 1 - i] = (c[m - 1 - i] - freq * b[m - 1 - i]) % mod b = tmp y = x else: for i in range(m): c[l - 1 - i] = (c[l - 1 - i] - freq * b[m - 1 - i]) % mod c.reverse() return c # https://nyaannyaan.github.io/library/fps/formal-power-series.hpp class FPS: ntt = None def __init__(self, mod: int) -> None: self.mod = mod @staticmethod def shrink(a: Vector) -> None: '''remove high degree coef == 0''' while a and not a[-1]: a.pop() @staticmethod def resize(a: Poly, length: int, val: int=0) -> None: a[length:] = [] a[len(a):] = [val] * (length - len(a)) def add(self, l: Poly, r: Union[Poly, int]) -> Poly: '''l += r''' if type(r) is int: res = l[:] res[0] = (res[0] + r) % self.mod return res mod = self.mod if type(r) is list: if len(l) < len(r): res = r[::] for i, x in enumerate(l): res[i] += x else: res = l[::] for i, x in enumerate(r): res[i] += x return [x % mod for x in res] raise TypeError() def sub(self, l: Poly, r: Union[Poly, int]) -> Poly: '''l -= r''' if type(r) is int: return self.add(l, -r) if type(r) is list: return self.add(l, self.neg(r)) raise TypeError() def neg(self, a: Poly) -> Poly: '''a *= -1''' mod = self.mod return [mod - x if x else 0 for x in a] def matmul(self, l: Poly, r: Poly) -> Poly: 'not verified' mod = self.mod return [x * r[i] % mod for i, x in enumerate(l)] def div(self, l: Poly, r: Poly) -> Poly: '''return: quo s.t. l = r*quo + rem''' if len(l) < len(r): return [] n = len(l) - len(r) + 1 if len(r) > 64: return self.mul(l[::-1][:n], self.inv(r[::-1], n))[:n][::-1] mod = self.mod f, g = l[::], r[::] cnt = 0 while g and not g[-1]: g.pop() cnt += 1 coef = modinv(g[-1], mod) g = self.mul(g, coef) deg = len(f) - len(g) + 1 gs = len(g) quo = [0] * deg for i in range(deg)[::-1]: quo[i] = x = f[i + gs - 1] % mod for j, y in enumerate(g): f[i + j] -= x * y return self.mul(quo, coef) + [0] * cnt def modulo(self, l: Poly, r: Poly) -> Poly: '''return: rem s.t. l = r*quo + rem''' res = self.sub(l, self.mul(self.div(l, r), r)) self.shrink(res) return res def divmod(self, l: Poly, r: Poly) -> Tuple[Poly, Poly]: '''return: quo, rem s.t. l = r*quo + rem''' quo = self.div(l, r) rem = self.sub(l, self.mul(quo, r)) self.shrink(rem) return quo, rem def eval(self, a: Poly, x: int) -> int: mod = self.mod r = 0; w = 1 for v in a: r += w * v % mod w = w * x % mod return r % mod def pow(self, f: Poly, k: int, deg=-1) -> Poly: '''return: g s.t. g == f**k (mod x**deg)''' n = len(f) if deg == -1: deg = n if k == 0: if not deg: return [] ret = [0] * deg ret[0] = 1 return ret mod = self.mod for i, x in enumerate(f): if x: rev = modinv(x, mod) ret = self.mul(self.exp(self.mul(self.log(self.mul(f, rev)[i:], deg), k), deg), pow(x, k, mod)) ret[:0] = [0] * (i * k) if len(ret) < deg: self.resize(ret, deg) return ret return ret[:deg] if (i + 1) * k >= deg: break return [0] * deg def log(self, f: Poly, deg=-1) -> Poly: '''return: g s.t. g == log(f) (mod x**deg)''' # assert(a[0] == 1) if deg == -1: deg = len(f) return self.integral(self.mul(self.diff(f), self.inv(f, deg))[:deg - 1]) def integral(self, a: Poly) -> Poly: mod = self.mod n = len(a) res = [0] * (n + 1) if n: res[1] = 1 for i in range(2, n + 1): j, k = divmod(mod, i) res[i] = (-res[k] * j) % mod for i, x in enumerate(a): res[i + 1] = res[i + 1] * x % mod return res def diff(self, f: Poly) -> Poly: '''return: dfdx''' mod = self.mod return [i * x % mod for i, x in enumerate(f) if i] def mul(self, l: list, r: Union[list, int]) -> list: raise AttributeError("'mul' method is virtual function. please override") def mul2(self, l: list) -> list: raise AttributeError("'mul2' method is virtual function. please override") def inv(self, a: list, deg: int=-1) -> list: raise AttributeError("'inv' method is virtual function. please override") def exp(self, a: list, deg: int=-1) -> list: raise AttributeError("'exp' method is virtual function. please override") # arbitrary_nttだとTLEするのでCooleyTukeyを使いたい # 2倍くらい速いっぽい # https://nyaannyaan.github.io/library/fps/arbitrary-fps.hpp def set_fft(self: FPS) -> None: self.ntt = CooleyTukey() FPS.set_fft = set_fft def mul(self: FPS, l: Poly, r: Union[Poly, int]) -> Poly: ''' if r is int: l *= r if r is Polynomial: convolve l and r ''' mod = self.mod if type(r) is int: return [x * r % mod for x in l] if type(r) is list: if not l or not r: return [] if self.ntt is None: self.set_fft() return self.ntt.karatsuba(l, r, mod) raise TypeError() FPS.mul = mul def mul2(self: FPS, l: Poly) -> Poly: ''' convolve l and l ''' mod = self.mod if self.ntt is None: self.set_fft() return self.ntt.karatsuba_pow2(l, mod) FPS.mul2 = mul2 def inv(self: FPS, a: Poly, deg: int=-1) -> Poly: '''return: g s.t. a*g == 1 (mod x**deg)''' # assert(self[0] != 0) if deg == -1: deg = len(a) mod = self.mod res = [pow(a[0], mod - 2, mod)] i = 1 while i < deg: i <<= 1 res = self.add(res, self.sub(res, self.mul(self.mul2(res), a[:i])[:i])) return res[:deg] FPS.inv = inv def exp(self: FPS, f: Poly, deg: int=-1) -> Poly: '''return: g s.t. log(g) == f (mod x ** deg)''' # assert(not self or self[0] == 0) if deg == -1: deg = len(f) ret = [1] i = 1 while i < deg: i <<= 1 ret = self.mul(ret, self.sub(self.add(f[:i], 1), self.log(ret, i)))[:i] return ret[:deg] FPS.exp = exp # https://nyaannyaan.github.io/library/fps/mod-pow.hpp def mod_pow(k: int, base: Poly, d: Poly, mod: int, fps: FPS=None) -> Poly: assert(d) if not fps: fps = FPS(mod) inv = fps.inv(d[::-1]) def quo(poly: Poly) -> Poly: if len(poly) < len(d): return [] n = len(poly) - len(d) + 1 return fps.mul(poly[:len(poly) - n - 1:-1], inv[:n])[n - 1::-1] res = [1] b = base[:] while k: if k & 1: res = fps.mul(res, b) res = fps.sub(res, fps.mul(quo(res), d)) FPS.shrink(res) b = fps.mul2(b) b = fps.sub(b, fps.mul(quo(b), d)) FPS.shrink(b) k >>= 1 # assert(len(b) + 1 <= len(d)) # assert(len(res) + 1 <= len(d)) return res # https://nyaannyaan.github.io/library/matrix/black-box-linear-algebra.hpp def inner_product(a: Poly, b: Poly, mod: int) -> int: res = 0 n = len(a) assert(n == len(b)) for i in range(n): res += a[i] * b[i] % mod return res % mod def random_poly(n: int, mod: int) -> Poly: return [randint(0, mod - 1) for _ in range(n)] class ModMatrix: def __init__(self, n: int, mod: int) -> None: self.mat = [[0] * n for _ in range(n)] self.mod = mod def add(self, i: int, j: int, x: int) -> None: self.mat[i][j] += x def __mul__(self, r: Poly) -> Poly: assert(len(self.mat) == len(r)) mod = self.mod return [sum(matij * r[j] % mod for j, matij in enumerate(mati)) % mod for mati in self.mat] def apply(self, i: int, r: int) -> None: mod = self.mod for j, matij in enumerate(self.mat[i]): self.mat[i][j] = matij * r % mod class SparseMatrix: def __init__(self, n: int, mod: int) -> None: self.mat: List[List[int]] = [[] for _ in range(n)] self.mod = mod def add(self, i: int, j: int, x: int) -> None: self.mat[i].append(j << 30 | x) def __mul__(self, r: Poly) -> Poly: assert(len(self.mat) == len(r)) mod = self.mod return [sum((jx & 0x3fffffff) * r[jx >> 30] % mod for jx in mati) % mod for mati in self.mat] def apply(self, i: int, r: int) -> None: mod = self.mod for idx, jx in enumerate(self.mat[i]): self.mat[i][idx] = (jx >> 30) << 30 | ((jx & 0x3fffffff) * r % mod) def vector_minpoly(b: List[Poly], mod: int) -> Poly: assert(b) n = len(b); m = len(b[0]) u = random_poly(m, mod) a = [0] * n for i, bi in enumerate(b): a[i] = inner_product(bi, u, mod) return berlekamp_massey(a, mod) def mat_minpoly(A: Union[ModMatrix, SparseMatrix]) -> Poly: n = len(A.mat) u = random_poly(n, A.mod) b: List[Poly] = [0] * (n << 1 | 1) for i in range(len(b)): b[i] = u u = A * u return vector_minpoly(b, A.mod) def fast_pow(A: Union[ModMatrix, SparseMatrix], b: Poly, k: int) -> Poly: n = len(b) mp = mat_minpoly(A) fps = FPS(A.mod) c = mod_pow(k, [0, 1], mp[::-1], A.mod, fps) res = [0] * n for ci in c: res = fps.add(res, fps.mul(b, ci)) b = A * b return res def fast_det(A: Union[ModMatrix, SparseMatrix]) -> int: n = len(A.mat) assert(n == len(A.mat)) mod = A.mod D = random_poly(n, mod) while 1: while any([not x for x in D]): D = random_poly(n, mod) AD = deepcopy(A) for i, d in enumerate(D): AD.apply(i, d) mp = mat_minpoly(AD) if mp[-1] == 0: return 0 if len(mp) != n + 1: continue det = -mp[-1] if n & 1 else mp[-1] Ddet = 1 for d in D: Ddet = Ddet * d % mod return det * modinv(Ddet, mod) % mod exit(1) # https://yukicoder.me/problems/no/1112 mod = 1_000_000_007 K, M, N = rd(), rd(), rd() m = ModMatrix(K * K, mod) for i in range(M): p, q, r = rd() - 1, rd() - 1, rd() - 1 m.add(p * K + q, q * K + r, 1) b = [0] * (K * K) for i in range(K): b[i * K] = 1 res = fast_pow(m, b, N - 2) wtn(sum(res[:K]))