Half GCD Algorithm

2023-10-05) 초안 작성

 

내용이 길어질것 같아 내용을 뺍니다.

 

1) Euclidean Algorithm과 Bezout's Identity

2) Half GCD, Full GCD
3) Half GCD, Full GCD on Polynomial

4) Half GCD Algorithm \( \mathcal{O}(n \log^2 n) \)

5) Full GCD Algorithm \( \mathcal{O}(n \log^2 n) \)

6) Inv of Polynomial(2) \( \mathcal{O}(n \log^2 n) \)

 

1) Euclidean Algorithm과 Bezout's Identity

 

두 음이 아닌 정수 \(a,b\) 가 존재한다고 가정합시다. 이 둘의 최대공약수 \( \gcd(a,b) \)를 구하고 싶습니다. 이 경우, 유클리드 알고리즘을 사용하면 \( \mathcal{O}( \log \max(a,b) ) \) 시간에 \( \gcd \)를 구할 수 있다는 것이 알려져 있습니다. 유클리드 알고리즘을 조금 더 엄밀하게 표현하면 다음과 같습니다.

 

1. \( a \geq b \)라고 가정합니다. 만약, \( a < b\)라면, \( \gcd(a,b) = \gcd(b,a) \)가 성립하기 때문에 \( \gcd(b,a) \)를 계산하면 됩니다.

2. 나머지 정리에 의하여 \(a = qb + r, r<b\)를 만족하는 음이 아닌 정수 \(q,r\)가 존재합니다.

3. \( \gcd(a,b) = \gcd(b,r) \)이 성립합니다. 이제 \( \gcd(b,r) \)에 대하여 2.를 반복해주면 됩니다. 만약, \(r=0\)이라면 \(\gcd(b,r) = b\)로 정의하고 알고리즘을 종료합니다.

 

유클리드 알고리즘은 베주의 정의를 보이는데도 쓰입니다. 베주의 정의란 다음 정리를 의미합니다:

 

두 정수 \(a,b\)의 공약수가 \(g\)라는 것은, 적당한 정수 \(s,t\)가 존재하여, \(as+bt=g\)꼴로 적을 수 있다는 것을 의미합니다. 이러한 \(g\)중 절대값이 가장 작으며, 0이 아닌 값을 최대공약수로 정의합니다.

 

\(s,t,g\)를 어떻게 계산할 수 있을까요? 유클리드 알고리즘을 행렬식으로 표현하면 쉽게 계산할 수 있습니다.

 

1. \(a = qb+r,r<b\)라고 가정합시다. 그렇다면, \( \begin{pmatrix} b \\ r \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & -q \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} \) 가 성립합니다.

2. \(b,r\)에 대하여 1.을 반복하다보면, \( \begin{pmatrix} g \\ 0 \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & -q_k \end{pmatrix} \cdots \begin{pmatrix} 0 & 1 \\ 1 & -q_1 \end{pmatrix}  \begin{pmatrix} a \\ b \end{pmatrix} \) 꼴로 표현할 수 있습니다.

3. \( M = \begin{pmatrix} m_{11} & m_{12} \\ m_{21} & m_{22} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & -q_k \end{pmatrix} \cdots \begin{pmatrix} 0 & 1 \\ 1 & -q_1 \end{pmatrix} \) 라고 정의합시다. 그러면, \(g = m_{11} a + m_{12} b\)가 성립합니다.

 

2) Half GCD, Full GCD

 

PS에서 문제를 풀때 최대공약수 \(g\)던지, \(as+bt=g\)를 만족시키는 \(s,t\)라던지를 요구하는 경우가 종종 있습니다. 뭐가 되었든 행렬 \(M\)을 계산할수만 있다면 쉽게 구할 수 있는 값들입니다. \(M\)은 유클리드 알고리즘을 계속 적용해서 구할 수 있지만, 한가지 문제가 존재합니다. 바로 행렬 곱셈을 \(k\)번 해야한다는 것입니다(여기서 \(k\)란 유클리드 알고리즘의 호출 횟수입니다). 이 문제를 해결하기 위하여 Half GCD라는 개념을 생각해 볼 수 있습니다.

 

 \( M_{HGCD}(a,b) = \begin{pmatrix} 0 & 1 \\ 1 & -q_{k/2} \end{pmatrix} \cdots \begin{pmatrix} 0 & 1 \\ 1 & -q_1 \end{pmatrix}  \) 로 정의합니다.

 

