5. 밀러-라빈 소수판정법과 폴라드-로 소인수분해

16563번 문제에서 말한 바 있듯이, 무식하게 전부 나눠 보는 방식인 에라토스테네스의 체는 시간복잡도가 \(O(N\log\log N)\)이고 조금 더 영리하게 계산한 linear seive는 시간복잡도가 \(O(N)\)이다.

이 두 가지 모두 소수를 판정하고 소인수분해하기는 확실하지만, 매우 심각한 하자가 있다. 시간이 너무 오래 걸린다는 것이다.

충분히 빠른 컴퓨터라면 대여섯자리 수가 소수인지 검사하는 것은 그렇게 오래 걸리지 않는다. 하지만 10자리라면? 20자리라면? 4149번 문제에서는 최악의 경우 19자리 수를 계산해야 한다.

영리한 linear seive를 사용하더라도 걸리는 시간이 N에 비례하기 때문에 시간초과가 날 것이다.

이를 정말 빠르게 계산할 수 있는 것이 밀러-라빈 소수판정법과 폴라드-로 소인수분해이다. 밀러-라빈은 적당히 작은 N에 대해서 시간복잡도가 \(O(\log^3 N)\)에 불과하고, 폴라드-로는 \(O(n^{1/4})\)밖에 되지 않는다.

 

먼저 밀러-라빈 소수판정법을 살펴 보자.

2 빼면 모든 짝수는 소수가 아니므로, 그 경우만 예외처리하면 소수인지 아닌지 판별할 수 N은 홀수라고 볼 수 있다. 이 N에 대해서, N-1을 짝수부와 홀수부로 분리해 보자. 이렇게 말이다.

$$N-1=2^s\cdot d\,\,\,\,(0<s\in\mathbb{Z}, 2\not\mid d)$$

만약에 N이 소수라고 치면, 페르마의 소정리에 의해서 N의 배수가 아닌 a에 대해서 \(a^{N-1}\equiv a\pmod{N}\)이 성립할 것이다.

그럼 식을 변형하여 \(\left(a^{2^{s-1}\cdot d}\right)^2\equiv 1\pmod{N}\)이 될 것이고, \(a^{2^{s-1}\cdot d}\equiv -1\pmod{N}\)이거나 \(a^{2^{s-1}\cdot d}\equiv 1\pmod{N}\)이 성립할 것이다.

후자의 경우에는 다시 똑같은 방법으로 \(a^{2^{s-2}\cdot d}\equiv -1\pmod{N}\)이거나 \(a^{2^{s-2}\cdot d}\equiv 1\pmod{N}\)이 성립할 것이고... 계속 반복한다면 다음 둘 중 하나는 성립을 하게 된다.

\(a^d\equiv1\pmod{N}\)이거나, 적당한 \(0\le r < s\)이 존재해서 \(a^{2^r\cdot d}\equiv -1 \pmod{N}\)이거나.

밀러-라빈 소수판정법은 이걸 이용해서, N이 충분히 많은 a에 대해서 위의 관계가 성립하면 N을 아마 소수이지 않을까?라고 생각해 보자는 것이다.

이를 실제로 소수인지 아닌지 판정해 보려면 적당히 많은 a에 대해서만 테스트해 보면 된다. Sorenson과 Webster는 N이 3,317,011,064,679,887,385,961,981보다 작을 때, a가 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41일 때만 검증하면 된다고 밝힌 바 있다.

이를 Python3 코드로 짜면 다음과 같다.

pli = [2,3,5,7,11,13,17,19,23,29,31,37,41]

def power(x,y,p):
    if y<2: return (x**y)%p
    else:
        d = y//2
        if y%2 == 0: return (power(x,d,p)**2)%p 
        else: return (x*(power(x,d,p))**2)%p
        
def miller_rabin(n,a):
    s,d = 0,n-1
    while d%2 == 0: s += 1; d //= 2
    x = power(a,d,n)
    if x == 1 or x+1 == n: return True
    for i in range(0, s-1):
        x = power(x,2,n)
        if x+1 == n: return True
    return False
 
def isprime(n):
    if n in pli: return True
    if n == 1 or n%2 == 0: return False
    for p in pli:
        if not miller_rabin(n,p): return False
    return True

pli는 a로 두어 줄 2부터 41까지의 소수이다. 그리고 power(x,y,p)는 \(x^y\)를 p로 나눈 나머지를 분할 정복으로 빠르게 구해 주는 함수. 18291번 문제의 함수를 개량한 것.

miller_rabin 함수가 바로 \(a^d\equiv1\pmod{N}\)이거나, 적당한 \(0\le r < s\)이 존재해서 \(a^{2^r\cdot d}\equiv -1 \pmod{N}\)이거나가 진짜로 성립하는지를 검증해 보는 것이다. 저 조건을 ★(a)라고 하자.

함수 속의 while문에서 N을 홀수 부분인 d와 짝수 부분인 \(2^s\)로 쪼개 주고, 그 다음 \(x \equiv a^d \pmod{N}\)으로 두었다.

만약 x가 1 혹은 N-1이라면 조건이 만족하는 것이므로 ★(a)가 성립하니까 True를 뱉는다. 아니라면 x를 계속 제곱해 가면서 ★(a)가 성립하는지 테스트해보는 것.

