서론
きたまさ法은 일본의 경기 프로그래머 wata가 블로그에서 소개한 점화식을 빠르게 계산하는 방법이다. 사실상 companion matrix의 멱승을 구하는 방법과 같지만, 선형 점화식이라는 특성을 활용하여 FFT를 이용, 선형 점화식의 항을 약간 더 빠르게 계산할 수 있다. ( \( O(k^2 \log N) \)에서 \( O(k \log k \log N )\). 따라서 \(k\)가 작은 점화식의 경우 별 의미가 없다. ) 이 테크닉의 이름은 블로그에서 소개한 wata가 붙였는데 きたまさ가 이 방법을 강력하게 주장했다고 한다. 이 きたまさ라는 사람이 한국에서는 잘 알려져 있지 않지만, kitamasa란 핸들을 이용하는 北川宜稔라는 사람으로, 일본에서는 '개미책'이라고 불리는 「プログラミングコンテストチャレンジブック」라는 책을 저술했다. (wata씨는 공동 저자 중 한 명이다.)
Companion matrix를 이용한 선형 점화식 계산
일단, 다항식에 대한 companion matrix에 관한 내용은 wikipedia나 linear algebra 관련 책을 찾아보길 바란다. 이 챕터에서는 companion matrix에 대해 알고 있다고 생각하고, 이를 통해 \(m\)개의 항으로 표현되는 선형 점화식의 \(n\)번째 항을 \(\Theta(m^3 \log(n) ) \)만에 구하는 방법을 설명한다.
먼저, \(m\)개의 항으로 이루어진 선형 점화식을 다항식으로 간주하여 companion matrix를 구해야 한다. 양의 정수 \(m\), \(n\)에 대해 \(a_{n+m}=c_0 a_n + c_1 a_{n+1} + ... + c_{m-2} a_{n+m-2}+ c_{m-1} a_{n+m-1} \)로 표현되는 선형 점화식이 있다고 하자. 이 점화식을 다항식 \( p(x)= c_0 x^0 +c_1 x^1 + ... +c_{m-2} x^{m-2} + c_{m-1} x^{m-1} \)으로 간주하여 이 \(p(x)\)의 companion matrix \(C(p(x))\)를 구할 수 있다. 이때, 행 벡터 \(A_n=(a_n,a_{n+1},...,a_{n+m-2},a_{n+m-1})\)이라 하자. 그러면, \(A_{n+1}^T = C(p(x))A_{n}^T\)이다.
따라서, \(A_{n+1}^T = C(p(x))C(p(x))A_{n-1}^T = ... = C(p(x))^n A_1 ^T\) 을 만족한다.
행렬의 거듭제곱을 구할 때 exponentiation by squaring을 이용하여 \(\Theta(m^3 \log(n) )\)만에 할 수 있으므로 총 시간복잡도는 점근적으로 \((\Theta(m^3 \log(n)) + \Theta(m)) = \Theta( m^3 \log(n) )\) 이다.
きたまさ法 Kitamasa method
\(n \times n\) 행렬 \(M\)이 있다고 하자. 이 행렬의 제곱을 구하는 데에는 \(\Theta(n^3)\)만큼의 연산이 필요하다. 물론 약간 더 빠른 방법도 있지만, 현재까지 알려진 방법의 시간복잡도는 모두 \(\omega(n^2)\)이며, 점근적 표기법만으로 보면 꽤 빠른 것처럼 보이는 알고리즘 또한, 점근적 표기에 의해 삭제된 상수가 크기 때문에 일반적인 경우에서는 naive 한 방법보다 느리다.
우리가 목표로 하는 것을 생각해 보자. 선형 점화식으로 정의된 수열의 \(n\)번째 항을 빠르게 구하는 것이다. 보통 선형 점화식을 이루는 항의 수는 크지 않아 companion matrix를 이용하는 것이 충분히 효율적이다. 그렇지 않은 경우도 확실히 존재한다. きたまさ法은 이런 부분을 계산하기 위한 기법이다.
이전 챕터에서 사용한 점화식의 형태를 보자.
$$ a_{n+m}=c_0 a_n + c_1 a_{n+1} + ... + c_{m-2} a_{n+m-2}+ c_{m-1} a_{n+m-1} $$
이 점화식의 \(n\)번째 항을 아래와 같이 표현하자.
$$ a_n = \sum_{i=0}^{m-1} c_{n,i}a_i $$
\(c_{m,i}\)를 통해 \(c_{n,i}\)를 빠르게 계산할 수 있다면, \(a_n\) 도 빠르게 계산할 수 있을 것이다.
\(0\) 이상의 정수 \(k\)에 대해 \(a_{n+k}\)를 아래와 같이 표현할 수 있다.
$$ a_{n+k}=\sum_{i=0}^{m-1} c_{n,i}a_{i+k} $$
이 부분이 잘 이해가 안 된다면, companion matrix를 이용한 방법을 생각해 보면 된다. \(C(p(x))^{n} \times A_{k}^T = A_{k+n}^T\)를 이용하여 \(a_{n+k}\)를 계산하는 것과 같다고 생각해도 좋다.
이제 \(c_{n,i}\)를 \(c_{m,i}\)로 표현하기 위해 \(a_{n+1}\)의 경우를 고려하겠다.
$$ \begin{eqnarray} a_{n+1} &=& \sum_{i=0}^{m-1}c_{n+1,i} a_{i} \\ &=& \sum_{i=0}^{m-1}c_{n,i}a_{i+1} \\ &=& \sum_{i=0}^{m-2}c_{n,i}a_{i+1} + c_{n,m-1}a_m \\ &=& \sum_{i=1}^{m-1}c_{n,i-1}a_{i} + c_{n,m-1} \sum_{i=0}^{m-1} c_{m,i}a_i \\ &=& \sum_{i=1}^{m-1} \left[ \left( c_{n,i-1} + c_{n,m-1} c_{m,i} \right) a_i \right] + c_{n,m-1}c_{m,0}a_{0} \end{eqnarray} $$
따라서,
$$ c_{n+1,i} = \begin{cases} c_{n,i-1} + c_{n,m-1}c_{m,i} & \mbox{if } i>0 \\ c_{n,m-1}c_{m,0} & \mbox{if } i=0 \end{cases} $$
\(c_{n+1,i}\)를 \(c_{n,i}\)와 \(c_{m,i}\)를 이용하여 표현하였다. 이를 이용하여 \(c_{n,i}\)를 \(c_{m,i}\)를 통해 표현할 수는 있다만 이래서는 naive 하게 점화식을 계산하는 것보다 느릴 수 있다. 이제 약간 더 효율적인 방법을 찾기 위해 앞서 확인한 사실을 통해 \(a_{2n}\)을 \(a_{n}\)으로 표현해 보겠다.
$$ \begin{eqnarray} a_{2n} &=& \sum_{i=0}^{m-1}c_{2n,i}a_i \\ &=& \sum_{i=0}^{m-1} \left[ c_{n,i} \left( \sum_{j=0}^{m-1}c_{n+i,j}a_{j} \right) \right] \\ &=& \sum_{i=0}^{m-1}\left[ \left(\sum_{j=0}^{m-1} c_{n,j}c_{n+j,i}\right)a_{i} \right] \end{eqnarray} $$
세 번째 줄의 식으로부터,
$$ c_{2n,i} = \sum_{j=0}^{m-1}c_{n,j}c_{n+j,i} $$
앞서 구한 두 가지 식을 이용하면 \(c_{n,i}\)를 \(c_{m,i}\)로 \(O(m^2 \log{n})\)번의 연산을 이용해서 나타낼 수 있다. 그러면, \(a_n\)을 \(\Theta(m)\)번의 연산을 이용하여 계산할 수 있다.
Source
https://github.com/Sora-Sugiyama/Libs/blob/main/Linear-Algebra/kitamasa.h
'Algorithms' 카테고리의 다른 글
Bowyer-Watson Algorithm : Online Delaunay Triangulation (1) | 2023.10.31 |
---|---|
Knuth-Morris-Pratt Algorithm (1) | 2023.10.02 |
Dijkstra's Algorithm (1) | 2023.09.11 |
Segment Tree with Bitset (0) | 2023.08.09 |
이분 탐색 (0) | 2023.05.08 |