이에 대응되는, 원래 구하려고 하는 행렬 \(M\)을 Full GCD로 정의합시다.

 

\( M_{FGCD}(a,b) = \begin{pmatrix} 0 & 1 \\ 1 & -q_{k} \end{pmatrix} \cdots \begin{pmatrix} 0 & 1 \\ 1 & -q_1 \end{pmatrix}  \) 로 정의합니다.

 

이렇게 정의하면, Sparse Table이나 빠른 거듭제곱을 할 때처럼, HGCD를 \( \mathcal{O}( \log k ) \)번 호출하는 것으로 FGCD를 계산할 수 있습니다.

 

1. \( M_{HGCD}(a,b) \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} c \\ d \end{pmatrix} \)라고 정의합시다.

2. \( M_{HGCD}(c,d) \begin{pmatrix} c \\ d \end{pmatrix} = \begin{pmatrix} e \\ f \end{pmatrix} \)라고 정의합시다.

...

3. 과정을 \( \mathcal{O}( \log k ) \)번 반복하면 \( \begin{pmatrix} g \\ 0 \end{pmatrix} \) 꼴이 됩니다. 이제 등장한 모든 \( \mathcal{O}( \log k ) \)개의 행렬을 전부 곱하면 구하고자 하는 FGCD를 계산할 수 있습니다.

 

하지만 애석하게도, 이러한 저희의 노력은 수포로 돌아갔습니다(...)

왜냐하면, 어짜피 \(k\)가 \( \mathcal{O}( \log \max(a,b) ) \)밖에 안되기 때문에 \(k\)개의 행렬을 일일히 곱해도 시간이 그렇게 느리지 않기 때문입니다.

 

그런데 이렇게 아무짝에도 쓸모없는 개념이라면 블로그 글로 적지는 않았겠죠?

 

3) Half GCD, Full GCD on Polynomial

 

사실 유클리드 알고리즘은 정수 위에서만 사용할 수 있는 알고리즘이 아닙니다. Euclidean Domain이라 불리는 공간 위에서라면 항상 유클리드 알고리즘을 적용할 수 있습니다. 대표적인 Euclidean Domain이 바로 Polynomial Ring으로, 만약 모든 계수들이 Field의 원소라면 다음 나머지 정리가 성립합니다.

 

\(f(x),g(x)\)가 존재할 때, \( f(x) = q(x)g(x) + r(x), \deg r < \deg g \)를 만족하는 다항식 \(q(x),r(x)\)가 존재합니다.

 

유클리드 알고리즘은 나머지 정리를 반복해서 사용하는것 뿐이므로, 두 다항식 간에 유클리드 알고리즘이 존재하며, 따라서 \(f(x),g(x)\)의 최대공약수인 다항식 \(h(x)\)역시 존재하며, 계산이 가능합니다. 하지만, 이때 큰 문제가 생기게 되는데, \( \deg g = n\)이라고 정의하면,

 

1. \( q(x),r(x) \)를 계산하는데 \( \mathcal{O}(n \log n) \)이 소요됩니다. \( \mathcal{O}(1) \)이 소모되던 정수에서의 유클리드 알고리즘과 상황이 달라져서, 한번의 나머지 정리 조차 엄청난 시간이 소요됩니다.

2. 최악의 경우 나머지 정리를 최대 \(n\)번 사용해야 합니다. 즉, \(k=n\)이 된다는 뜻으로, 만약 \(g\)가 \(1000000 \)차 다항식이라고 한다면, 나머지 정리를 \(1000000 \)번 사용해야 합니다. 아무리 큰수가 들어와도 최대 \(100\)번이면 결과를 얻어낼 수 있는 64bit 정수에서의 유클리드 알고리즘과 다릅니다.

 

이처럼 두 다항식의 최대공약수를 계산하는 것은 \( \mathcal{O}(n \log n) \)가 소요되는 나머지 정리를 최대 \( n\)번 사용해야 하기에 \( \mathcal{O}(n^2 \log n) \)의 시간이 소모됩니다. \(n\)이 10만을 넘어가는 경우도 매우 빈번하기 때문에, 정수에서 하듯 유클리드 알고리즘을 적용하면 최대공약수를 절대 구할 수 없습니다. 그렇기에, 빠르게 Full GCD를 계산하기 위하여, Half GCD라는 개념이 도입되게 됩니다.

 

Half GCD와 Full GCD를 더 엄밀하게 정의해봅시다. 두 다항식 \(f(x),g(x)\)가 주어질 때,

 

