August 18, 2022 00:00

Introduction

$n \times n$ 정사각행렬 $A = (a _ {ij})$가 주어져 있을 때, $A$의 determinant $\det A$는 아래와 같이 계산할 수 있습니다.

오늘은 $\det A$ 와, 그를 일반화하는 특성다항식(Characteristic polynomial) $\phi _ {A}(x) = \det(x I - A) = x^{n} - \mathrm{tr} (A) x^{n-1} + \cdots + (-1)^{n-1} \det A$를 계산하는 방법을 알아봅니다. 보다 정확히는, $A$의 원소가 단순히 실수나 $\mod p$ 가 아닌 일반적인 경우를 포괄하는 방법을 알아봅니다.

독자가 Gaussian Elimination을 이해하고 있다고 가정합니다. 학부 대수학을 간단히 알면 이 글에서 사용하는 모든 알고리즘의 등장 맥락을 깊이 이해할 수 있으나, 모르는 단어나 cf)가 나오면 무시하셔도 이해에 지장이 없습니다.

Algebraic conditions on $A$

$A$의 원소가 될 수 있는 수들은 실수 $\mathbb{R}$, 복소수 $\mathbb{C}$, modulo $M$으로 나타낸 수들의 모임 $\mathbb{Z} _ {M}$, 정수 $\mathbb{Z}$, 정수 계수 다항식들의 모임 $\mathbb{Z}[x]$를 비롯하여 무궁무진하게 많습니다. 만약 원소가 복소수인 행렬, 실수인 행렬 등에 대해 $\det A$를 계산하는 알고리즘이 모두 다르다면 효율적이지 않을 것입니다.

따라서 우리는 $A$의 원소들이 속한 집합 $R$의 실체를 들여다보지 않고, $R$에서 “어떤 연산이 가능한가?” 라는 조건에 따라 알고리즘을 설계합니다. $\det A$의 계산식을 들여다보면, 최소한 덧셈과 곱셈에 대한 교환법칙은 성립해야 할 것 같습니다.

구체적으로, 다음 조건을 만족하는 수(원소)들의 집합 $R$을 Commutative Ring (가환환)이라고 합니다.

• 덧셈, 곱셈에 대해 닫혀 있다.
• 덧셈, 곱셈에 대해 항등원(0, 1)이 있고 교환법칙/결합법칙이 성립한다.
• 덧셈에 대한 역원 $-a$가 존재한다.
• 분배법칙 $a(b+c) = ab + ac$가 성립한다.

엄밀한 정의는 이와 다를 수 있으나, 최소한 $\det A$를 구하기 위해선 수들을 곱하고 더하고 빼야 하니 이것만큼은 포기할 수 없는 조건들입니다. 주목할 것은 곱셈에 대한 역원 $1 / a$의 존재를 가정하지 않았다는 것입니다. 즉, determinant를 계산하는 데 있어 곱셈에 대한 역원은 필수적인 조건이 아닙니다.

• $\mathbb{Z} _ {6}$의 경우 $2 \cdot 3 = 6 = 0$이므로, $2$의 곱셈에 대한 역원 $x$는 존재하지 않습니다. 그런 $x$가 존재한다면 $2x = 1$이어야 하는데, $3 = 3\cdot 2x = 6x = 0$이기 때문입니다.

어떤 Commutative Ring (이하 CRing) $R$에 대해 $A$의 원소들이 $R$에 속한다고 하고, 이 경우에 $A \in R^{n \times n}$이라고 씁시다. 이 때 $\det A \in R$이 성립하는 것도 알 수 있습니다.

$R$을 계수로 갖는 다항식들로 이루어진 CRing을 $R[x]$라고 쓰면, Matrix $xI - A$는 $R[x]^{n \times n}$이기도 합니다. 이 때 $\phi _ {A}(x) = \det(xI - A) \in R[x]$가 되겠네요.

앞으로 3개의 챕터로 나누어, 우리는 $R$에 어떤 조건이 주어졌을 때 $\det A$와 $\phi _ {A}(x)$를 계산하는 방법을 알아볼 것입니다. 각 챕터는

