結果
| 問題 |
No.1112 冥界の音楽
|
| コンテスト | |
| ユーザー |
 abUma
|
| 提出日時 | 2023-06-14 03:04:34 |
| 言語 | PyPy3 (7.3.15) |
| 結果 |
AC
|
| 実行時間 | 297 ms / 2,000 ms |
| コード長 | 23,643 bytes |
| コンパイル時間 | 423 ms |
| コンパイル使用メモリ | 82,280 KB |
| 実行使用メモリ | 112,784 KB |
| 最終ジャッジ日時 | 2024-06-22 15:08:19 |
| 合計ジャッジ時間 | 8,082 ms |
|
ジャッジサーバーID （参考情報） |
judge4 / judge5 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 3 |
| other | AC * 34 |
ソースコード
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]) % mod)