0. 일반성을 잃지않고 \( \deg f \geq \deg g\)라고 가정합니다.

1. \( M_{HGCD}(f,g) = \begin{pmatrix} m_{11}(x) & m_{12}(x) \\ m_{21}(x) & m_{22}(x) \end{pmatrix} \)를 다음을 만족하는 행렬 \(M\)으로 정의합니다: \( M \begin{pmatrix} f(x) \\ g(x) \end{pmatrix} = \begin{pmatrix} a(x) \\ b(x) \end{pmatrix} \)라고 표현할 때, \( \frac{ \deg g}{2} > \deg b\)가 성립한다.

2. \( M_{FGCD}(f,g) = \begin{pmatrix} m_{11}(x) & m_{12}(x) \\ m_{21}(x) & m_{22}(x) \end{pmatrix} \)를 다음을 만족하는 행렬 \(M\)으로 정의합니다: \( M \begin{pmatrix} f(x) \\ g(x) \end{pmatrix} = \begin{pmatrix} a(x) \\ b(x) \end{pmatrix} \)라고 표현할 때, \( b = 0\)가 성립하며, \( \deg a > 0 \)이면서 가능한한 \( \deg \)가 작아야 한다.

 

\( M_{HGCD} \)의 정의가 well-defined 되어있지 않음에 유의해주세요. 즉, 같은 \(f,g\)를 넣어도 결과값인 \(M_{HGCD}(f,g)\)가 달라질 수 있습니다.

 

4) Half GCD Algorithm \( \mathcal{O}(n \log^2 n) \)

 

\( M_{HGCD} \)를 어떻게 계산할 수 있을까요? \(n = \deg g, m = \lceil \frac{n}{2} \rceil \) 라고 정의하겠습니다. 이때, \(f(x),g(x)\)를 다음과 같이 분리해봅시다.

 

\(f(x) = f_0(x) + x^m f_1(x), \deg f_0 < m \)

\(g(x) = g_0(x) + x^m g_1(x), \deg g_0 < m \)

 

이때, \(n,m\)의 정의상 \(n \leq 2m\)이 성립하기 때문에 \( \deg g_1 \leq m \)이 성립합니다.

 

재귀적으로 생각하여, \( M_{HGCD}(f_1,g_1) \)을 이미 알고있다고 가정합시다. 이 행렬을 \(M^1\)으로 정의하고, 다음이 성립한다고 가정합시다.

 

\( M^1 \begin{pmatrix} f_1(x) \\ g_1(x) \end{pmatrix} = \begin{pmatrix} s_1(x) \\ t_1(x) \end{pmatrix} \), \( \frac{ \deg g_1}{2} > \deg t_1 \)

 

그러면 다음이 성립합니다.

 

\( M^1 \begin{pmatrix} f(x) \\ g(x) \end{pmatrix} = \begin{pmatrix} (M^1_{11}(x) f_0(x) + M^1_{12}(x) g_0(x)) + x^m s_1(x) \\ (M^1_{21}(x) f_0(x) + M^1_{22}(x) g_0(x))  + x^m t_1(x) \end{pmatrix} = \begin{pmatrix} s_2(x) \\ t_2(x) \end{pmatrix} \)

 

만약 운이 아주 좋아서 \( t_2(x) = 0 \) 이 성립한다고 가정합시다. 그러면, \( M_{HGCD} = M^1 \)로 정의하면 충분하기 때문에, \( M^1 \)을 return해주고 알고리즘을 종료해줍시다.

만약, \( t_2(x) \neq 0\)이 성립한다고 가정합시다. 그렇다면, 여기서 나머지 정리를 사용할 수 있습니다.

 

\( s_2(x) = q(x) t_2(x) + r(x), \deg r < \deg t_2 \)가 성립하는 \(q(x),r(x)\)가 존재합니다.

 

따라서, \( \begin{pmatrix} 0 & 1 \\ 1 & -q(x) \end{pmatrix} \begin{pmatrix} s_2(x) \\ t_2(x) \end{pmatrix}  = \begin{pmatrix} t_2(x) \\ r(x) \end{pmatrix} \) 가 성립합니다.

 

만약 운이 아주 좋아서 \( r(x) = 0 \)이 성립한다고 가정합시다. 그러면, \( \begin{pmatrix} 0 & 1 \\ 1 & -q(x) \end{pmatrix} M^1 \)을 return해주고 알고리즘을 종료해줍시다.

 