• $R$에서 나눗셈이 가능한 경우 (Gauss, Hessenberg)
• $R$의 두 원소를 몫과 나머지꼴로 표현할 수 있는 경우
• $R$에 아무 조건도 주어져 있지 않은 경우 (Namajan)

에 대해 다룹니다.

cf) Determinant의 정의에서 곱셈의 교환법칙은 필수적이지 않다고 합니다. Non-commutative determinant를 계산하는 것은 어렵다는 것이 많이 알려져 있으니 (Nisan91, Arvind09), 이것까진 고려하지 않기로 합니다.

When $R$ is a field

CRing $R$의 모든 원소 $x \neq 0$에 대해 역원 $1 / x$가 있으면, $R$을 체 (Field)라고 합니다.

• CRing’s which are fields: $\mathbb{R}$, $\mathbb{C}$, 소수 $p$에 대해 $\mathbb{Z} _ {p}$, Galois field $\mathrm{GF}(p^n)$ 등
• CRing’s which are not fields: $\mathbb{Z}$, Polynomial ring $\mathbb{C}[x]$, 합성수 $n$에 대해 $\mathbb{Z} _ {n}$ 등

Determinant

Field에서 덧셈, 뺄셈, 나눗셈, 곱셈을 하는 시간이 모두 상수 시간이라고 가정하면, Gaussian Elimination을 이용해 $O(n^3)$ 시간에 손쉽게 Determinant를 계산할 수 있습니다. $n \times n$ 행렬곱을 하는 시간 $O(n^\omega) \simeq O(n^{2.38})$ 시간에도 구하는 방법이 알려져 있으나 이 글에서 $O(n^3)$ 아래의 복잡도에 큰 의미는 두지 않기로 합니다.

사실 Gaussian elimination은 upper triangular matrix를 만들기 위한 하나의 수단에 불과합니다. Upper triangular matrix의 determinant는 대각선 entry의 합이고, Gaussian elimination의 과정에서 determinant는 (up to sign) 변하지 않기 때문입니다. 물론 언제 부호가 바뀌는지도 추적할 수 있습니다.

# n: size of matrix a
# pseudo-code not guaranteed to be run correctly
def get_det(a):
    det = 1
    for i in range(n):
        for j in range(i, n):
            while a[j][i]:
                swap(a[j], a[i])
                det *= -1
                a[j] -= (a[j][i] / a[i][i]) * a[i]  
        # now a[0..i, 0..i] is upper triangular  
        det *= a[i][i]
    return det

pseudo-code를 보시면 gaussian elimination 과정에서 나눗셈이 사용된 걸 선명하게 확인하실 수 있습니다.

Characteristic polynomial

이 문단의 대부분은 Cohen과 https://rkm0959.tistory.com/141 를 참고하여 작성하였습니다.

Characteristic polynomial $\phi _ {A}(x)$는 일반적으로 eigenvalue problem에 많이 등장하지만, 각각의 계수들도 조합적인 의미를 가집니다. 한 예로, $i$차항의 계수 $[x^{i}]\phi _ {A}(x)$는 $\lbrace1, \cdots, n\rbrace$의 크기 $n-i$짜리 부분집합 $I$에 대해, $I$에 속한 행, 열만 추려낸 submatrix $A _ {I, I}$의 determinant를 적절히 더하고 뺀 것과 같습니다. 실제로 $n-1$차항의 계수 $-\mathrm{tr}(A)$, 상수항 계수 $(-1)^{n} \det A$가 이에 들어맞는 것을 확인하실 수 있습니다.

$xI - A \in R[x]^{n \times n}$에도 Gaussian elimination이 되면 참 좋겠지만, 당장 $xI - A$의 diagonal element인 $x - a \in R[x]$의 곱셈에 대한 역원도 없는 판에 나눗셈을 시도할 수는 없습니다.

하지만 Field의 또다른 좋은 성질 중 하나는, 다항식의 interpolation이 가능하다는 점입니다.

