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

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

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

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

(1) 페르마의 소정리

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

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

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

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

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

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

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

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

Corollary

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

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

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

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

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

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

(2-1) 위수

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

정수 $a$가 법 $n$에 대해 위수 $k$를 가진다고 하자. 이때, 다음이 성립한다.

정리 1
$a^h \equiv 1 (\mod{n}) \iff k|h$

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

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

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

다음의 정리도 염두에 두자

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

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

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

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

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

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

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

(2-3) 개선

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

발전된 정리
$n-1$을 나누는 각 소수 $p_i$에 대해, $a_i^{n-1} \equiv 1 (\mod{n})$이면서 $a_i^{(n-1)/p_i} \not\equiv 1 (\mod{n})$인 $a_i$가 존재하면, $n$은 소수다.

$n-1$ = $p_1^{k_1}p_2^{k_2}...p_r^{k_r}$이라 하고(각 $p_i$는 서로 다른 소수), $h_i$를 $a_i$의 법 $n$에 대한 위수라 하자. 앞서 논의된 바에 따라 $h_i|n-1$이다. $h_i \mid (n-1)/p_i$라 해보자. 이 경우 $h_i j_i = (n-1)/p_i$가 되는 $j_i$가 존재하게 된다. 그런데 $a^{h_i} \equiv 1 (\mod{n})$이므로 $a^{h_i j_i} = a^{(n-1)/p_i} \equiv 1 (\mod{n})$이 된다. 이는 모순이므로 $h_i \not\mid (n-1)/p_i$여야 한다. $h_i | (n-1)$와 $h_i \not\mid (n - 1)/p_i$로부터 $p_i^{k_i} \mid h_i$를 얻을 수 있디.

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

(2-4) 문제점

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

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

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

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

2. 밀러-라빈 소수 판별

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

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

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

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

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

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

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

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

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

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

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

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

• 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");
}