6. Miller-Rabin 소수 판별 알고리즘과 Pollard-Rho 소인수분해

관련 코드는 github에서 찾아볼 수 있다.

이 글에서 다룰 내용의 많은 부분에서 완벽한 수학적 논의는 생략할 것이다.

Miller-Rabin 소수 판별법

이 부분에 대해서는 이미 좋은 자료가 있습니다. https://casterian.net/archives/396

우리는 주어진 $n$이 소수임을 확인하고자 한다. 먼저 $n-1 = d \cdot 2^s$인 정수 $s, d$를 찾는다.

만약 $n$이 소수라면, 각 자연수 $1 \le a < n$에 대하여 다음이 성립한다.

• [P] $a^d \equiv 1 \pmod{n}$이거나, $r=0, 1, \cdots , s-1$ 중 적어도 하나에 대해 $a^{2^rd} \equiv -1 \pmod{n}$

증명을 해보자. $a^{n-1}-1 = a^{d \cdot 2^s}-1 = (a^d-1)(a^d+1)(a^{2d}+1)(a^{2^2d}+1) \cdots (a^{2^{s-1}d}+1)$이다.

그런데 페르마 소정리에 의하여 좌변은 $0 \pmod{n}$이므로, 우변 역시 $n$의 배수임을 알 수 있다.

$n$이 소수이므로, 여러 자연수의 곱이 $n$의 배수라는 것은 그 중 하나 이상이 $n$의 배수임을 의미한다. 

그러므로, $a^d \equiv 1 \pmod{n}$이거나 적당한 $0 \le r <s$가 있어 $a^{2^rd} + 1 \equiv 0 \pmod{n}$이 성립한다.

$n$이 주어졌을 때, $s, d$의 값을 계산하는 것은 쉽게 할 수 있다. 

또한, $1 \le a < n$을 만족하는 $a$가 주어졌을 때, [P]가 성립하는지 여부도 역시 빠르게 확인할 수 있다. 

이제 다음과 같은 알고리즘을 하나 고안할 수 있다.

• 테스트 횟수 $k$를 하나 정한다. $k$를 어떻게 잡는가에 대한 이야기는 나중에 한다.
• $1 \le a < n$을 만족하는 $a$를 랜덤하게 하나 잡자.
• [P]가 성립하는지 확인하자. 만약 성립하지 않는다면, $n$이 소수가 아니라는 뜻이고 알고리즘 종료.
• [P]가 성립한다면, 다시 처음으로 돌아가서 알고리즘을 반복한다. $k$번에 걸쳐 항상 [P]가 성립한다면 소수라고 판단하고 종료.

만약 $n$이 소수라면, $a$를 어떻게 잡더라도 [P]가 성립할 것이니, 무조건 소수라고 판단할 것이다.

$n$이 합성수라면 어떨까? $n$이 소수라면 [P]가 성립한다고 했지, $n$이 합성수면 무조건 [P]가 거짓이라고 하지는 않았다.

위 알고리즘이 정당한 소수 판별 알고리즘이 되려면, $n$이 합성수일 경우 높을 확률로 $n$이 소수가 아니라는 결론을 내려야 한다.

즉, $k$번의 테스트를 거쳤을 때 [P]가 거짓이 되는 $a$의 값이 매우 높은 확률로 한 번 이상 뽑혀야 한다. 이는 결국

• $n$이 합성수인 경우, [P]가 거짓이 되는 $a$가 얼마나 많은가? 

라는 문제와 동치이다. 이제 중요한 결론을 증명하지 않고 제시한다. 증명은 구글 ㄱㄱ

• $n$이 합성수인 경우, [P]가 거짓이 되는 $a$는 적어도 $3n/4$개 존재한다.
• 즉, $n$이 합성수인 경우, 한 번의 테스트를 통과할 확률은 $1/4$ 이하이다. 

이는 $k$번의 독립적인 테스트를 거쳤을 때, 합성수가 모든 테스트를 통과할 확률이 $4^{-k}$ 이하라는 뜻이다. 

그러니 $k=20$ 정도로 잡아주면 매우 높은 확률로 정답을 반환하는 좋은 소수 판별 알고리즘이 될 수 있다.

하지만 랜덤이 싫을 수 있다. 랜덤을 제거한 Miller-Rabin 알고리즘의 여러 variant가 존재한다.

• Generalized Riemann Hypothesis를 가정하면, $2 \le a < 2(\ln n)^2$를 전부 시도하면 정확하게 소수 판별을 할 수 있다.
• $n < 2^{32}$라면, $a = 2, 7, 61$인 경우만 따져보면 된다. 이때, $n = 2, 7, 61$인 경우를 따로 처리해야 함에 주의하라.
• $n < 2^{64}$라면, $a$로 $2$ 이상 $40$ 이하의 소수만 따져보면 된다. 이때, $n$ 역시 $2$ 이상 $40$ 이하의 소수인 경우에 주의하라.    

랜덤을 제거한 variant가 아마 실제 성능도 좋을 것이다. $n < 2^{64}$의 variant에서도 테스트를 12번만 하기 때문이다.