Lemma. (Lagrange interpolation) Field $R$의 서로 다른 $n$개의 원소 $\lambda _ {1}, \cdots, \lambda _ {n} \in R$과 $\mu _ {1}, \cdots, \mu _ {n} \in R$에 대해, $p(\lambda _ {i}) = \mu _ {i}$를 만족하는 유일한 $n-1$차 이하의 다항식 $p \in R[x]$를 $O(n \log^{2} n)$ 시간에 찾을 수 있다.

$\lambda \in R$에 대해 $\lambda I - A \in R^{n \times n}$이고, $\phi _ {A}(\lambda) = \det(\lambda I - A)$이니 서로 다른 $n$개의 $\lambda$에 대해 $\det$을 계산하면 $O(n^{4})$에 Characteristic polynomial을 복원할 수 있습니다.

서로 다른 $n$개의 $\lambda _ {1}, \cdots, \lambda _ {n}$에 대해 $\det(\lambda _ {k}I - A)$를 계산하는 건 좀 더 빠르게 할 수 있는데, $A$를 좋은 행렬로 similar transform하는 것입니다.

어떤 matrix $P$가 있어서 $A = PBP^{-1}$이라면 $A$는 $B$와 similar하다고 합니다. 두 matrix가 similar하면 두 matrix의 determinant와 char.poly는 같게 됩니다.

가령 upper triangular matrix $D \in R^{n \times n}$에 대해 $A = PDP^{-1}$로 나타낼 수 있다면, $\det(\lambda I - D)$는 $O(n)$ 시간에 계산해버릴 수 있습니다.

사실 이 경우는 interpolation을 적용할 필요도 없는 것이 $D$의 diagonal entry를 $d _ {1}, \cdots, d _ {n}$이라고 하면, $\phi _ {A}(x) = (x - d _ {1}) \cdots (x - d _ {n})$으로 바로 나타내어지기 때문입니다. 하지만 $A = \begin{pmatrix} 0 & 2 \ 1 & 0 \end{pmatrix} \in \mathbb{Q}^{n \times n}$과 같은 경우를 보면 $\phi _ {A}(x) = x^{2} - 2 = (x - d _ {1})(x - d _ {2})$를 만족하는 $d _ {1}, d _ {2} \in \mathbb{Q}$는 존재하지 않습니다. 즉, 이러한 접근은 항상 성공하진 못합니다.

놀랍게도, 모든 matrix는 similar transform만을 이용해 Hessenberg form을 만들 수 있습니다. Hessenberg form이란 아래와 같이 diagonal 아래 한 칸까지 nonzero element가 존재할 수 있는 matrix를 말합니다.

Hessenberg matrix $H \in R^{n \times n}$에 대해, $\det H$는 Gaussian elimination으로 $O(n^2)$ 시간에 계산할 수 있습니다. Gaussian elimination에서 시간복잡도는 대략 $n \times (\text{Row operation의 횟수})$정도로 결정되는데, row operation의 횟수가 $O(n)$회에 불과하기 때문입니다.

따라서 $A = PHP^{-1}$이라고 할 때, $\det(\lambda I - A) = \det(\lambda I - H)$ 역시 $O(n^2)$에 계산할 수 있고, 전체 시간복잡도는 $O(n^3)$입니다.

이제 임의의 matrix를 $O(n^3)$ 시간에 Hessenberg form으로 만들 수만 있으면 됩니다. 이 역시 Gaussian elimination과 같은 Elementary row/column operation으로 구현할 수 있습니다.

두 행을 swap하는 row operation $S _ {ij}$에 대해, similarity transform $S _ {ij} A S _ {ij}^{-1}$은 $A$의 $i$행과 $j$행을 swap한 뒤, $i$열과 $j$열 또한 swap하는 연산으로 나타납니다.

행 $j$에 행 $i$의 $x$배를 더하는 row operation $T _ {ij}^{x}$에 대해, similarity transform $T _ {ij}^{x} A (T _ {ij}^{x})^{-1}$ 은 $A$의 $j$행에 $i$행의 $x$배를 더한 뒤, $i$열에 $j$열의 $-x$배를 더하는 연산으로 나타납니다.

