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

Park Yeongseo·2024년 4월 9일
2

Algorithm

목록 보기
8/19
post-thumbnail

David.M.Burton의 Elementary Number Theory 번역 및 정리 글입니다.

소수를 판별하는 문제는 오랫동안 수학자들의 관심을 끌어온 문제지만 어렵다. 가장 흔히 쓰고 잘 알려진 방법은 "에라토스테네스의 체"를 이용하는 것으로, O(N)O(\sqrt{N})의 시간복잡도를 가진다. 문제는 long long 범위의 수의 경우, O(N)O(\sqrt{N})의 시간복잡도도 너무 크다는 것이다.

밀러-라빈 소수 판별법은 확률적으로 소수를 판별하는 알고리즘으로, O(log3N)O(\log^3 N)의 시간복잡도를 가지며, long long 범위 정도에 대해서는 결정론적으로도 소수를 판별할 수 있다.

1. 페르마의 소정리와 소수 판별법

밀러-라빈 소수 판별법에 앞서 페르마의 소정리부터 차근차근 논의를 시작해보자.

(1) 페르마의 소정리

페르마의 소정리는 다음과 같다.

페르마의 소정리
pp가 소수이고 p∤ap \not\mid a라 하면, ap11(modp)a^{p-1} \equiv 1 (\mod{p})다.

증명은 다음과 같이 할 수 있다.

aa부터 (p1)a(p-1)a까지 p1p-1개의 aa의 배수들을 생각해보자.

a,2a,3a,...,(p1)aa, 2a, 3a, ..., (p-1)a

이 수들 중 어떤 임의의 두 수를 뽑더라도 법 pp에 대해 합동인 것은 없다. 1r<sp11 \leq r < s \leq p-1를 만족하는 두 정수 r,sr, s에 대해, rasa(modp)ra \equiv sa (\mod{p})를 만족하는 것이 있다고 해보자.

rasa(modp)    p(rasa)    pa(rs)\begin{aligned} ra &\equiv sa(\mod{p})\\ \iff p &\mid (ra - sa)\\ \iff p &\mid a(r- s) \end{aligned}

pp가 소수이고 p∤ap \not\mid a이므로 p(rs)p | (r-s)여야 한다. 즉 rs(modp)r \equiv s (\mod{p})이다. 하지만 pp보다 작은 자연수 1,2,...,p11, 2, ..., p-1 중에는 법 pp에 대해 합동인 수가 있을 수 없다. 따라서 a,2a,3a,...,(p1)aa, 2a, 3a, ..., (p-1)a에서 어떤 임의의 두 수를 뽑더라도 법 pp에 대해 합동인 것은 있을 수 없다.

그렇다면 a,2a,3a,...,(p1)aa, 2a, 3a, ..., (p-1)amodp\mod {p}를 해서 얻을 수 있는 결과들의 집합{1,2,3,...,p1}\{1, 2, 3, ..., p-1 \}과도 같게 될 것이다.(둘 모두 pp 미만의 모든 자연수를 원소로 갖는 집합이므로) 따라서 다음이 성립한다.

a2a...(p1)a12...(p1)(modp)    (p1)!ap1(p1)!(modp)\begin{aligned} a \cdot 2a \cdot ... \cdot (p-1)a &\equiv 1 \cdot 2 \cdot ... (p-1) (\mod{p})\\ \iff (p-1)!a^{p-1} &\equiv (p-1)! (\mod{p}) \end{aligned}

p∤(p1)!p \not\mid (p-1)!이므로, 앞의 rasa(modp)ra \equiv sa (\mod p)에서 aa를 제거한 것과 같이, (p1)!ap1(p1)!(modp)(p-1)!a^{p-1} \equiv (p-1)! (\mod {p})에서도 (p1)!(p-1)!을 제거할 수 있다. 즉 ap11(modp)a^{p-1} \equiv 1 (\mod{p})이다. \square

Corollary

위의 정리는 아래의 따름정리로 일반화할 수 있다.

따름정리
pp가 소수이면, 정수 aa에 대해 apa(modp)a^p \equiv a (\mod{p}) 이다.

만약 pap \mid a이면, ap0a(modp)a^p \equiv 0 \equiv a (\mod{p})이므로 당연히 성립한다. 한편 p∤ap \not\mid a인 경우, 위의 정리에 따라 ap11(modp)a^{p-1} \equiv 1 (\mod{p})임을 알고 있고, 이로부터 apa(modp)a^p \equiv a (mod p)를 얻을 수 있다. \square(여기의 Theorem 1.4 - ii)

