백준 17315번: Matrix Game (matrix interpretation)

 

https://www.acmicpc.net/problem/17315

17315번: Matrix Game

---

National Olympiad in Informatics, China, 2013, Day2에 Problem1로 출제된 문제라고 합니다.

 

이 문제에 대해 행렬을 사용한 점화식, 일반적인 정수 계수 점화식 두 방법을 이용하여 풀었는데, 복잡도는 비슷하지만 상수 차이로 실행시간 차이가 컸습니다. 그리고 두 풀이가 꽤 상이한 듯하여 둘 다 소개하고자 합니다.

 

정수 거듭제곱 풀이에 대해 궁금하신 분은 이쪽 링크를 이용하시면 됩니다: Matrix Game (integer interpretation)

---

행렬 풀이는 몇 번만 점화식을 이어가보면 금방 행렬 점화식을 찾을 수 있다는 게 장점이며, 아쉬운 건 행렬이 갖는 몇 가지 제약이 있어 코딩의 난이도가 올라간다는 것입니다. 개인적으로는 수식으로는 행렬이 더 깔끔하지 않나라고 생각합니다.

 

우선 주어진 규칙에서 $\left ( j \ne 1 \right )$일 때를 이용하여, 충분히 큰 $n$과 $m$에 대하여 일반적인 $F[n][m]$에 관한 점화식을 유도합니다.

$$ F[n][m] = a \times F[n][m-1] + b$$

이므로 이런 주제에 익숙하다면, 다음을 쉽게 관측할 수 있습니다.

$$ \begin{bmatrix} F[n][m] \\ 1 \end{bmatrix} = \begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix} \begin{bmatrix} F[n][m-1] \\ 1 \end{bmatrix} $$

따라서, $m>1$에서 다음이 성립합니다.

$$ \begin{bmatrix} F[n][m] \\ 1 \end{bmatrix} = \begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix}^{m-1} \begin{bmatrix} F[n][1] \\ 1 \end{bmatrix} $$

주어진 규칙에서 $\left ( i \ne 1 \right )$일 때를 이용하여 다음을 얻습니다. (저희는 충분히 큰 $n$, $m$에 대해 생각하고 있습니다!)

$$ F[n][1] = c \times F[n-1][m] + d$$

그러므로 다음 식이 나옵니다.

$$ \begin{bmatrix} F[n][1] \\ 1 \end{bmatrix} = \begin{bmatrix} c & d \\ 0 & 1 \end{bmatrix} \begin{bmatrix} F[n-1][m] \\ 1 \end{bmatrix} $$

이를 위에서 구했던 식에 대입합니다.

$$ \begin{bmatrix} F[n][m] \\ 1 \end{bmatrix} = \begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix}^{m-1} \begin{bmatrix} c & d \\ 0 & 1 \end{bmatrix} \begin{bmatrix} F[n-1][m] \\ 1 \end{bmatrix} $$

이제 $n>1 \land m>1$에서 다음이 성립함을 알 수 있습니다!

$$ \begin{bmatrix} F[n][m] \\ 1 \end{bmatrix} = \left \{ \begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix}^{m-1} \begin{bmatrix} c & d \\ 0 & 1 \end{bmatrix} \right \} ^{n-1} \begin{bmatrix} F[1][m] \\ 1 \end{bmatrix} $$

마지막으로 $\left ( i \ne 1 \right )$ 일 때의 규칙을 이용하여 일반항을 얻습니다.

$$ \begin{bmatrix} F[n][m] \\ 1 \end{bmatrix} = \left \{ \begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix}^{m-1} \begin{bmatrix} c & d \\ 0 & 1 \end{bmatrix} \right \} ^{n-1} \begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix}^{m-1} \begin{bmatrix} F[1][1] \\ 1 \end{bmatrix} $$

마무리로 정의에서 $F[1][1]=1$을 가져오면 저희가 계산할 것은 다음과 같습니다. (계산결과의 $[0][0] + [0][1]$ 이 답!)

$$\left \{ \begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix}^{m-1} \begin{bmatrix} c & d \\ 0 & 1 \end{bmatrix} \right \} ^{n-1} \begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix}^{m-1} $$

물론 위에서 논의했던 대로, 이것은 $n>1 \land m>1$에서 성립합니다.

---

이제 $n=1 \lor m=1$일 경우를 생각해보겠습니다.

 

$n=1$이고 $m>1$인 경우, 계산해야할 식은 다음과 같습니다.

$$\begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix}^{m-1} $$