이 두 가지 similarity transform만을 이용하여 hessenberg form을 만들 수 있고, row/column operation은 $O(n^2)$ 회 일어나므로 시간복잡도는 $O(n^3)$입니다.

# n: size of matrix a
# pseudo-code not guaranteed to be run correctly
def make_hessenberg(a):
    for i in range(n-1):
        for j in range(i+1, n):
            while a[j][i]:
                swap(a[j], a[i+1])
                swap(a[:, j], a[:, i+1]) # to guarantee similarity transform

                quotient = (a[j][i] / a[i+1][i])
                a[j] -= quotient * a[i+1]
                a[:, i+1] += quotient * a[:, j]

        # now a[0..i+1, 0..i+1] is upper hessenberg
    return a

Similarity transform이 기존에 만들어둔 hessenberg matrix를 깨뜨리지 않는다는 것 또한 주의깊게 들여다보면 알 수 있습니다.

연습문제: BOJ 19562 Matrix and Queries

Pitfall: What if field is too small?

웬만한 경우에 Lagrange interpolation은 잘 동작하지만, Field가 너무 작아서 ($\mathbb{Z} _ {2}$ 등) distinct한 point를 $n$개씩이나 잡지 못할 수도 있습니다.

cf) $\mathbb{Z} _ {p}$에서 가장 먼저 생각할 수 있는 방법은, $p^{d} > n$인 $d$를 잡은 뒤, $\mathrm{GF}(p^{d})$로 체를 확장하는 것입니다. $\mathrm{GF}(p^{d})$에서도 lagrange interpolation이 유일하게 결정되므로, $\phi _ {A}(x) \in \mathbb{Z} _ {p}[x] \subseteq GF(p^{d})[x]$를 구할 수 있습니다.

하지만 Determinant의 성질을 이용하면 hessenberg matrix에서 characteristic polynomial을 $O(n^3)$ 만에 계산할 수 있습니다.

Hessenberg matrix $H = (h _ {ij})$에 대해, $p _ {k}(x)$를 $xI - H$의 row $1, \cdots, k$, column $1, \cdots, k$로 만든 $k \times k$ 행렬의 determinant라고 하면 아래의 점화식이 성립하기 때문입니다.

$p _ {n}(x)$가 $\det(xI - H)$가 되고, 이는 $O(n^3)$ 만에 계산할 수 있습니다.

번외: Characteristic Polynomial in Subcubic time?

이 글의 맥락과는 무관하지만, 다음 글의 주제를 조사하다가 알게 된 사실입니다.

Theorem (Labahn16). Finite field $\mathbb{F}$와 Non-singular matrix $B \in \mathbb{F}[x]$에 대해, $\det B$는 $\tilde{O}(n^{\omega} \lceil \mu \rceil)$ 시간에 계산할 수 있다. 이 때 $\mu$는 각 column별 max degree를 모든 column에 대해 평균한 값, 각 row별 max degree를 모든 row에 대해 평균한 값보다 작거나 같다.

이를 characteristic polynomial의 경우에 대입하면 $\lceil \mu \rceil = 1$이 되므로, $\tilde{O}(n^{\omega})$에 계산할 수 있다는 뜻이 됩니다. 보다 자세한 내용은 다음 글에서 다루도록 하겠습니다.

When $R$ is Euclidean

이 경우는 “나눗셈이 있다” 보다는 좀 더 추상적이지만, 우리에게 너무나 친숙한 예시 $\mathbb{Z}$가 있습니다.

두 정수 $a, b$에 대해, 항상 $a = bq + r$ 꼴로 쓸 수 있습니다. 여기서 중요한 것은, $0 \le r < b$가 되도록 쓸 수 있습니다.

정수 집합 $\mathbb{Z}$는 Field가 아니기에 $i$행에 $j$행의 $-a _ {ik} / a _ {jk}$배를 더해서 $a _ {ik}$를 완전히 소거해버릴 수는 없으나, $a _ {ik}$의 크기를 줄일 수 있다는 말이 됩니다.

특히나 정수에서는 연산 2번으로 $a _ {ik}, a _ {jk}$를 절반으로 줄일 수 있으므로, $O(n^{2} \log X)$ 번의 row operation으로 gaussian elimination과 동등하게 upper triangular matrix를 만들 수 있습니다.

