밀러 라빈 소수 판별법(Miller-Rabin primality test) 배우기

1. 밀러 라빈 소수 판별법(Miller-Rabin primality test)

 

어떤 홀수 n이 소수인지 확률적으로 판별해주는 알고리즘

 

확률적이라는 말은, 합성수를 소수로 잘못 판별할 수 있다는 뜻이다.

 

주어진 n이 "합성수이다" 또는 "아마도 소수일 것이다"라는 대답을 내놓는다는 의미

 

 

2. 핵심 아이디어

 

알고리즘의 아이디어는 페르마의 소정리에서 시작한다.

 

https://deepdata.tistory.com/573

페르마의 소정리 문제 풀어보면서 익히기

 

 

p가 소수이고, a가 p의 배수가 아니면, $$a^{p-1} \equiv 1 (mod p)$$가 성립한다.

 

n이 홀수라고 했으므로, n-1은 짝수이다.

 

그러므로 적당한 홀수 d와 정수 s에 대하여.. 홀수와 짝수의 곱은 짝수니까.. $$n-1 = d \times 2^{s}$$

 

(근데 s = 0이면 d는 짝수네...)

 

만약 n이 정말 소수라면, n보다 작은 임의의 양의 정수 a에 대하여..

 

$$a^{n-1} = a^{d \times 2^{s}} \equiv 1 (mod n)$$

 

합동식의 양변에 1을 빼도, 성립하므로 $$a^{d \times 2^{s}} -1 \equiv 0 (mod n)$$

 

이 말은 $a^{d \times 2^{s}} - 1$은 n의 배수이거나, $a^{d \times 2^{s}} = 1$이라는 의미다.

 

여기서 $a^{d \times 2^{s}} - 1$을 인수분해하면... $(a+1)(a-1) = (a^{2} - 1)$을 이용해서...

 

$$a^{d \times 2^{s}} - 1 = (a^{d \times 2^{s-1}} + 1)(a^{d \times 2^{s-2}} + 1) ... (a^{d} + 1)(a^{d} - 1)$$

 

좌변이 n으로 나눈 나머지가 0이므로... n이 정말로 소수라면.

 

적당한 정수 r = 0,1,2,...,s-1에 대하여.. $a^{d \times 2^{r}} +1 \equiv 0 (mod n)$이거나,

 

$a^{d} \equiv 1(mod n)$이어야 한다.

 

따라서, 적당한 정수 a < n을 골랐을때, 우변의 인수 s+1개 중 n의 배수가 존재한다면 n은 아마도 소수일 것이라고 유추할 수 있다.

 

합성수를 소수로 잘못 판단할 확률을 줄이고 싶다면, 최대한 많은 a에 대해 시험해보면 된다.

 

 

3. 예시로 이해하기

 

n = 221을 예로 들어 생각해보자.

 

n - 1 = 220이고, 적당한 홀수 d와 정수 s에 대하여.. d = 55, s = 2로 한다면 55*4 = 220이다.

 

만약 a = 174를 선택한다면.. 생각해야될 인수 s+1 = 3개는..

 

$174^{55 \times 2} + 1$, $174^{55} + 1$, $174^{55} - 1$로 3개이다.

 

실제로 3개 각각을 221로 나눈 나머지를 구해보면 0,48,46이 나옴

 

인수 중 적어도 하나가 221의 배수이므로, 221은 아마도 소수일 것이다.

 

이번엔 a = 137을 선택한다면..

 

$137^{55 \times 2} + 1$, $137^{55} + 1$, $137^{55} - 1$로 3개

 

각각 나머지를 구해본다면? 206,189,187이 나온다

 

인수중 적어도 하나가 221로 나눈 나머지가 0인 경우가 없으므로 221은 확실히 합성수이다. 

 

실제로 221은 13과 17의 곱으로 합성수이다.

 

----------------------------------------------------------------------------------------------------------------------------

 

4. 몇가지 중요한 사실

 

n이 정말로 소수라면, a를 어떻게 선택하더라도 s+1개의 인수중 적어도 하나는 n의 배수일 것이다.

 

그런데 n이 합성수라면 위의 221처럼 적어도 하나가 n의 배수일 수도 있지만, 아닐수도 있다.

 

좋은 알고리즘이 될려면, n이 합성수일때, 높은 확률로 "n이 합성수이다"라는 결론을 내려야한다.

 

k번의 테스트를 수행했을때, 아주 높은 확률로, s+1개의 인수가 모두 n의 배수가 아닌 a 값을 1번 이상 뽑아야한다.

 

