library_for_python
This documentation is automatically generated by online-judge-tools/verification-helper
Modulo_Matrix
(Modulo_Matrix/Modulo_Matrix.py)
- View this file on GitHub
- Last update: 2025-02-24 11:37:33+09:00
Outline
$m$ を素数とする. $\mathbb{Z}/m \mathbb{Z}$ を要素に持つ行列に対する様々な計算を行うクラス, 関数.
ただし, たいていの場合は $m$ が素数であることを要求する. これにより, $\mathbb{Z}/m \mathbb{Z}$ は体になる.
以降, $p$ が素数の時には体 $\mathbb{Z}/p \mathbb{Z}$ を $\mathbb{F}_p$ と書く.
Theory
逆行列
行列に対する次の操作を行基本変形という.
- $i$ 行目と $j$ 行目を交換する $(i \neq j)$
- $i$ 行目を $\alpha$ 倍 $(\alpha \in \mathbb{F}_p^\times)$ する.
- $i$ 行目に $j$ 行目の $\beta$ 倍 $(\beta \in \mathbb{F}_p)$ を加える.
次のような行列を基本行列という.
- $P_{N,i,j}:=(\bm{e}1, \dots, \bm{e}{i-1}, \bm{e}j, \bm{e}{i+1}, \dots, \bm{e}{j-1}, \bm{e}_i, \bm{e}{j+1}, \dots \bm{e}_N)$
- $Q_{N,i,\alpha}:=(\bm{e}1, \dots, \bm{e}{i-1}, \alpha \bm{e}i, \bm{e}{i+1}, \dots, \bm{e}_N)$
- $R_{N,i,j,\beta}:=(\bm{e}1, \dots, \bm{e}{j-1}, \bm{e}j+\beta \bm{e}_i, \bm{e}{j+1}, \dots \bm{e}_N)$
ここで, $A \in M_N(\mathbb{F}_p)$ に対して,
- $P_{N,i,j}A$ は $A$ の第 $i$ 行と第 $j$ 行を入れ替えた行列
- $Q_{N,i,\alpha} A$ は $A$ の第 $i$ 行を $\alpha$ 倍した行列
- $R_{N,i,j,\beta} A$ は $A$ の第 $i$ 行に第 $j$ 行の $\beta$ 倍を加えた行列
になり, 左から基本行列を掛けることと, 行基本変形を行うことが1対1に対応する.
ここで, 任意の $A \in M_N(\mathbb{F}_p)$ に対して, 基本行列 $S_1, \dots, S_L$ と $0 \leq K \leq N$ が存在して,
\[S_L S_{L-1} \dots S_1 A=(\bm{e}_1, \dots, \bm{e}_K, \bm{0}, \dots \bm{0})\]となる. ここで, $K$ については $S_1, \dots, S_L$ の取り方に依らず, $A$ のみによって定まる. よって, この $K$ のことを行列 $A$ の Rank (階数) といい, $\operatorname{rank} A$ と表す.
そして, $A \in M_N(\mathbb{F}_p)$ に対して, 以下は同値になる.
- $A \in M_N(\mathbb{F}_p)^\times$
- $\operatorname{rank} A=N$
つまり, $A \in M_N(\mathbb{F}_p)^\times$ のとき,
\[S_L S_{L-1} \dots S_1 A=(\bm{e}_1, \dots, \bm{e}_N)=I_N\]である. $S:=L S{L-1} \dots S_1$ とすると, $SA=I_N$ となる. このとき, $A^{-1}=S$ であることも導ける.
行列式
行列式 $\det: M_N(\mathbb{F}_p) \to \mathbb{F}_p$ を以下で定義する. ただし,$\mathfrak{S}_N$ で $N$ 次対称群を表すとする.
\[\det A:=\sum_{\sigma \in \mathfrak{S}_N} \operatorname{sgn} \sigma \prod_{i=1}^N A_{i, \sigma(i)}\]この行列式は次のような性質を満たす.
- 各行に対する多重線形性
- 交代性
- $\det I_N=1$
実はこの3条件を満たすような写像は行列式のみである.
また, 行列式は積に関して準同型を成す. つまり, $A,B \in M_N(\mathbb{F}_p)$ に対して,
\[\det (AB)=(\det A)(\det B)\]となる.
このとき, 逆行列の求め方と同様にして, 行基本行列 $S_1, \dots, S_L$ と上三角行列が存在して,
\[S_1 S_2 \dots S_K A=U\]となる.
そして,
- $\det P_{i,j}=-1$
- $\det Q_{i,\alpha}=\alpha$
- $\det R_{i,j,\beta}=1$
- $\det U$ は $N$ 個の対角成分の積
である.
これによって,
\[\det A=\dfrac{\det U}{(\det S_1) \dots (\det S_K)}\]として求められることができる. 計算量は $O(N^3)$ である.
Verified with
- test_verify/yosupo_library_checker/Matrix/Adjugate_Matrix.test.py
- test_verify/yosupo_library_checker/Matrix/Determinant.test.py
- test_verify/yosupo_library_checker/Matrix/Inverse.test.py
- test_verify/yosupo_library_checker/Matrix/Matrix_Rank.test.py
- test_verify/yosupo_library_checker/Matrix/Power.test.py
- test_verify/yosupo_library_checker/Matrix/Product.test.py
Code
from copy import deepcopy
class SingularMatrixError(Exception):
def __str__(self):
return "非正則行列の逆行列を求めようとしました."
class Modulo_Matrix():
__slots__ = ("ele", "__row", "__col")
# property
@property
def row(self):
return self.__row
@property
def col(self):
return self.__col
@property
def size(self):
return (self.row, self.col)
#入力
def __init__(self, M: list[list[int]]):
""" 行列 M を生成する.
※ Mod: 法はグローバル変数から指定
Args:
M (list[int]): 行列
"""
self.ele=[[x%Mod for x in X] for X in M]
self.__row = len(M)
self.__col = len(M[0]) if self.row > 0 else 0
#出力
def __str__(self):
return "["+"\n".join(map(str,self.ele))+"]"
def __repr__(self):
return str(self)
# 零行列, 単位行列
@classmethod
def Zero_Matrix(cls, row: int, col: int) -> "Modulo_Matrix":
""" row 行 col 列のゼロ行列を生成する.
Args:
row (int): 行
col (int): 列
Returns:
Modulo_Matrix: row 行 col 列のゼロ行列
"""
return Modulo_Matrix([[0] * col for _ in range(row)])
@classmethod
def Identity_Matrix(cls, N: int) -> "Modulo_Matrix":
""" N 次の単位行列を生成する.
Args:
N (int): 次数
Returns:
Modulo_Matrix: N 次単位行列
"""
return Modulo_Matrix([[1 if i==j else 0 for j in range(N)] for i in range(N)])
#+,-
def __pos__(self):
return self
def __neg__(self):
return self.__scale__(-1)
#加法
def __add__(self, other):
C = [None] * self.row
for i, (Ai, Bi) in enumerate(zip(self.ele, other.ele)):
C[i] = [Ai[j] + Bi[j] for j in range(self.col)]
return Modulo_Matrix(C)
def __iadd__(self,other):
M=self.ele; N=other.ele
for i in range(self.row):
Mi,Ni=M[i],N[i]
for j in range(self.col):
Mi[j]+=Ni[j]
Mi[j]%=Mod
return self
#減法
def __sub__(self,other):
C = [None] * self.row
for i, (Ai, Bi) in enumerate(zip(self.ele, other.ele)):
C[i] = [Ai[j] - Bi[j] for j in range(self.col)]
return Modulo_Matrix(C)
def __isub__(self,other):
M=self.ele; N=other.ele
for i in range(self.row):
Mi,Ni=M[i],N[i]
for j in range(self.col):
Mi[j]-=Ni[j]
Mi[j]%=Mod
return self
#乗法
def __mul__(self, other):
if isinstance(other, int):
return self.__scale__(other)
if not isinstance(other, Modulo_Matrix):
raise TypeError
assert self.col == other.row, f"左側の列と右側の行が一致しません (left: {self.col}, right:{other.row})."
A = self.ele; B = other.ele
C = [[0] * other.col for _ in range(self.row)]
for i, Ci in enumerate(C):
for k, a_ik in enumerate(A[i]):
for j, b_kj in enumerate(B[k]):
Ci[j] = (Ci[j] + a_ik * b_kj) % Mod
return Modulo_Matrix(C)
def __rmul__(self,other):
if isinstance(other,int):
return self.__scale__(other)
def inverse(self) -> "Modulo_Matrix":
""" 逆行列を求める
Raises:
SingularMatrixError: 非正則行列の逆行列を求めようとしたときに発生
Returns:
Modulo_Matrix: 逆行列
"""
inverse, _ = self.inverse_with_determinant()
if inverse is None:
raise SingularMatrixError()
return inverse
def inverse_with_determinant(self) -> tuple["Modulo_Matrix", int] | tuple[None, int]:
""" self の逆行列と self の行列式を求める.
Returns:
tuple["Modulo_Matrix", int] | tuple[None, int]: (逆行列 (非正則の場合は None), 行列式)
"""
assert self.row == self.col,"正方行列ではありません."
M = self
N = M.row
R = [[1 if i == j else 0 for j in range(N)] for i in range(N)]
T = deepcopy(M.ele)
det = 1
for j in range(N):
if T[j][j] == 0:
for i in range(j+1,N):
if T[i][j]:
break
else:
return None, 0
T[j], T[i] = T[i], T[j]
R[j], R[i] = R[i], R[j]
det = -det % Mod
Tj, Rj = T[j] ,R[j]
inv = pow(Tj[j], -1, Mod)
det = (Tj[j] * det) % Mod
for k in range(N):
Tj[k] *=inv; Tj[k] %= Mod
Rj[k] *=inv; Rj[k] %= Mod
for i in range(N):
if i == j:
continue
c = T[i][j]
Ti, Ri = T[i], R[i]
for k in range(N):
Ti[k] -= Tj[k] * c; Ti[k] %= Mod
Ri[k] -= Rj[k] * c; Ri[k] %= Mod
for i in range(N):
det = (T[i][i] * det) % Mod
return Modulo_Matrix(R), det
#スカラー倍
def __scale__(self, r: int) -> "Modulo_Matrix":
""" r 倍する
Args:
r (int): スカラー倍
Returns:
Modulo_Matrix: r 倍
"""
r %= Mod
return Modulo_Matrix([[r * m_ij for m_ij in Mi] for Mi in self.ele])
#累乗
def __pow__(self, n):
assert self.row==self.col, "正方行列ではありません."
sgn = 1 if n >= 0 else -1
n = abs(n)
C = Modulo_Matrix.Identity_Matrix(self.row)
tmp = self
while n:
if n & 1:
C = C * tmp
tmp = tmp * tmp
n >>= 1
return C if sgn == 1 else C.inverse()
#等号
def __eq__(self,other):
return self.ele==other.ele
#不等号
def __neq__(self,other):
return not(self==other)
#転置
def transpose(self) -> "Modulo_Matrix":
""" 転置行列を求める.
Returns:
Modulo_Matrix: 転置行列
"""
return Modulo_Matrix(list(map(list,zip(*self.ele))))
#行基本変形
def row_reduce(self) -> "Modulo_Matrix":
""" 行基本変形をできるだけ施した後の行列を求める.
Returns:
Modulo_Matrix: 行基本変形をできるだけ施した後の行列
"""
(row, col) = self.size
T = deepcopy(self.ele)
I = 0
for J in range(col):
if T[I][J] == 0:
for i in range(I + 1, row):
if T[i][J] != 0:
T[i], T[I] = T[I], T[i]
break
else:
continue
u = T[I][J]
u_inv = pow(u, -1, Mod)
for j in range(col):
T[I][j] *= u_inv
T[I][j] %= Mod
for i in range(row):
if i == I:
continue
v = T[i][J]
for j in range(col):
T[i][j] -= v * T[I][j]
T[i][j] %= Mod
I += 1
if I == row:
break
return Modulo_Matrix(T)
#列基本変形
def column_reduce(self) -> "Modulo_Matrix":
""" 列基本変形をできるだけ施した後の行列を求める.
Returns:
Modulo_Matrix: 列基本変形をできるだけ施した後の行列
"""
(row, col) = self.size
T = deepcopy(self.ele)
J = 0
for I in range(row):
if T[I][J] ==0 :
for j in range(J + 1, col):
if T[I][j] != 0:
for k in range(row):
T[k][j], T[k][J] = T[k][J], T[k][j]
break
else:
continue
u = T[I][J]
u_inv = pow(u, -1, Mod)
for i in range(row):
T[i][J] *= u_inv
T[i][J] %= Mod
for j in range(col):
if j != J:
v = T[I][j]
for i in range(row):
T[i][j] -= v * T[i][J]
T[i][j] %= Mod
J += 1
if J == col:
break
return Modulo_Matrix(T)
#行列の階数
def rank(self) -> int:
""" 行列のランクを求める
Returns:
int: ランク
"""
row_reduced = self.row_reduce()
(row, col) = row_reduced.size
rnk = 0
for i in range(row):
Ti = row_reduced.ele[i]
if any(Ti[j] for j in range(col)):
rnk += 1
return rnk
# 単射 ?
def is_injection(self) -> bool:
""" 行列が表す線形写像は単射?
Returns:
bool: 単射 ?
"""
return self.rank() == self.col
# 全射 ?
def is_surjective(self) -> bool:
""" 行列が表す線形写像は全射?
Returns:
bool: 全射 ?
"""
return self.rank() == self.row
# 全単射 ?
def is_bijection(self) -> bool:
""" 行列が表す線形写像は全単射?
Returns:
bool: 全単射 ?
"""
return self.col == self.row == self.rank()
#行の結合
def row_union(self,other):
return Modulo_Matrix(self.ele+other.ele)
#列の結合
def column_union(self,other):
E=[]
for i in range(self.row):
E.append(self.ele[i]+other.ele[i])
return Modulo_Matrix(E)
def __getitem__(self,index):
if isinstance(index, int):
return self.ele[index]
else:
return self.ele[index[0]][index[1]]
def __setitem__(self,index,val):
assert isinstance(index,tuple) and len(index)==2
self.ele[index[0]][index[1]]=val
#=================================================
#正方行列?
def Is_Square(M: Modulo_Matrix) -> bool:
""" M は正方行列か?
Args:
M (Modulo_Matrix): 行列
Returns:
bool: 正方行列 ?
"""
return M.row == M.col
#対角行列
def Diagonal_Matrix(D: list[int]) -> Modulo_Matrix:
""" D の第 i 成分が (i, i) 成分になる対角行列を生成する.
Args:
D (list[int]): 対角成分のリスト
Returns:
Modulo_Matrix: 対角行列
"""
N=len(D)
return Modulo_Matrix([[D[i] if i==j else 0 for j in range(N)] for i in range(N)])
#行列の直和
def Direct_Sum(*A):
""" A=[A_0, A_1, ..., A_{N-1}] に対する直和行列を求める.
"""
r=c=0
for a in A:
r+=a.row
c+=a.col
M=[[0]*c for _ in range(r)]
x=y=0
for p in range(len(A)):
a=A[p]
for i in range(a.row):
b=A[p].ele[i]
m=M[x+i]
for j in range(a.col):
m[y+j]=b[j]
x+=a.row; y+=a.col
return Modulo_Matrix(M)
#クロネッカー積
def Kronecker_Product(*X):
A=[[1]]
for B in X:
A=[[A[i//B.row][j//B.col]*B[i%B.row][j%B.col]%Mod for j in range(len(A[0])*B.col)] for i in range(len(A)*B.row)]
return Modulo_Matrix(A)
#クロネッカー和
def Kronecker_Sum(*X):
A=Modulo_Matrix([[0]])
for B in X:
A=Kronecker_Product(A, Modulo_Matrix.Identity_Matrix(B.row))+Kronecker_Product(Modulo_Matrix.Identity_Matrix(A.row),B)
return A
#跡
def Trace(M: Modulo_Matrix) -> int:
""" 正方行列 M の跡 (対角成分の和) を求める.
Args:
M (Modulo_Matrix): 正方行列
Returns:
int: 跡
"""
assert Is_Square(M)
return sum(M.ele[i][i] for i in range(M.row)) % Mod
def Determinant(M: Modulo_Matrix) -> int:
""" 正方行列 M の行列式 (素数 mod) を求める.
Args:
M (Modulo_Matrix): 正方行列
Returns:
int: 行列式 (mod 素数)
"""
assert Is_Square(M)
N=M.row
T=deepcopy(M.ele)
det=1
for j in range(N):
if T[j][j]==0:
for i in range(j+1,N):
if T[i][j]:
break
else:
return 0
T[j],T[i]=T[i],T[j]
det=-det
Tj=T[j]
inv=pow(Tj[j], -1, Mod)
for i in range(j+1,N):
Ti=T[i]
c=-inv*Ti[j]%Mod
for k in range(N):
Ti[k]+=c*Tj[k]
Ti[k]%=Mod
for i in range(N):
det*=T[i][i]
det%=Mod
return det
def Determinant_Arbitrary_Mod(A: Modulo_Matrix) -> int:
""" 正方行列 M の行列式 (任意 mod) を求める.
Args:
M (Modulo_Matrix): 正方行列
Returns:
int: 行列式 (mod 任意)
"""
N=A.row
A=deepcopy(A.ele)
det=1
for i in range(N):
Ai=A[i]
for j in range(i+1, N):
Aj=A[j]
while Aj[i]:
alpha=Ai[i]//Aj[i]
if alpha:
for k in range(i, N):
Ai[k]-=alpha*Aj[k]
Ai[k]%=Mod
A[i], A[j]=A[j], A[i]
Ai=A[i]; Aj=A[j]
det*=-1
det*=Ai[i]
det%=Mod
if det==0:
break
return det
def Characteristic_Polynomial(M: Modulo_Matrix) -> list[int]:
""" 正方行列 M の固有多項式を sum(P[i] X^i) としたとき, P を求める.
Args:
M (Modulo_Matrix): 正方行列
Returns:
list[int]: 固有多項式を sum(P[i] X^i) としたときの P を求める.
"""
T=deepcopy(M.ele)
N=M.row
for j in range(N-2):
for i in range(j+1, N):
if T[i][j]:
break
else:
continue
T[j+1],T[i]=T[i],T[j+1]
for k in range(N):
T[k][j+1],T[k][i]=T[k][i],T[k][j+1]
if T[j+1][j]:
Tjj=T[j+1]
inv=pow(Tjj[j], -1, Mod)
for i in range(j+2, N):
Ti=T[i]
c=inv*Ti[j]%Mod
for k in range(j,N):
Ti[k]-=c*Tjj[k]
Ti[k]%=Mod
for k in range(N):
T[k][j+1]+=c*T[k][i]
T[k][j+1]%=Mod
dp=[[0]*(i+1) for i in range(N+1)]; dp[0][0]=1
for i in range(N):
P=dp[i+1]
for k in range(i+1):
P[k+1]=dp[i][k]
for k in range(i+1):
P[k]+=T[i][i]*dp[i][k]
P[k]%=Mod
p=1
for j in range(i-1,-1,-1):
p*=-T[j+1][j]; p%=Mod
c=p*T[j][i]%Mod
for k in range(j+1):
P[k]+=c*dp[j][k]
P[k]%=Mod
P=dp[-1]
for i in range(N+1):
if i%2:
P[~i]*=-1; P[~i]%=Mod
return P
def Adjugate_Matrix(A: Modulo_Matrix) -> Modulo_Matrix:
""" 正方行列 A の余因子行列 adj A := ((-1)^(i+j) det A_{i,j}) を求める.
Args:
A (Modulo_Matrix): 正方行列
Returns:
Modulo_Matrix: 余因子行列
"""
from random import randint
N = A.row
A_ext = [[0] * (N + 1) for _ in range(N + 1)]
for i in range(N):
for j in range(N):
A_ext[i][j] = A[i][j]
for i in range(N):
A_ext[i][N] = A_ext[N][i] = randint(0, Mod - 1)
A_ext_inv, det = Modulo_Matrix(A_ext).inverse_with_determinant()
if A_ext_inv is None:
return Modulo_Matrix.Zero_Matrix(N, N)
adj = [[det * ((A_ext_inv[N][N] * A_ext_inv[i][j] - A_ext_inv[i][N] * A_ext_inv[N][j]) % Mod) for j in range(N)] for i in range(N)]
return Modulo_Matrix(adj)
#===
Mod=998244353Traceback (most recent call last):
File "/opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages/onlinejudge_verify/documentation/build.py", line 71, in _render_source_code_stat
bundled_code = language.bundle(stat.path, basedir=basedir, options={'include_paths': [basedir]}).decode()
~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages/onlinejudge_verify/languages/python.py", line 96, in bundle
raise NotImplementedError
NotImplementedError