이 때 $\mathbb{Z}$에서는 $X$가 bounded되어 있지 않음에 주의하세요. 실제로는 동일한 방법이 적용되는 $\mathbb{Z} _ {m}$에서 $O(n^{2} \log m)$ 번의 row operation을 사용하여 $\det A \mod m$을 구할 수 있게 됩니다. 실제로는 row operation의 횟수가 $O(n^{2} + n \log m)$ 정도가 되는 것도 보일 수 있지만, 정확한 증명은 생략합니다.

그래서 의사코드도 바뀌는 것이 없습니다.

# n: size of matrix a
# pseudo-code not guaranteed to be run correctly
def get_det(a):
    det = 1
    for i in range(n):
        for j in range(i, n):
            while a[j][i]:
                swap(a[j], a[i])
                det *= -1
                # x // y := an element q minimizing |x - qy|
                a[j] -= (a[j][i] // a[i][i]) * a[i]  
        # now a[0..i, 0..i] is upper triangular  
        det *= a[i][i]
    return det

Hessenberg form 또한 동일하게 만들 수 있고, 앞서 설명한 recursion 때문에 char.poly도 $O(n^3)$ 시간에 구할 수 있습니다.

Euclidean Domain

$\mathbb{Z}$와 같이, 몫-나머지 표현으로 나머지의 “크기”를 더 줄어들게 만들 수 있는 종류의 domain을 Euclidean domain이라고 부릅니다.

추상적으로는 함수 $\delta : R \to \mathbb{Z} _ {\ge 0}$을 잘 잡으면, 모든 $a, b$에 대해 $a = bq + r$이고 $r = 0$ 혹은 $\delta(r) < \delta(b)$를 만족하는 어떤 $q, r$이 있으면 됩니다. 이 $\delta$를 Valuation, 혹은 Euclidean Norm이라고 합니다.

$\mathbb{Z}$(나 $\mathbb{Z} _ {m}$)의 경우, $\delta(x) = \lvert x \rvert$가 됩니다. Field $F$에 대해 $F[x]$는 Euclidean domain이 되는데, “다항식의 나눗셈”이 되기 때문에 $f \in F[x]$에 대해 $\delta(f) = \deg f$로 두면 됩니다. 물론 이런 경우에 $\delta(f)$가 Gaussian elimination에서 $\mathbb{Z}$만큼 빠르게 줄어들지 않기 때문에, 경우에 따라 $O(n^3)$보다 느린 시간 복잡도에 돌게 될 수 있습니다.

특이한 것은, $\mathbb{C}$의 “격자점”들을 모아놓은 Gaussian integer $\mathbb{Z}[i] = \lbracea + bi : a, b \in \mathbb{Z}\rbrace$ 또한 Euclidean domain이 된다는 것입니다. 따라서 우리는 Gaussian integer matrix의 determinant 또한 빠르게 구할 수 있습니다. 이 영상에서 이에 대한 기하적인 설명을 볼 수 있습니다.

Principal Ideal Domain

cf) 이 문단은 PID에 대한 지식을 가정하는, 글의 전체 맥락과 무관한 문단입니다.

Theorem. (Dedekind-Hasse) Integral domain $R$에 대해, $R$이 PID인 것과 다음 조건을 만족하는 Dedekind-Hasse norm $N : R \to \mathbb{Z} _ {\ge 0}$이 존재하는 것은 동치이다.

• 모든 $a, b$에 대해 $a$가 $b$의 배수이거나, $pa = qb + r$, $N(r) < N(b)$인 $p,q,r \in R$이 존재한다.

우리가 알고 있는 대부분의 PID는 Euclidean Domain이고, 그렇지 않은 예시는 $\mathbb{Z}[\frac{1 + \sqrt{-19}}{2}]$와 같은 극히 특수한 예시뿐입니다.

하지만 PID $R$에 대해

• Dedekind-Hasse norm $N$을 알고 있고
• $a, b$에 대해 $N(sa - tb) < N(b)$인 $s, t$를 찾을 수 있으며
• $N(sa - tb)$가 충분히 빠르게 decay한다면

