[Algorithm] 소수 판별 - 밀러-라빈 판정법

* Python을 기준으로 합니다

---

소수 판별 - 밀러-라빈 판정법 (Miller-Rabin Primality Test)

- 개념

• 밀러-라빈 소수 판정법
  • 간단한 소수 판별은 $O(\sqrt{n})$의 time complexity를 가지고, 에라토스테네스의 체를 사용하면 $O(n \log \log n)$으로 개선 가능 하지만, 여전히 상당히 큰 수에 대해서는 소수 판별을 적용하기 어려움
  • 밀러-라빈 판정법은 기존 소수 판별법들과 달리 결정론적으로 동작하지 않고 확률적으로 동작하지만, 상당히 빠른 time complexity를 달성할 수 있음
    - 이때 $n < 2^{64}$의 모든 수들에 대해 $a=2,3,5,7,11,13,17,19,23,29,31,37$에 대해서만 밀러-라빈 판정을 적용하는 것 만으로도 결정론적으로 소수를 판별할 수 있다는 것이 알려져 있음
  • 이론적 유도
    1. 먼저 모든 짝수는 2를 제외하면 소수가 아니므로 $n$을 홀수인 소수라고 가정하자
    2. 그러면 $n-1$을 정수 $s$에 대해 짝수 $2^{s}$와 홀수 $d$의 곱으로 나타낼 수 있다:
       (Eq. 1) $n-1=2^{s}d$ 
    3. 이때 어떤 $a\in (\mathbb{Z}/n\mathbb{Z})$에 대해 $n$이 소수라고 가정하면, 다음 중 하나가 성립한다:
       (Eq. 2) $a^{d}\equiv 1 \, (\text{mod } n)$
       (Eq. 3) $a^{2^{r}d} \equiv (\text{mod } n) \,\,\, \text{for some} \, 0\leq r \leq s-1$
    4. 그러면 페르마의 소정리에 의해 다음이 성립한다:
       (Eq. 4) $a^{n-1}\equiv a \, (\text{mod } n)$
    5. 즉, $a^{n-1}$의 제곱근에 대한 $n$ 나머지를 재귀적으로 계속 계산했을 때, 그 결과는 항상 1 또는 -1이 된다
       - 이때 -1인 경우, (Eq. 3)가 성립한다는 것을 보일 수 있으므로 $n$은 소수이다
       - 1인 경우, (Eq. 3)이 성립하지 않으면 (Eq. 2)는 반드시 성립한다
    6. 결과적으로 밀러-라빈 판정법은 다음의 식이 모두 성립하지 않은 경우, $n$을 소수일 것 같다(probable prime)로 판별한다:
       (Eq. 5) $a^{d}\not\equiv1 \,(\text{mod } n)$
       (Eq. 6) $a^{2^{r}d} \not\equiv -1 \, (\text{mod } n) \,\,\, \forall 0\leq r\leq s-1$ 
  • 밀러-라빈 판정법의 동작 
    1. 임의의 수 $a$를 선정
       - $n < 2^{64}$의 모든 수들에 대해 $a=2,3,5,7,11,13,17,19,23,29,31,37$만 검사해도 충분함
    2. $x = a^{d} \mod n$을 계산
    3. $x$가 1 or -1 이면 `True`를 return
    4. $d$가 $n-1$이 아닌 동안 다음을 반복
       - $x = x^{2} \mod n$를 계산
       - 만약 $x$가 1이면 `False` / $x$가 $n-1$이면 `True`

밀러-라빈 판정법

- 구현

1. 밀러-라빈 판정법의 구현

• 모듈러 지수 연산을 활용하면 $O(\log^{3}n)$의 time complexity로 단축할 수 있음

#모듈러 거듭제곱
def power(x, y, p):
    res = 1

    x = x % p
    while y > 0:
        if y & 1: #y가 홀수면 
            res = res * x % p #x와 곱하고 나머지 연산

        #y가 짝수이면
        y = y>>1 #y /= 2
        x = (x * x) % p #거듭제곱 계산
     
    return res

#밀러-라빈 판정
def miller_rabin(n, a):
    d = n - 1 #홀수
    while d % 2 == 0:
        d //= 2 #2^d*r+1

    x = pow(a, d, n) #a^{d} % n 계산

    if x == 1 or x == n-1: #x가 1 또는 n-1이면 소수
        return True
    
    while d != n-1: #n-1이 아닌 동안 
        x = pow(x,2,n) #x^{2} % n 계산
        d *= 2

        #x가 1이면 합성수이고 아니면 소수
        if x == 1:
          return False
        if x == n-1:
          return True 
    
    return False


def isPrime(n):
    if n in [2,3]: #기본적으로 소수인 경우
        return True

    if n == 1 or n%2 == 0: #1과 2의 배수는 소수가 아님
        return False
    
    for a in [2,3]: #n < 1,373,653의 수에 대해서는 a = 2, 3만 검사해도 충분함
        #밀러-라빈이 성립하면 소수
        #성립하지 않아 false를 반환하면 소수가 아님 (not False = True)
        if not miller_rabin(n,a): 
            return False
 
    return True;

#-----
primes = []
for n in range(1, 100):
    if isPrime(n):
        primes.append(n)

print(primes) 
#[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]