(2) 페르마의 소정리를 이용한 결정론적 소수 판별

페르마의 소정리가 가지는 함의는 해당 정리를 통과하지 못하는 모든 정수에 대해 합성수라는 결론을 내릴 수 있다는 점이다. 그렇다면 위 테스트를 통과하는 수들은 어떨까? 어떤 정수가 위 테스트를 통과했다 하더라도 그 수가 정확히 소수인지는 알 수 없지만, "아마도 소수"일지도 모른다는 결론은 내릴 수 있다. an11(modn)a^{n-1} \equiv 1 (\mod{n})을 만족하는 정수 nn을 가리켜, 밑 aa에 대한 유사소수(pseudoprime)이라 부른다. 유사소수인 합성수는 희박하지만, 무한히 많다.

aa에 제한을 두면, nn정확히 소수인지 아닌지를 판단할 수 있게 되는데, 예전에 포스팅했던 뤼카의 정리의 에두아르 뤼카가 제시한 다음의 정리가 대표적이다. 해당 정리에 들어가기 전에, 우선 위수(order)와 관련 정리에 대해 알아보자.

(2-1) 위수

위수
n>1n>1이고 gcd(a,n)=1gcd(a,n)=1일 때, 어떤 정수 aa가 법 nn에 대해 kk의 위수를 가진다는 것은,  kk 가 ak1(modn)a^k \equiv 1(\mod { n}) 을 만족하는 최소의 자연수임을 말한다.

정수 aa가 법 nn에 대해 위수 kk를 가진다고 하자. 이때, 다음이 성립한다.

정리 1
ah1(modn)    kha^h \equiv 1 (\mod{n}) \iff k|h

우선 khk|h라 하면, h=jkh = jk가 되는 정수 jj가 있다. ak1(modn)a^k \equiv 1 (\mod{n})이므로 (ak)j1j(modn)(a^k)^j \equiv \equiv 1^j(\mod{n})이고, 따라서 ah1(modn)a^h \equiv 1 (\mod{n})이다.

반대로, hhah1(modn)a^h \equiv 1 (\mod{n}) 을 만족하는 임의의 양의 정수라 해보자. h=qk+r(0r<k)h = qk + r (0 \leq r < k)를 만족하는 정수 q,rq, r이 존재하며, 다음과 같이 쓸 수 있다.

ah=aqk+r=(ak)qara^h = a^{qk + r} = (a^k)^qa^r

ah1(modn)a^h \equiv 1 (\mod{n})이라 가정했고, ak1(modn)a^k \equiv 1 (\mod{n})이므로, ar1(modn)a^r \equiv 1 (\mod{n})이다. kkak1(modn)a^k \equiv 1 (\mod{n})을 만족하는 가장 작은 양의 정수이고, 0r<k0 \leq r < k이므로 r=0r = 0이어야 한다. 따라서 h=qkh = qk이고, khk |h다. \square

다음의 정리도 염두에 두자

정리 2(오일러 정리)
a,na, n이 서로소인 양의 정수일 때,
aϕ(n)1(modn)a^{\phi(n)} \equiv 1 (\mod{n})이다.

ϕ(n)\phi(n)nn 이하의 양의 정수 중 nn과 서로소인 것의 개수이며, ϕ(n)=n1\phi(n) = n-1이라는 말은 nn과 서로소인 nn 이하 양의 정수가 n1n-1개라는 것, 즉 nn이 소수임을 가리킨다. 위 정리의 증명은 나무위키를 참고하자. 페르마 소정리 증명과 같은 방식으로 증명이 이루어지니 이해하기 쉽다. 또한 위의 정리 1과 함께, kϕ(n)k | \phi(n)임도 알 수 있다.

(2-2) 뤼카가 제시한 정리

뤼카가 제시한 정리는 다음과 같다.

뤼카
(n1)(n-1)을 나누는 모든 소수 pp에 대해, an11(modn)a^{n-1} \equiv 1 (\mod{n})이고 a(n1)/p≢1(modn)a^{(n-1)/p} \not\equiv 1 (\mod{n})인 정수 aa가 존재하면 nn은 소수다.

정수 aa가 법 nn에 대해 kk의 위수를 가진다고 하자(존재성 증명은 건너뛰었지만, 모든 정수 aa에 대해 위수가 존재함은 증명되어 있다). 정리 1에 따라 an11(modn)a^{n-1} \equiv 1 (\mod{n})이면 kn1k | n-1이다. 어떤 정수 jj에 대해 n1=kjn-1 = kj라 하자. 만약 j>1j > 1이면 jj는 소수인 약수 qq를 가지게 되며, 따라서 j=qhj = qh를 만족하는 정수 hh가 존재해 다음과 같이 된다.