$a \leftarrow sa - tb$ 또한 elementary row operation으로 표현되는 영역이므로, Gaussian elimination을 이용하여 Euclidean domain과 같이 determinant를 구할 수 있습니다. 이 때 한 row를 $s$배 scale하는 연산이 있으므로 당연히 determinant도 $s$배 scale해줘야 합니다.

Hessenberg form에 대해서는 이야기가 좀 더 어려워집니다. 한 row를 $s$배하는 연산을 similarity transform으로 바꿔야 하는데, $s$의 역원이 정의되지 않으면 similarity transform을 정의하기 곤란해지기 때문입니다. 따라서 이러한 특이 케이스에는 determinant는 쉽고, char poly의 계산은 까다로울 것이라 예상할 수 있습니다.

이 경우에 지문을 따로 할애하는 이유는, Gaussian Elimination으로 determinant, char poly를 구하는 기법에 상당히 강력한 한계를 제공하기 때문입니다. Ring이 PID가 아닌 경우 어떤 elementary row operation을 사용하더라도, off-diagonal element의 norm을 $0$으로 만들지 못하는 경우가 발생할 수 있습니다. 이러한 결점이, 드디어 등장하는 오늘의 메인 주제인 division-free determinant algorithm의 필요성을 대두시키는 요인 중 하나입니다.

$R$이 Generic Ring일 때

이 단원에서는 $R$에 어떠한 부가적 조건도 없는 경우에 determinant, characteristic polynomial을 $O(n^4)$ 번의 ring operation (덧셈, 곱셈)에 계산하는 알고리즘을 소개할 것입니다. 해결법은 놀랍게도 “simple”하다고 칭할 법한 DP입니다.

Determinant의 계산식의 각 항을 Sign $\mathrm{sgn} \sigma$와, 가중치 $a _ {1\sigma _ {1}} \cdots a _ {n \sigma _ {n}}$으로 분리해서 봅시다.

이 때, 모든 $i$에 대해 정점 $i$에서 정점 $\sigma _ {i}$로 가중치 $a _ {i \sigma _ {i}}$의 간선을 잇는다고 생각하면 각 순열은 정점 $\lbrace1, \cdots, n\rbrace$의 Disjoint (directed) cycle cover와 일대일대응됩니다. Cycle cover에 사용된 disjoint cycle들을 $C _ {1}, \cdots, C _ {k}$라고 둡시다.

Cycle cover $\sigma$의 weight $w(\sigma)$를 $\prod _ {i} a _ {i \sigma _ {i}}$라고 정의하고, $\mathrm{sgn}(\sigma) = (-1)^{n-k}$가 되는 것을 이용해 아래와 같이 쓸 수 있습니다.

이제 Cycle cover를 더 generalize한 CloW sequence (ClosedWalk Sequence)를 정의합시다.

CloW란 $1$개 이상의 정점으로 구성된 닫힌 경로 (정점 중복을 허용) $C$ 중 다음 조건을 만족하는 것을 말합니다.

• $C$에 속한 정점 중 번호가 가장 낮은 정점 $\mathrm{head}(C)$는 outdegree, indegree 모두 $1$이다. (단 한 번만 방문한다)

닫힌 경로 $(1), (1, 2, 1), (1, 3, 3, 1), (1, 4, 6, 4, 1)$ 등은 CloW이지만, $(1, 1, 1), (1, 2, 1, 2, 1)$ 등은 CloW가 아닙니다.

$\lbrace1, \cdots, n\rbrace$의 CloW sequence $\mathcal{C} = (C _ {1}, \cdots, C _ {k})$란 $\mathrm{head}(C _ {1}) < \cdots < \mathrm{head}(C _ {k})$를 만족하고, 사용된 간선이 정확히 $n$개인 clow $C _ {1}, \cdots, C _ {k}$를 말합니다. Cycle cover는 당연히 CloW sequence가 되고, Cycle cover의 것을 일반화하여 $\mathrm{sgn}(\mathcal{C}) = (-1)^{n-k}$, $w(\mathcal{C}) = \prod _ {e \in E(\mathcal{C})} w _ {e}$로 정의합니다. (regarding multiplicity)