miller_rabin 함수가 True를 뱉었다면 N은 ★(a)가 성립한다는 사실을 말해주고 있는 것이다.

이제 isprime 함수 내에서, pli 속의 모든 a에 대해 miller_rabin이 성립하는지. 곧 ★(a)가 성립하는지를 검증한다. 만약 그렇다면 N이 소수인 것이고, True를 뱉어주는 것.

 

밀러-라빈만 있으면 소인수분해를 하기 충분할까? 전혀 아니다.

밀러-라빈은 어떤 수가 소수인지 아닌지만 알려주는 알고리즘이지, 소인수를 모두 구해 주는 것은 아니다. 소인수분해까지 고려를 해 보려면 폴라드-로 알고리즘을 이용해야 한다.

정확하게는 소인수분해가 아니라 인수분해이고 더 정확히 말하자면 자명하지 않은 인수(1과 N을 뺀 인수)를 하나 뱉어주는 알고리즘이지만, 어차피 인수분해를 적당히 여러 번 반복해 주면 소인수로 분해가 될 테니 거기서 거기이다.

우선 폴라드-로 알고리즘이 효과를 발휘하기 위해서는 N이 소수가 아니라는 전제조건이 필요하다.

N의 약수를 구하기 위해서 적당한 수 x를 잡아 \(d=\gcd(x,N)\)를 구할 것이고, d가 1이 아니라면 N의 자명하지 않은 인수가 될 것인데 N이 소수라면 매우매우 높은 확률로 d=1이 될 것이기 때문이다.

그러니 방금 미리 구한 밀러-라빈으로 N이 소수일 때와 짝수일 때 미리 걸러 주자.

그 다음으로는 \(a_n \equiv a_{n-1}^2 + c \pmod{N}\)이, c가 N의 약수가 아닐 경우 충분히 훌륭한 유사난수 생성 수열임을 알아야 한다.

그러니 비둘기집의 원리에 의해서, N의 약수 d에 대해서는 \(a_n \mod d\)가 원소의 수 d개 이하의 사이클을 가지게 된다.

따라서 \(a_n \mod d\)에서 \(a_{n+1} \mod d\)로 가는 화살표를 그어 본다면 예쁜 그래프가 형성된다. 이렇게 말이다.

폴라드-로 그래프

이 알고리즘을 폴라드-로 알고리즘이라고 부르는 이유가 있다. 그래프의 모양이 로(\(\rho\))를 닮았기 때문.

이제, N의 약수 d에 대해서 수열이 사이클을 이룬다면, 어떠한 두 정수 x, y가 존재해서 \(a_x\equiv a_y \pmod{d}\)이 성립을 하게 될 것. 이 말이 무엇이냐. \(\gcd(|a_x-a_y|,N)\)이 1이 아니라는 것이다.

x=y를 2에서 N-1 중 랜덤으로 하나 뽑은 값으로 두고 이걸 \(a_0\)이라고 하자. \(x=a_n, y=a_{2n}\)일 경우에, \(d=\gcd(x,y)\)로 둔다면, d가 처음으로 1이 아닌 값을 가질 때 이 d가 N의 약수 중 하나일 것이다.

만약 d가 N이 되면 c를 잘못 잡은 것이니 처음부터 다시 시작해서 c를 다른 값으로 채택하면 된다. 그렇지 않다면 성공적으로 N의 약수를 하나 잡은 것이다.

만약 d가 소인수라면? N의 소인수를 하나 잡은 것이니 빙고. 그렇지 않다면 d를 N으로 두고 처음부터 폴라드-로 알고리즘을 돌리면 된다. 그러면 d의 소인수를 뱉을 텐데, d는 N의 약수이므로 뱉어지는 값은 N의 약수이기도 할 것이다.

from random import randrange

def gcd(a,b):
    if a<b: a,b = b,a
    while b != 0: a,b = b,a%b
    return a
    
def rho(n):
    if isprime(n): return n
    if n == 1: return 1
    if n%2 == 0: return 2
    x,c,d = randrange(2,n),randrange(1,n),1
    y = x
    while d == 1:
        x = ((x**2%n)+c)%n
        y = ((y**2%n)+c)%n
        y = ((y**2%n)+c)%n
        d = gcd(n,abs(x-y))
        if d == n: return rho(n)
    if isprime(d): return d
    else: return rho(d)

위에서 논한 모든 것을 Python3 코드로 짠 결과물이 바로 이것이다.

gcd는 math 모듈을 가져올 수도 있지만, 그렇게 하지 않고 직접 유클리드 호제법으로 만드는 방법을 선택했다.

그리고 rho 함수가 바로 N의 임의의 소인수를 뱉는 코드.

 

이것을 잘 이용한다면, 막대하게 큰 수의 소인수분해도 무리 없이 단시간 내에 수행이 가능하다.

그 활용법을 4149번 문제의 해설 포스트에서 확인할 수가 있다. 아니면 큰 수의 소인수분해를 요구하는 다른 문제들이라거나.

기본적으로, BOJ에서는 다이아몬드 5 이상으로 랭크된 고난이도 문제들이다.

 

이 포스트는 이쯤에서 마치도록 하자.