구현 노트

• 구현이 까다로운 알고리즘은 아니지만, C++의 int 범위를 넘어가는 modulus에서 여러 연산을 해야한다. 
• 특히 까다로운 것은 곱셈이다. $a, b, n$이 모두 64bit 정수일 때, $ab \pmod{n}$을 계산해야 한다.
• 이를 위해서 쓸 수 있는 첫 번째 방법은 이집트 곱셈, 즉 이진법을 사용한 곱셈이다. 이를 위해서 로그 시간이 필요하다.
• 두 번째 방법은 __int128을 동원하는 것이다. $ab$의 값을 __int128에서 계산한 뒤, $n$에 대한 나머지를 구한 후 64bit 정수로 반환한다. 
• $n$이 작으면, $\mathcal{O}(\sqrt{n})$ 알고리즘이 Miller-Rabin 알고리즘보다 빠를 수 있다. 최적화가 필요한 경우 참고하자.
• 더 적은 테스트 횟수로 $2^{64}$ 미만 자연수의 소수 여부를 판별할 수 있다. 이때, $a$ 값이 소수가 아닌 경우 "$a$의 약수인 소수"들을 특별히 고려해야 한다. 
• 위 두 가지 포인트를 추가로 고려하여 작성한 Miller-Rabin 구현체를 원한다면, 맨 위에 소개한 casterian님의 블로그를 참고하라.

Pollard-Rho 소인수분해

이 부분 역시 좋은 자료가 있습니다. https://aruz.tistory.com/140

핵심 아이디어 1. Floyd's Cycle Detection Algorithm

잠깐 전형적인 CP/PS 토픽인 그래프로 넘어오자. functional graph는 각 정점이 outdegree 1을 갖는 directed graph다.

정점 $v$가 있으면 그 정점에서 나가는 간선이 정확히 하나이니, 그 간선을 따라 도달하는 정점을 $f(v)$라 정의할 수 있다. 

즉, $f(v)$는 $v$의 "다음 정점"이다. 이제 $v \rightarrow f(v) \rightarrow f(f(v)) \rightarrow \cdots$라는 path를 생각하자.

결국 이 경로는 하나의 일직선을 그리다가, 이미 방문한 정점을 다시 방문하면서 cycle에 들어가게 될 것이다. 이 cycle을 어떻게 찾을까?

물론 여러 방법이 있겠으나, Floyd의 "토끼와 거북이" 알고리즘이 매우 효율적이다. 

• 잠시 $f^n(v) = f(f^{n-1}(v))$로 정의하자. 즉, $f^n(v)$는 $v$의 "$n$번째 뒤"의 정점이다.
• 토끼 한 마리와 거북이 한 마리를 $v$에 두자. 각 동물들의 위치를 변수 $r, t$로 명명하고 관리하겠다.
• 각 단계마다, 토끼는 두 칸, 거북이는 한 칸 이동시킨다. 즉, $r \leftarrow f(f(r))$, $t \leftarrow f(t)$를 반복한다.
• 만약 $r = t$인 시점이 오면, 사이클을 찾은 것이다. 물론, 맨 처음에 $r=t=v$인 경우는 제외하고 생각한다.

• 왜 사이클을 찾은 것인가? $k$번의 단계를 거친 후에 $r = t$임을 확인했다고 가정하자.
• 이는 $f^k(v) = f^{2k}(v)$임을 의미한다. 즉, $k$번째 뒤 정점과 $2k$번째 뒤 정점이 같다는 것을 의미한다.
• 이제 $k$번째 정점인 $f^k(v)$부터 시작하여, 다음 정점으로 이동하는 것을 반복한다. 
• 다시 시작한 $f^k(v)$에 도착할 때까지 이를 반복하면, 지금까지의 경로가 완벽하게 찾고자 하는 사이클이 된다.

• 왜 사이클이 찾아지는가? $i < j$이고 $f^i(v) = f^j(v)$인 $i, j$가 존재할 때, (즉 사이클이 있을 때) 자연수 $k$가 있어 $f^k(v) = f^{2k}(v)$임을 보이면 된다.
• "한 칸씩 격차가 벌어지는 걸 결국 사이클 길이만큼 격차가 벌어져서 다시 만났다"고 생각하면 직관적이며, 증명도 똑같다.
• 속도 차이가 나는 두 사람이 400m 트랙에서 달렸을 때, 실제로 달린 거리가 400m 차이가 나서 만날 수 있는 것과 같은 원리.
• 자세한 증명을 원한다면, $t(j-i) \ge i$인 자연수 $t$를 하나 잡고, $k=t(j-i)$라고 정의해보면 증명됨을 확인할 수 있다.

• 왜 효율적인가? 일단 시간복잡도가 (사이클 도달에 필요한 거리) + (사이클의 길이) 정도로 매우 효율적이다.
• 공간복잡도도 효율적이다. 각 정점 $v$에 대해서 $f(v)$를 계산하는데 드는 메모리를 제외하면, $\mathcal{O}(1)$만큼의 메모리가 추가적으로 필요하다.
• 특히, 여러 경우에서 $f(v)$는 진짜 간선으로 표현되는 "다음 정점"이 아니라, 말 그대로 $v$에 대한 함수가 된다. 
• 이는 아래의 경우도 마찬가지. 애초에 이름이 functional graph인 이유가 있다. 이 경우, $f(v)$를 계산하는데 드는 메모리도 $\mathcal{O}(1)$.