그러니까 단 1번이라도 통과 못하면 n은 확실하게 합성수

 

실제 알려진 결론이 다음과 같다고 한다.

 

1) "n이 합성수인 경우, 's+1개의 인수 모두가 n의 배수가 아닌 a'는 적어도 3n/4개가 존재한다"

 

따라서 n이 합성수라면, 1번의 밀러 라빈 테스트를 통과할 확률이 1/4이하라는 소리다.

 

k번의 독립 테스트를 수행했을때, 모든 테스트를 통과할 확률은 $(1/4)^{k}$이하이다.

 

k = 20정도라면 매우 높은 확률로 소수를 정확히 판별하는 알고리즘이 된다.

 

2) 다행히 n이 $2^{32}$보다 작으면 a = 2, 7, 61만 검사해봐도 잘못 판별하는 일이 없으며, 

 

$2^{64}$보다 작으면 a = 2, 325, 9375, 28178, 450775, 9780504, 1795265022만 넣어봐도 충분하다는 사실이 알려져있다.

 

더 작게는 a = 2이상 40이하의 소수만 따져봐도 된다고 한다..

 

참고로 위키피디아피셜 n의 범위에 따라 a를 어디를 잡아야하는지 알려진 사실들이 있다고 

 

 

 

3) n이 매우 작으면 오히려 무식하게 기본형 $O(\sqrt{n})$ 소수 알고리즘이 더 빠를 수 있다.

 

 

5. 구현 예시

 

실제 구현에서 d와 s를 어떻게 선택하는지 고민해보자면...

 

$$n-1 = d \times 2^{s}$$를 만족하는 d와 s를 골라야하는데

 

s = 0이면 d = n - 1

 

s = 1이면 d는 무엇을 선택해야할까? d/2를 고르면 된다.

 

s = 2이면 d는? d/4를 고르면 된다.

 

s = 3이면? d/8

 

...

 

그러므로, d = n-1부터 시작해서 1/2로 나눠보면서 test를 하면 될 것이다.

 

d는 짝수이기 때문에 2로 나누다가 홀수가 된다면 반복문을 탈출한다

 

d가 홀수가 된다는 의미는 무슨말인가?

 

다음 인수분해 식에서..

 

$$a^{d \times 2^{s}} - 1 = (a^{d \times 2^{s-1}} + 1)(a^{d \times 2^{s-2}} + 1) ... (a^{d} + 1)(a^{d} - 1)$$

 

 

$d \times 2^{s-1}$부터 시도해서, $a^{d} + 1$이랑 $a^{d} - 1$까지 도착했다는 의미다.

 

그러므로, 반복문을 탈출하고 power(a,d,n)을 구해보는데..

 

이것이 n - 1이라는 의미는? $a^{d} + 1$이 n의 배수(n으로 나눈 나머지가 0)

 

이것이 1이라는 의미는? $a^{d} - 1$이 n의 배수라는 의미(n으로 나눈 나머지가 0)

 

---------------------------------------------------------------------------------------------------------------------

 

다음은 n이 $2^{32}$보다 작을때, 소수를 확실히 판별하는 알고리즘

 

실제로 에라토스테네스의 체로 13222123은 소수이고 11221222는 소수가 아니더라..

 

파이썬은 오버플로우가 없다보니 확실히 유리하다..

 

참고로 직접 구현한 power() 대신에 파이썬에서는 x의 y승을 m으로 나눈 나머지를 pow(x,y,m)함수로도 구할 수 있다.

 

a_list = [2,7,61]

#modulo exponentiation

def power(x,y,m):
    
    result = 1

    x = x % m

    while y > 0:
        
        if y & 1:
            
            result *= x
            result %= m
        
        y = y >> 1
        x *= x
        x %= m
    
    return result

#miller_rabin test
def miller_rabin(n, a):
    
    d = n - 1

    while (d % 2 == 0):

        if power(a, d, n) == n-1:

            return True

        d = d >> 1

    temp = power(a, d, n)

    return temp == n-1 or temp == 1 #a^d = 1(mod n)인 경우

#primality test
def is_prime(n):

    if n <= 1:

        return False

    if n <= 10000: #n is very small
        
        for i in range(2,int((n+1)**(1/2))+1):
            
            if n % i == 0:
                
                return False
        
        return True
    
    for a in a_list:
        
        if miller_rabin(n,a) == False:
            
            return False
    
    return True


print(is_prime(13222123)) #True
print(is_prime(11221222)) #False

 

일반적으로 알려진 n의 범위에 따른 a값들에 대해 소수임을 확실히 판별하는 알고리즘

 