만약 \(r(x) \neq 0\)이면 어떻게 할까요? HGCD를 구하기 위해 처음 했던 행동을 그대로 반복하면 됩니다.

 

적당한 \( \alpha \) (스포일러 하자면, \( \alpha = 2m - \deg r \)입니다.) 를 설정하고, 다음과 같이 다항식을 분할합니다.

 

\( t_2(x) = t_3(x) + x^\alpha t_4(x) \)

\( r(x) = r_3(x) + x^\alpha r_4(x) \)

 

\( M^2 = M_{HGCD}(t_4,r_4) \)로 정의합시다. \( M^2 \begin{pmatrix} t_4(x) \\ r_4(x) \end{pmatrix} = \begin{pmatrix} t_5(x) \\ r_5(x) \end{pmatrix} \)로 정의할 때, 다음이 성립합니다.

 

\( M^2 \begin{pmatrix} t_2(x) \\ r(x) \end{pmatrix} = \begin{pmatrix} (M^2_{11}(x) t_3(x) + M^2_{12}(x) r_3(x)) + x^\alpha t_5(x) \\ (M^2_{21}(x) t_3(x) + M^2_{22}(x) r_3(x))  + x^\alpha r_5(x) \end{pmatrix} \)

 

저희의 목적은 \( \deg (M^2_{21}(x) t_3(x) + M^2_{22}(x) r_3(x))  + x^\alpha r_5(x) < m \)으로 만드는 것입니다. 만약, 이것이 성립한다면, \( m-1 < \frac{n}{2} = \frac{\deg g}{2} \)가 성립하기 때문에 \( M_{HGCD}(f,g) = M^2 \begin{pmatrix} 0 & 1 \\ 1 & -q(x) \end{pmatrix} M^1 \)를 return해주면 됩니다. 하지만 어떻게 저 다항식의 degree를 \(m\)미만으로 만들 수 있을까요? 저희가 가지고 있는 정보를 정리해봅시다.

 

\[ \begin{align*} &M^2_{21}(x) t_4(x) + M^2_{22}(x) r_4(x)  = r_5(x) \Rightarrow \\ &\deg M^2_{21}(x) \leq \deg r_5 - \deg t_4 \leq \deg r_5 , \\ &\deg M^2_{22}(x) \leq \deg r_5 - \deg r_4 \leq \deg r_5  \end{align*} \]

\[ \deg r_5 < \frac{ \deg r_4}{2} \]

\[ r(x) = r_3(x) + x^\alpha r_4(x) \Rightarrow \deg r_4 \leq \deg r - \alpha \]

\[ \deg t_3, \deg r_3 < \alpha \]

 

따라서, \( \deg (M^2_{21}(x) t_3(x) + M^2_{22}(x) r_3(x))  + x^\alpha r_5(x) < \frac{ \deg r - \alpha}{2} + \alpha = \frac{ \deg r + \alpha}{2} \)가 성립합니다. 그러므로, \( \alpha = 2m - \deg r \)로 설정하면 조건을 만족시킬 수 있습니다.

 

이제 마지막으로 확인해주어야 할것은 \( \deg r_4 \)입니다. 만약, 이 값이 아주 크다면(예를 들어, \(n-1\)이라고 가정하면 어떨까요) 기껏 이렇게 분할정복을 하더라도 시간은 엄청난 상태 그대로일것입니다.

 

다음이 성립합니다.

\[ \deg r < \deg t_2  = \deg (M^1_{21}(x) f_0(x) + M^1_{22}(x) g_0(x))  + x^m t_1(x) ) \]

\[ \begin{align*} &M^1_{21}(x) f_1(x) + M^1_{22}(x) g_1(x) = t_1 \Rightarrow \\ &\deg M^1_{21} \leq \deg t_1 - \deg f_1 \leq \deg t_1 \\ & \deg M^1_{22} \leq \deg t_1 - \deg g_1 \leq \deg t_1 \end{align*} \]

\[ \deg t_1 < \frac{ \deg g_1}{2} \leq \frac{m}{2} \]

\[ \deg f_0, \deg g_0 < m \]

\[ \deg t_2  = \deg (M^1_{21}(x) f_0(x) + M^1_{22}(x) g_0(x))  + x^m t_1(x) ) < \frac{3}{2}m \]

\[ \deg r_4 = \deg r - \alpha = 2 \deg r - 2m < 2 \deg t_2 - 2m < 3m - 2m = m \]

따라서 \( \deg r_4 < m\)이 성립합니다.

 