$n>1$이고 $m=1$인 경우, 계산해야할 식은 다음과 같습니다.

$$ \begin{bmatrix} c & d \\ 0 & 1 \end{bmatrix} ^{n-1} $$

 

그런데 이것은 결국 $n>1 \land m>1$ 식에서 $n=1$ 또는 $m=1$을 넣었다고 가정했을 때와 같습니다.

물론 $n=m=1$인 경우에도 행렬의 0제곱이 identity matrix임을 이용하면 정의와 부합하는 답을 얻습니다.

 

따라서 코드의 분량을 줄이고 싶다면 분기를 파지 않고 그냥 저 위의 일반항을 그대로 구현해도 됩니다.

---

이제 이것을 코드로 구현하려고 할 때, 인풋의 범위를 확인하고 잠시 당황스러울 수 있습니다. $n,m \leq 10^{1000000}$으로 주어진 것은 무엇을 의도한 걸까요?

 

수학적으로는 아무 문제 없지만, C++ 또는 C를 쓰는 프로그래머 입장에서 행렬의 $10^{10^{6}}$ 승은 매우 곤란합니다. 이유는 $10^{10^{6}}$ 을 저장할 방법이 문자열 또는 복잡한 BigInteger 라이브러리 뿐이기 때문입니다.

 

백준과는 별도로, Boost 같은 상용 BigInteger 라이브러리를 가져온다고 해도 한계가 있습니다. 애초에 $\log _{2} 10^{10^{6}}$ 값을 생각해보면 $10^{10^{6}}$을 이진수로 저장하기 위해 3백만 비트 이상이 필요하기 때문입니다. 이정도 규모의 연산은 어떤 라이브러리를 가져와도 시간제한 1초에 들어갈 수 없을 것입니다. 다른 언어라면 상황은 더 악화될 것입니다.

 

지수의 크기를 유의미하게 줄이려고 해도 행렬의 거듭제곱이므로, 정수론에서의 Euler's Theorem 처럼 획기적으로 지수의 크기를 줄일 방법도 잘 알려진 게 없습니다.

 

따라서 언어를 막론하고, $n$과 $m$을 문자열로 저장했다고 가정한 뒤 새로운 알고리즘을 고안해야 합니다.

 

참고로 PS에서 쓰이곤 하는 BigInteger Template들을 활용하여도 기존의 exponentiation by squaring 전략이 가능은 하겠으나, 시간제한이 문제인지 AC 코드중에서는 그런 방법이 나타나진 않았습니다.

 

Exponentiation by squaring 알고리즘의 원리를 생각해보면, 지수를 binary로 변환했을 때 각 자릿수에는 1과 0이 있습니다.

 

그러므로 Exponentiation by squaring 알고리즘을 러프하게 다음처럼 생각할 수 있습니다.

$$ a^{p} = a^{ b_{l}b_{l-1} \cdots b_{2}b_{1} } = \prod_{i=1}^{l} a^{b_{i} \times 2^{i-1}} = \prod_{i=1}^{l} \left ( a^{2^{i-1}} \right )^{b_{i}} $$

여기서 $a^{2^{i-1}}$ 끼리는 $i$가 1씩 커짐에 따라 base case를 제곱해나가면 구할 수 있다는 게 핵심 아이디어였습니다.

 

그렇다면 $n$을 십진수 문자열로 저장했을 땐 어떻게 변할까요?

$$ a^{n} = a^{ d_{l}d_{l-1} \cdots d_{2}d_{1} } = \prod_{i=1}^{l} a^{d_{i} \times 10^{i-1}} = \prod_{i=1}^{l} \left ( a^{10^{i-1}} \right )^{d_{i}} $$

여기서는, $a^{10^{i-1}}$ 끼리는 $i$가 1씩 커짐에 따라 base case를 10제곱해나가면 구할 수 있다는 게 핵심이 됩니다.

---

그렇게 행렬을 십진수 문자열로 저장된 지수만큼 제곱하는 알고리즘을 구현한 순간, 또다른 문제에 직면합니다. 저희가 구현해야 하는 식에는 $n-1$이나 $m-1$ 제곱이 있습니다. 즉 $n$과 $m$이 단독으로 지수에 나타나는 일이 없습니다.

 

물론 이에 대해 마지막 자릿수에서만 1을 빼는 방법을 생각할 수 있습니다만, 아무래도 받아올림이나 받아내림을 생각하면 곤란해집니다. 물론 이에 관해서는 BigInteger Templates를 써도 되겠습니다.

 

BigInteger Template 없이 $n-1$ 제곱을 구현하는 방법은? 역행렬을 사용하는 것입니다!

 