a(n1)/q=(ak)h1h=1(modn)a^{(n-1)/q} = (a^k)^h \equiv 1^h = 1(\mod{n})

하지만 정리에서 가정하기를 어떤 (n1)(n-1)을 나누는 모든 소수 pp에 대해 a(n1)/p≢1(modn)a^{(n-1)/p} \not\equiv 1 (\mod{n})이어야 하므로 이는 모순이다. 따라서 jj는 1이어야 하고, n1=kn-1 = k이다. 정수 aa의 위수가 ϕ(n)\phi(n)을 초과하지 않으므로, n1=kϕ(n)n1n-1 = k \leq \phi(n) \leq n-1이고, ϕ(n)=n1\phi(n) = n-1이다. 즉 nn은 소수다.

(2-3) 개선

뤼카의 정리는 (n1)(n-1)를 나누는 소수 pp들 모두에 대해 성립하는 한 정수 aa를 찾아야 했다. 하지만 보다 발전된 정리에서는 (n1)(n-1)을 나누는 소수 pip_i 각각에 대해 성립하는 정수 aia_i를 찾을 수 있으면 됨을 보인다.

발전된 정리
n1n-1을 나누는 각 소수 pip_i에 대해, ain11(modn)a_i^{n-1} \equiv 1 (\mod{n})이면서 ai(n1)/pi≢1(modn)a_i^{(n-1)/p_i} \not\equiv 1 (\mod{n})aia_i가 존재하면, nn은 소수다.

n1n-1 = p1k1p2k2...prkrp_1^{k_1}p_2^{k_2}...p_r^{k_r}이라 하고(각 pip_i는 서로 다른 소수), hih_iaia_i의 법 nn에 대한 위수라 하자. 앞서 논의된 바에 따라 hin1h_i|n-1이다. hi(n1)/pih_i \mid (n-1)/p_i라 해보자. 이 경우 hiji=(n1)/pih_i j_i = (n-1)/p_i가 되는 jij_i가 존재하게 된다. 그런데 ahi1(modn)a^{h_i} \equiv 1 (\mod{n})이므로 ahiji=a(n1)/pi1(modn)a^{h_i j_i} = a^{(n-1)/p_i} \equiv 1 (\mod{n})이 된다. 이는 모순이므로 hi∤(n1)/pih_i \not\mid (n-1)/p_i여야 한다. hi(n1)h_i | (n-1)hi∤(n1)/pih_i \not\mid (n - 1)/p_i로부터 pikihip_i^{k_i} \mid h_i를 얻을 수 있디.

ii에 대해, hiϕ(n)h_i | \phi(n)이므로, 위와 함께 pikiϕ(n)p_i^{k_i} \mid \phi(n)임을 알 수 있다. 각 pip_i에 대해 pikiϕ(n)p_i^{k_i} | \phi(n)이므로 n1=p1k1p2k2...prkrϕ(n)n-1 = p_1^{k_1} p_2^{k_2} ... p_r^{k_r} | \phi(n)이다. n1ϕ(n)n-1 | \phi(n)이므로 nn은 소수가 된다.

(2-4) 문제점

문제는 위의 두 정리를 이용하려면 nn의 소수성을 검사하기 위해, n1n-1을 이루는 소인수가 무엇인지를 미리 알아야 한다는 점이다. n1n-1의 소인수를 구하는 것이나 nn의 소인수를 구하는 것이나 어렵기는 마찬가지다. 또한 각 소인수를 통해 검사를 진행해야하니 최악의 경우 너무 많은 소인수가 필요해질 수도 있다.

1914년 헨리 포클링턴은 다음의 정리를 제시했다.

포클링턴의 정리
n1=mj, m=p1k1p2k2...psks,mnn-1 = mj,\ m =p_1^{k_1}p_2^{k_2}...p_s^{k_s}, m \geq \sqrt{n}이고 gcd(m,j)=1gcd(m, j) = 1이라 하자. 각 pip_i에 대해, ain11(modn)a_i^{n-1} \equiv 1 (\mod{n})이고 gcd(ai(n1)/pi1,n)=1gcd(a_i^{(n-1)/p_i} - 1, n) = 1이면 nn은 소수다.