핵심 아이디어 2. Pseudorandom Sequence. Birthday Paradox

Birthday Paradox는 $\{1, 2, \cdots , n\}$에서 $\sqrt{n}$급 개수의 원소를 독립적으로 뽑으면, 같은 원소를 두 번 이상 뽑을 확률이 상당히 높다는 결과이다.

이는 $n=365$일 때, 대충 20~30명의 학급에서 생일이 겹치는 경우가 많다는 내용으로 알려진 결과일 것이다.

이제 우리의 목표는 주어진 합성수 $n$을 소인수분해 하는 것이다. $p$를 $n$의 가장 작은 소인수라고 하자. 우리는 $n$은 알고 있으나 $p$는 모른다.

이제 다음과 같은 방식으로 랜덤한 수열을 잡는다. $x, c$ 각각을 $1$ 이상 $n$ 미만 랜덤한 자연수라고 두자.

$x_0 = x$, $x_i = (x_{i-1}^2 + c) \pmod{n}$으로 수열을 잡자. 여기서 짚고 넘어가야 할 부분은

• 수열 $\langle x_k \rangle$은 pseudorandom. 즉, 완전 랜덤한 것은 아니지만 랜덤에 가깝다.
• $x_i$의 값은 $x_{i-1}$의 값에 의하여 완전히 결정된다. 그러니, $x_{i-1} \rightarrow x_i$를 간선으로 생각하는 functional graph를 생각할 수 있다.
• 즉, 이 수열의 값을 정점으로 갖는 functional graph를 생각할 수 있고, 마찬가지로 cycle 개념도 생각할 수 있다.
• $x_i \equiv (x_{i-1}^2 + c) \pmod{p}$ 역시 성립한다. 그러니, $y_i = x_i \pmod{p}$로 정의하면, $y_i \equiv (y_{i-1}^2 + c) \pmod{p}$이다.
• 현재 계산할 수 있는 것은 $x_i$이며, $y_i$는 실제로는 $p$를 모르기 때문에 계산할 수 없는 값이다.

이제 본격적으로 알고리즘을 설계해보자.

• $\text{gcd}(x_{2k} - x_k, n)$이라는 값을 생각해보자. 이 값을 계산하려면 $x_{2k}, x_k$의 값이 필요하며, 이는 "토끼와 거북이"의 방식으로 관리할 수 있다.
• $y_k = y_{2k}$라면, (즉 수열 $y$에서 cycle이 발견되었다면) $\text{gcd}(x_{2k} - x_k, n)$은 $p$의 배수가 된다. 
• Birthday Paradox에 의해, 이렇게 수열 $y$에서 cycle이 발견되려면 평균적으로 $\mathcal{O}(\sqrt{p})$ 정도 index가 되어야 한다. 
• 비슷하게, 수열 $x$에 대해서 cycle이 발견되려면 평균적으로 $\mathcal{O}(\sqrt{n})$ 정도 index가 되어야 한다.
• 즉, 평균적으로 $\mathcal{O}(\sqrt{p})$ 정도의 계산량 이후로 $\text{gcd}(x_{2k}-x_k, n)$이 $1$이 아닌 시점이 오며, 이때 $\text{gcd}(x_{2k}-x_k , n)$은 높은 확률로 $n$이 아니다.
• 이는 $\text{gcd}(x_{2k}-x_k, n)$이 $n$의 "nontrivial divisor"가 된다는 뜻이다. 이를 반복하면 소인수분해를 얻는다.

$n$의 non-trivial divisor를 얻기까지 평균적으로 $\mathcal{O}(\sqrt{p}) = \mathcal{O}(n^{1/4})$번의 GCD 연산이 필요하다.

구현 노트.

• 이 알고리즘에서는 random을 사용한다. C++의 rand()를 사용하기 보다는, mt19937 등을 사용하는 것을 추천한다.
• 앞에서도 언급했지만 $n$이 합성수라는 조건이 붙는다. 미리 Miller-Rabin 소수 판별을 통해 $n$이 소수인지 판별을 해야한다.
• $\text{gcd}(x_{2k}-x_k, n)=n$일 확률이 낮다고 했지, 불가능하다고 하지는 않았다. 이 경우, 시작하는 값 $x, c$를 다시 잡고 알고리즘을 다시 돌리자.
• $\text{gcd}(x_{2k}-x_k, n)=d$가 $1, n$이 아니라면, 문제를 $d$와 $n/d$를 각각 소인수분해 하는 것으로 바꿀 수 있다. 이제 재귀적으로 문제를 풀자.
• 앞에서 언급한 64bit 정수의 곱셈에 대한 내용이 여기서도 적용된다. 이집트 곱셈을 쓰거나, __int128을 사용하자.