물론 역행렬이라는 것이 그렇게 단순한 대상은 아니어서, square matrix가 invertible하기 위한 필요충분조건은 그 determinant가 0이 아닐 것입니다.

 

다행히 저희가 사용하는 행렬들에 대하여

$$ \textrm{det} \begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix} = a \geq 1$$

$$ \textrm{det} \begin{bmatrix} c & d \\ 0 & 1 \end{bmatrix} = c \geq 1$$

임이 보장되므로 두 행렬은 항상 invertible합니다.

 

여기서 살짝 맘에 걸리는 것은 구해야하는 식 중에서 다음 부분인데,

$$ \left \{ \begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix}^{m-1} \begin{bmatrix} c & d \\ 0 & 1 \end{bmatrix} \right \} ^{n-1} $$

두 invertible matrices A,B에 대해 성립하는 다음 식을 생각해봅시다.

$$ \left ( AB \right )^{-1} = B^{-1} A^{-1} $$

이 성질을 통하여, 위에서 행렬 $\left \{ \begin{bmatrix} a & b \\ 0 & 1 \end{bmatrix}^{m-1} \begin{bmatrix} c & d \\ 0 & 1 \end{bmatrix} \right \}$ 가 invertible함을 보일 수 있습니다.

 

여기서 선형대수학을 활용한 PS에 익숙한 분이라면 일반적인 square matrix의 역행렬을 구해주는 코드를 미리 갖고올 수 있겠지만, 저는 그냥 2 by 2 수준의 행렬이라면 하드코딩도 나쁘지 않을 거라 판단했습니다.

 

위키피디아에 따르면 invertible한 2 by 2 행렬의 역행렬은 다음과 같습니다.

$$ A = \begin{bmatrix} a & b \\ c & d \end{bmatrix} \rightarrow A^{-1} = \frac{1}{\textrm{det} A} \begin{bmatrix} d & -b \\ -c & a \end{bmatrix} = \frac{1}{ad-bc} \begin{bmatrix} d & -b \\ -c & a \end{bmatrix} $$

---

여기서 한 가지 더 주의할 점이 있습니다. 정수 $a$, $b$에 대하여

$$ \left ( a^{m-1}b \right )^{n-1} = a^{(m-1)(n-1)}b^{n-1} $$

이지만 행렬 $A$, $B$에서 일반적으로

$$ AB \ne BA $$

이므로

$$ \left ( A^{m-1}B \right )^{n-1} \ne A^{(m-1)(n-1)}B^{n-1} $$

입니다.

 

참고로, 행렬에서 $AB=BA$ 가 성립하는 $A$, $B$의 조건에는 이런 것들이 있습니다.

https://math.stackexchange.com/questions/2699672/when-are-two-matrices-a-and-b-ab-ba

 

물론 $a$, $b$, $c$, $d$가 저런 특수한 상황을 만족하게 주어질 수 있겠지만, 굳이 분기를 나눌 필요까진 없어보입니다...

 

이제 구현만 하면 됩니다! 다음은 제 코드입니다.

#include <stdio.h>
#include <string.h>
#define mod 1000000007

long long power_modular (long long X, long long Y);
 
typedef struct {
	long long array[2][2];
} Matrix;

Matrix matrix_inverse_modular (Matrix K);
Matrix matrix_multiply_modular (Matrix A, Matrix B);
Matrix matrix_power_modular (Matrix A, int P);
Matrix matrix_power_modular_log10 (Matrix A, char* P, int L);

int main(void) {
	
	char N[1000002], M[1000002];
	long long a, b, c, d;
	scanf("%s %s %lld %lld %lld %lld", N, M, &a, &b, &c, &d);
	
	Matrix Base1 =
	{{
		{a, b},
		{0, 1}
	}};
	Matrix inv_Base1 = matrix_inverse_modular(Base1);
	
	int lengthM = strlen(M);
	Base1 = matrix_power_modular_log10 (Base1, M, lengthM);
	Base1 = matrix_multiply_modular (Base1, inv_Base1);
	
	Matrix Temp =
	{{
		{c, d},
		{0, 1}
	}};
	
	Matrix Base2 = matrix_multiply_modular (Base1, Temp);
	Matrix inv_Base2 = matrix_inverse_modular(Base2);
	
	int lengthN = strlen(N);
	Base2 = matrix_power_modular_log10 (Base2, N, lengthN);
	Base2 = matrix_multiply_modular (Base2, inv_Base2);
	
	Matrix Answer = matrix_multiply_modular (Base2, Base1);
	printf("%lld", (Answer.array[0][0] + Answer.array[0][1]) % mod);
	return 0;
}