증명은 생략하고, 이 정리에 따르면 n1n-1의 모든 소인수를 알 필요 없이, 인수분해된 부분이 그렇지 않은 부분보다 커질 때까지만 인수분해를 진행하면 된다. 그럼에도 불구하고 문제는 그렇게 하기 위해서는 사전에 충분히 많은 n1n-1의 인수를 얻어야 한다는 점이다.

2. 밀러-라빈 소수 판별

위와 같이 페르마의 정리를 이용하면 nn의 자명하지 않은 인수를 직접 찾지 않고도 큰 홀수 정수 n>1n>1이 합성수인지의 여부를 결정할 수 있다(nn이 아닌 n1n-1의 인수를 찾아야 한다는 문제가 있기는 하지만).

밀러-라빈 소수 판별은 합성수를 직접 테스트할 수 있는 다른 방법이다. 우선 임의의 정수를 선택한다. 이 정수를 통해 nn에 대한 테스트를 진행하며, 만약에 테스트를 통과하지 못하면 nn은 확실히 합성수이고, 테스트를 통과하면 nn은 소수일 수도 있다.

임의의 양의 홀수 nn에 대해, n1=2hmn-1 = 2^hm이라 쓸 수 있다(mm은 홀수). 이제 1<a<n11 < a < n-1인 임의의 정수 aa를 선택하고, 다음의 수열을 만든다(각 항마다 modulo nmodulo\ n을 해줘야 한다).

am,a2m,a4m,...,a2h1m,a2hm=an1a^m, a^{2m},a^{4m}, ..., a^{2^{h-1}m}, a^{2^hm} = a^{n-1}

정수 nnam1(modp)a^m \equiv 1 (\mod{p}) 이거나, 어떤 j=1,2,...,h1j=1,2,...,h-1에 대해 a2jm1(modp)a^{2jm} \equiv -1 (\mod{p})를 만족하면 테스트를 통과한다. 임의의 홀수인 소수는 어떤 밑 aa에 대해서도 위의 테스트를 통과하며, 따라서 어떤 aa를 정하든, 아래의 테스트를 통과하지 못하면 그 수는 합성수라 할 수 있다.

정리 3
pp가 홀수인 소수고, p1=2hmp-1 = 2^hm이라 하자(mm은 홀수, h1h \geq 1). 이때 임의의 정수 a(1<a<p1)a(1 < a < p-1)am1(modp)a^m \equiv 1 (\mod{p}) 혹은 어떤 j=1,2,...,h1j = 1, 2,...,h-1에 대해 a2jm1(modp)a^{2^jm} \equiv -1 (\mod{p})를 만족한다.

aa가 법 pp에 대해 위수 kk를 가진다고 하자. kϕ(p)=p1=2kmk | \phi(p) = p-1 = 2^km이다. kk가 홀수이면, 유클리드의 보조 정리에 의해 kmk|m이다. 따라서 어떤 정수 rr에 대해 m=krm = kr이고 다음의 결과가 따른다.

am=(ak)r1r=1(modp)a^m = (a^k)^r \equiv 1^r = 1 (\mod{p})

이번에는 kk가 짝수라고 해보자. 이 경우 k=2j+1dk = 2^{j+1}d로 쓸 수 있다(j0j\geq 0, dd는 홀수). 2j+1d2hm2^{j+1}d | 2^hm이므로 j+1hj+1 \leq h, dmd |m이다. 또한 a2j+1d1(modp)a^{2^{j+1}d} \equiv 1 (\mod{p})로부터, a2jd±1(modp)a^{2^jd} \equiv \pm 1 (\mod{p})을 얻을 수 있다. aa의 위수가 kk이므로 a2jd1(modp)a^{2^jd} \equiv 1 (\mod{p})는 불가능하다. 따라서 a2jd1(modp)a^{2^jd} \equiv -1 (\mod{p})이다.

dmd|m이므로, 이제는 어떤 양의 홀수 tt에 대해 m=dtm = dt로 쓸 수 있다. 그 결과 다음이 성립하며, 증명이 완료된다. \square

a2jm=(a2jd)t(1)t=1(modp)a^{2^jm}= (a^{2^jd})^t \equiv (-1)^t = -1(\mod{p})

다만 어떤 정수 nn이 이 테스트를 한 번 통과했다고 해서 반드시 소수라는 보장이 주어지는 것은 아님에 주의하자. 예를 들어 n=2047=2389n = 2047 = 23 \cdot 89의 경우, n1=21023n-1 = 2 \cdot 1023이고, 210231(mod2047)2^{1023} \equiv 1 (\mod {2047})이므로 테스트를 통과한다.