#modulo exponentiation

def power(x,y,m):
    
    result = 1

    x = x % m

    while y > 0:
        
        if y & 1:
            
            result *= x
            result %= m
        
        y = y >> 1
        x *= x
        x %= m
    
    return result
    
def miller_rabin(n, a):
    
    d = n - 1

    while (d % 2 == 0):

        if power(a, d, n) == n-1:

            return True

        d = d >> 1

    temp = power(a, d, n)

    return temp == n-1 or temp == 1

def miller_rabin_test(n,a_list):
    
    for a in a_list:
        
        if miller_rabin(n,a) == False:
            
            return False
    
    return True

def is_prime(n):

    if n <= 1:

        return False

    if n <= 1000: #n is very small
        
        for i in range(2,int((n+1)**(1/2))+1):
            
            if n % i == 0:
                
                return False
        
        return True
    
    elif n < 2047:
        
        a_list = [2]
    
    elif n < 1373653:
        
        a_list = [2,3]
    
    elif n < 9080191:
        
        a_list = [31,73]
    
    elif n < 25326001:
        
        a_list = [2,3,5]
    
    elif n < 3215031751:
        
        a_list = [2,3,5,7]
    
    elif n < 4759123141:
        
        a_list = [2,7,61]
    
    elif n < 1122004669633:
        
        a_list = [2,13,23,1662803]
    
    elif n < 2152302898747:
        
        a_list = [2,3,5,7,11]
    
    elif n < 3474749660383:
        
        a_list = [2,3,5,7,11,13]
    
    elif n < 341550071728321:
        
        a_list = [2,3,5,7,11,13,17]
    
    elif n < 3825123056546413051:
        
        a_list = [2,3,5,7,11,13,17,19,23]
    
    elif n < 18446744073709551616:
        
        # n < 2**64
        
        #a_list = [2,3,5,7,11,13,17,19,23,29,31,37]

        a_list = [2,325,9375,28178,450775, 9780504, 1795265022]
    
    elif n < 318665857834031151167461:

        a_list = [2,3,5,7,11,13,17,19,23,29,31,37]
    
    elif n < 3317044064679887385961981:
        
        a_list = [2,3,5,7,11,13,17,19,23,29,31,37,41]

    return miller_rabin_test(n,a_list)


print(is_prime(13222123))
print(is_prime(11221222))

 

 

6. 연습문제

 

5615번: 아파트 임대 (acmicpc.net)

5615번: 아파트 임대

 

 

7. 풀이

 

넓이 2xy + x + y를 인수분해한다면..

 

2xy + x + y = k라고 할때, x(2y+1) + y = k인데... 양변에 2를 곱하면..

 

2x(2y+1) + 2y = 2k

 

양변에 1을 더하면.. 2x(2y+1) + 2y + 1 = 2k + 1이므로... 2k+1 = (2x+1)(2y+1)이다.

 

x,y가 자연수이므로, 2k+1은 소수가 아니다.

 

따라서 입력으로 자연수 k가 주어질텐데.. 홀수 2k+1이 소수인지 아닌지 판별하면 된다.

 

소수가 아니라면 자연수 x,y에 대해 (2x+1)(2y+1)인 x,y가 존재하고 소수라면 자연수 x,y가 존재하지 않는다.

 

입력이 2^31 - 1이하이므로, a = [2,7,61]에 대해서 밀러 라빈 소수 판별법으로 빠르게 검사하면 된다

 

from sys import stdin

def miller_rabin(n, a):
    
    d = n - 1

    while d % 2 == 0:
        
        if pow(a,d,n) == n-1:
            
            return True
        
        d = d >> 1
    
    temp = pow(a,d,n)

    return temp == n-1 or temp == 1

def is_prime(n,a_list):
    
    if n <= 1:
        
        return False
    
    elif n <= 10000:
        
        for i in range(2,int((n+1)**(1/2))+1):
            
            if n % i == 0:
                
                return False
        
        return True
    
    else:
        
        for a in a_list:
            
            if miller_rabin(n,a) == False:
                
                return False
        
        return True

a_list = [2,7,61]

n = int(stdin.readline())

answer = 0

for _ in range(n):
    
    area = int(stdin.readline())

    p = 2*area+1

    answer += is_prime(p,a_list)

print(answer)

 

참조

 

 

https://casterian.net/algo/miller-rabin.html

casterian.net

 

 

https://rkm0959.tistory.com/183

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

 

https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test

Miller–Rabin primality test - Wikipedia