long long power_modular (long long X, long long Y) {
	long long result = 1;
	while(Y) {
		if (Y % 2) {
			result = result * X % mod;
		}
		X = X * X % mod;
		Y /= 2;
	}
	return result;
}

Matrix matrix_inverse_modular (Matrix K) {
	for (int i = 0; i < 2; i++) {
		for (int j = 0; j < 2; j++) {
			K.array[i][j] %= mod;
		}
	}
	/* inverse of the 2x2 matrix: https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_2_%C3%97_2_matrices
	 * a b
	 * c d
	*/
	long long det = ((K.array[0][0] * K.array[1][1] - K.array[0][1] * K.array[1][0]) % mod + mod) % mod;
	long long invdet = power_modular(det, mod - 2);
	
	Matrix result =
	{{
		{K.array[1][1], mod - K.array[0][1]},
		{mod - K.array[1][0], K.array[0][0]}
	}};
	for (int i = 0; i < 2; i++) {
		for (int j = 0; j < 2; j++) {
			result.array[i][j] = result.array[i][j] * invdet % mod;
		}
	}
	
	return result;
}
 
Matrix matrix_multiply_modular (Matrix A, Matrix B) {
	Matrix result = {{{}}};
	for (int i = 0; i < 2; i++) {
		for (int j = 0; j < 2; j++) {
			for (int k = 0; k < 2; k++) {
				result.array[i][j] += A.array[i][k] * B.array[k][j];
			}
            result.array[i][j] %= mod;
		}
	}
	return result;
}
 
Matrix matrix_power_modular (Matrix A, int P) {
	Matrix result =
	{{
		{1, 0},
		{0, 1}
	}};
	while(P) {
		if (P % 2) {
			result = matrix_multiply_modular (result, A);
		}
		A = matrix_multiply_modular (A, A);
		P /= 2;
	}
	return result;
}

Matrix matrix_power_modular_log10 (Matrix A, char* P, int L) {
	Matrix result =
	{{
		{1, 0},
		{0, 1}
	}};
	for (int i = L-1; i >= 0; i--) {
		result = matrix_multiply_modular(result, matrix_power_modular(A, P[i] - '1' + 1));
		A = matrix_power_modular(A, 10);
	}
	return result;
}

(C11, 408ms, 2952KB, 제출번호 54484757)

---

처음에 풀었을 때는 행렬의 크기를 3 by 3으로 했었는데, 이 글을 쓰던 중 다른 분들의 코드를 보고 2 by 2가 가능하다는 걸 알아 수정했습니다. 3 by 3은 무려 800ms가 소요됐었습니다.

 

unsigned를 써서 380ms로 줄여보기도 했는데 제가 사용 중인 matrix 구조체 관련 함수가 전반적으로 오버헤드가 큰 것 같기도 합니다. 2 by 2 matrix multiplication을 하드코딩했을 때는 172ms까지 나오기도 했습니다. (제출번호 54493568)

 

현재 수식에서 사용하는 matrix가 2 by 2 크기이고, bottom row가 0, 1임을 이용하여 C++의 pair로 matrix multiplication을 exploit한 분들도 있었습니다. 이런건 볼 때마다 참 신기한 것 같습니다.

 

---

구글 검색을 해본 결과 Euler's Theorem을 2 by 2 matrix에 대하여 확장하려한 논문이 있었습니다. 링크

 

나름대로 해석해보니 2 by 2 matrix $A$와 modulo $m=10^{9}+7$에 대해 다음을 얻을 수 있었습니다.

$$ A^{m(m^{2}-1)} \equiv I $$

물론 determinant가 m으로 나눠떨어져야 하지 않는 조건이 있습니다만 $a,b,c,d$의 범위 상 $\textrm{det} = 0$이 아니면 신경 쓸 필요가 없고 거기에 이 문제에서 쓰는 행렬이 invertible하다는 걸 이미 보였으니 쓸 수 있는 것 같습니다. __int128 같은 자료형을 쓰면 C++로도 가능은 할 텐데, 해보지는 않아서 잘 모르겠습니다.

 

(2023-02-03 추가: hardcoded matrix multiplication + unsigned __int128 type 사용 + register keyword 사용으로 140ms까진 떨궜는데 (제출번호:55216364) 행렬로는 여기까지가 한계인지 더 줄일 수 있는 건지... 가장 큰 bottleneck이 $n,m$의 modulo $(P+1)(P-1)P$ 값을 구하는 과정 같은데, 개선이 쉽지 않아보입니다.)