밀러-라빈 소수 판별 코드는 종종 "확률론적 소수 판별법"이라 불린다. 0<ai<n0 < a_i < nkk개의 정수 a1,a2,...,aka_1, a_2, ..., a_k를 임의로 고르자. 만약 nn이 이것들을 이용한 테스트에 한 번이라도 실패한다면, nn은 확실히 합성수다. 물론 kk개의 테스트를 모두 통과했다 하더라도 nn이 합성수라는 보장은 없고, nn이 "아마도 소수"일 것이라는 강한 추측만을 할 수 있다.

합성수가 kk개의 테스트를 모두 통과할 확률은 최대 (14)k(\frac{1}{4})^k 정도임이 알려져 있기에, 현대 컴퓨터들은 100개 정도의 테스트를 둬서 nn이 소수일 확률을 최소 1(14)1001 - (\frac{1}{4}) ^{100}까지 끌어올리곤 한다.

다만 intlong long과 같이 수학적으로는 작은 수들의 경우, 특정 aa 들에 대해서만 검사를 통과하면 결정론적으로 소수를 판별할 수 있음이 알려져있다.

  • int 범위의 경우 2, 7, 61을 가지고서 테스트를 통과하면 소수이다.
  • long long 범위의 경우 37 이하의 소수를 이용해 테스트를 통과하면 소수이다.

4. 밀러-라빈 소수 판별 코드

아래는 위를 구현한 코드다. 오버플로우를 방지하며 long long끼리 곱하는 함수, 거듭제곱을 위한 함수, , 밀러-밀러 라빈 테스트 및 소수 판정을 위한 함수가 있다. long long 범위 내에서는 잘 동작하는 것 같다. 아마도.

#include <iostream>
using namespace std;

int primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41};//혹시 몰라서 41까지 테스트해보자

/**
 * 오버플로우 없이 (p * q) % mod를 구하기 위한 함수 
*/
long long mult(long long p, long long q, long long mod){
    p %= mod;//일단 나머지 연산을 한다.
    q %= mod;

    unsigned long long r = 0;
    unsigned long long w = p;

    while (q){//q의 비트를 보면서, 덧셈을 통해 곱을 만든다.
        if (q % 2) r = (r + w) % mod;//0번째 비트가 1인 경우 결과에 더해주고
        w = (2 * w) % mod;//다음으로 넣을 수도 있을 값을 계산해준다.
        q >>= 1;//다음 비트를 볼 것
    }

    return (long long)r;
}

/**
 * a^m % p를 구하는 함수. 동작 방식은 위의 mult()와 동일하다.
*/
long long pow_mod(long long a, long long m, long long p){
    long long ret = 1;
    a %= p;

    while (m){
        if (m % 2) ret = mult(ret, a, p);
        a = mult(a, a, p);
        m >>= 1;
    }

    return ret;
}

/**
 * 임의의 정수 n에 대한, 밑 a의 밀러-라빈 테스트
*/
bool millerRabin(long long n, long long a){
    long long k = n - 1;//k = 2^h * m으로 나타낼 수 있다. 이때 h >= 1, m은 홀수다.

    while (true) {//a^{2^h * m}에서부터 거꾸로 진행한다.
        long long d = pow_mod(a, k, n);//a^k가 통과하나?
        if (k % 2) return (d == 1 || d == n - 1);//만약 k가 홀수라면, a^m을 확인하고 있다는 말이다. 이 경우 a^m === 1 (mod n)인지를 확인한다.
        if (d == n - 1) return true;//k가 짝수인 경우로, a^k === -1 (mod n)인 경우다.
        k >>= 1;//2로 나눠보자
    }
}

/**
 * 임의의 정수 N에 대한 소수 판별
*/
bool isPrime(long long N){
    if (N == 1) return false;

    for (int i = 0; i < 13; i++){
        if (N == primes[i]) return true;//소수만 들어있는 배열 a의 원소와 같은 경우, 소수다.
        if (N % primes[i] == 0) return false;//소수의 배수이므로 합성수다.
        if (!millerRabin(N, primes[i])) return false;//밀러-라빈 검사를 통과하지 못한 경우 합성수다.
    }
    return true;
}

int main(){
    long long N;
    scanf("%lld", &N);

    printf("%s", isPrime(N) ? "Yes" : "No");
}

1개의 댓글

comment-user-thumbnail
2024년 4월 9일

수학자 꿈나무 박영서군 앞으로도 정진하세요

답글 달기