Theorem (Mahajan97).

Proof. CloW sequence 사이에 involution $\phi$를 잡는 방법을 이용합니다. $\phi^{2} = \mathrm{id}$이되, $\phi$는 Cycle cover는 보존하고 cycle cover가 아닌 clow sequence는 sign term만 반대인 다른 CloW sequence로 보내는 일대일대응입니다. $\phi$의 existence만 보이면 non-cycle cover끼리 서로 cancel out하는 것이 당연합니다.

$\mathcal{C} = (C _ {1}, \cdots, C _ {k})$를 고정합시다. $C _ {i+1}, \cdots, C _ {k}$가 disjoint cycle이 되는 $i$의 최솟값을 잡습니다. $i = 0$이라면 $\mathcal{C}$는 cycle cover이므로 $\phi(\mathcal{C}) = \mathcal{C}$. $i > 0$인 경우 $C _ {i}$를 $\mathrm{head}(C _ {i})$부터 traverse하다보면 다음 두 경우 중 하나를 무조건 만나게 됩니다.

• $C _ {i+1}, \cdots, C _ {k}$에 속하는 정점 중 하나를 만난다.

  • 이 경우 만나는 cycle을 $C _ {j}$라고 했을 때 $C _ {j}$와 $C _ {i}$를 merge해서 만든 CloW $D$를 넣고 $C _ {i}, C _ {j}$를 제거한 것을 $\phi(\mathrm{C})$로 둡니다. 간선 개수는 보존되고, $\mathrm{head}(D) = \mathrm{head}(C _ {i})$입니다. CloW 개수가 하나 줄어드니 sign만 바뀝니다.

• $C _ {i}$의 정점 $v$를 두 번째로 방문하면서 $v$의 첫방문 $\sim$ 두번째 방문 구간이 simple cycle을 이룬다.

  • 이 경우 $v$로 인해 생긴 사이클을 잘라내서 CloW sequence에 새로 집어넣습니다.

첫 번째 경우의 construction이 두 번째 경우로 가고, 두 번째 경우의 construction이 첫 번째 경우로 가므로 non-cycle-cover CloW sequence 중 $\phi$의 fixed point는 없습니다. 또 $\phi$가 involution이 되는 것 또한 어렵지 않게 알 수 있습니다.

따라서 CloW sequence들을 DP로 잘 세어 주면 됩니다.

$D _ {i}(h, u)$를 $i$번째 edge까지 써서 만든 CloW의 head가 $h$이고, $u$까지 traverse한 경우의 weight 합이라고 합시다.

이 때 $i$번째 edge를 사용할 수 있는 경우는 다음 두 가지 뿐입니다.

• $u$에서 다른 정점 $v > h$를 탐색해서 CloW를 늘리는 경우, $D _ {i}(h, v) += a _ {uv}D _ {i-1}(h, u)$
• $u$에서 head $h$로 돌아와 CloW를 닫는 경우, CloW를 닫고 새 head를 찾아주어야 합니다.

이 때 $D _ {n-1}$에서 Closed CloW sequence들에 대한 값을 잘 더하면 determinant를 counting하게 됩니다.

Characteristic polynomial의 경우, 매우 비슷한 DP 풀이가 됩니다. 정확히 우리가 카운팅해야 하는 것은 loop의 개수가 $k$개인 cycle cover의 개수인데, 이는 기존의 관계식 및 DP를 이용하여 구할 수 있으니 시도해보시기 바랍니다.

References

• Labahn16: Labahn, George and Neiger, Vincent and Zhou, Wei, Fast, deterministic computation of the Hermite normal form and determinant of a polynomial matrix.
• rkm0959 - Fast Computation of Characteristic Polynomial
• Henri Cohen, A Course in Computational Algebraic Number Theory
• Arvind09: On the Hardness of the Noncommutative Determinant
• Mahajan97: Determinant: Combinatorics, Algorithms, and Complexity (PDF 다운로드 주의)