진짜진짜 마지막으로, Half GCD의 시간복잡도를 계산해봅시다. \(n = \deg g \)일 때, \( f(x),g(x) \)의 Half GCD를 구하는데 걸리는 연산수를 \(T(n)\)이라고 정의하겠습니다.

이때, 재귀적으로 부르는 HGCD는 \( M_{HGCD}(f_1,g_1), M_{HGCD}(t_4,r_4) \)가 존재하며, \( \deg g_1, \deg r_4 < m\)이기 때문에, \(\deg g_1, \deg r_4 < \frac{n}{2}\)가 성립합니다. 그러므로, 대략 이 작은 Half GCD 2개를 연산하는데 \(2T(\frac{n}{2})\)가 소요됩니다. 그리고 중간에 \( s_2(x) = q(x) t_2(x) + r(x), \deg r < \deg t_2 \) 부분을 계산할때 몫과 나머지를 계산하는 부분이 필요한데, 이는 https://hyperbolic.tistory.com/4에 적힌 방법을 이용하면 \( \mathcal{O}(n \log n) \)이 소요됩니다.

 

따라서, 연산량은 \( T(n) = 2T(\frac{n}{2}) + \mathcal{O}(n \log n) \)이 성립하며, 따라서 \( T(n) = \mathcal{O}(n \log^2 n) \)이 성립합니다.

 

5) Full GCD algorithm \( \mathcal{O}(n \log^2 n) \)

 

Half GCD를 계산하는 방법에 대해서 이미 알고있다고 가정합니다.

 

1. \( M^1 = M_{HGCD}(f,g)\)을 계산한 뒤, \( \begin{pmatrix} f_1(x) \\ g_1(x) \end{pmatrix} = M^1 \begin{pmatrix} f(x) \\ g(x) \end{pmatrix} \)를 계산합니다. 만약, \( g_1(x) = 0 \)이라면 \( M_{FGCD}(f,g) = M^1 \)을 return해주면 됩니다.

2. \( g_1(x) \neq 0 \)이라고 하면, 나머지 정리를 적용합니다. \(f_1(x)  = q(x)g_1(x) + r(x) \)라고 하면, \( \begin{pmatrix} 0 & 1 \\ 1 & -q(x) \end{pmatrix} M^1 \begin{pmatrix} f(x) \\ g(x) \end{pmatrix} = \begin{pmatrix} g_1(x) \\ r(x) \end{pmatrix} \) 가 성립합니다. 만약 \(r(x) = 0\)이라면 \( M_{FGCD}(f,g) = \begin{pmatrix} 0 & 1 \\ 1 & -q(x) \end{pmatrix}  M^1 \)을 return해주면 됩니다.

3. 만약 \(r(x) \neq 0\)라면, \( M^2 = M_{FGCD}(g_1,r) \)을 계산합니다. 그러면, \( M_{FGCD}(f,g) = M^2 \begin{pmatrix} 0 & 1 \\ 1 & -q(x) \end{pmatrix}  M^1 \)가 성립합니다.

 

시간복잡도는 \( \mathcal{O}(n \log^2 n) \)이며, 증명은 생략합니다.

 

6) Inv of Polynomial(2) \( \mathcal{O}(n \log^2 n) \)

 

\(f(x),h(x)\)가 주어질때, mod \(h(x)\) 위에서 \(f(x)\)의 inverse \(g(x)\)를 찾는 방법입니다. 즉, \(f(x)g(x) \equiv 1 \) mod \(h(x)\)를 만족하는 \(g(x)\)를 찾아야합니다. \( h(x) = x^T \) 꼴인 경우, Newton Method를 사용해서 \( \mathcal{O}(n \log n) \)에 할 수 있다는 것을 전 게시물에 설명하였습니다. 하지만 \(h(x) = x^T\)꼴이 아닌 경우, 특수한 방법이 필요합니다.

 

\(M_{FGCD}(f,h) = \begin{pmatrix} m_{11}(x) & m_{12}(x) \\ m_{21}(x) & m_{22}(x) \end{pmatrix} \)를 계산합니다. 만약, \( m_{11}(x) f(x) + m_{12}(x) h(x) = 1 \)이 성립한다면, 답은 \(m_{11}(x)\)가 되며, 그 이외의 경우는 inverse가 존재하지 않습니다.

 

마치며

 

출처는 https://www.csd.uwo.ca/~mmorenom/CS424/Lectures/FastDivisionAndGcd.html/node6.html 입니다.