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

밀러-라빈 소수 판별법은 어떤 홀수 $n$이 소수인지 확률적으로 판별해주는 알고리즘입니다. 여기서 확률적이라는 것은, 이 알고리즘은 주어진 $n$이 “합성수이다” 또는 “아마도 소수일 것이다”는 대답을 내놓는다는 뜻입니다. 그러므로 이 알고리즘은 합성수를 소수라고 잘못 판별할 수 있습니다.

$n$이 홀수라고 했으니 $n-1$은 짝수이고, 따라서 적당한 홀수 $d$와 정수 $s$가 존재하여 $n-1 = d \cdot 2^s$입니다. 만약 $n$이 소수라면, $n$보다 작은 양의 정수 $a$에 대해 다음 중 하나가 성립합니다.

  1. $a^d \equiv 1 \pmod n$
  2. $r=0,1,2,\cdots,s-1$중 적어도 하나에 대해 $a^{d \cdot 2^r} \equiv -1 \pmod n$

증명) 먼저 보조정리 하나를 증명하겠습니다.

보조정리) 소수 $p$에 대해 $x^2 \equiv 1 \pmod p$이면 $x \equiv 1 \pmod p$이거나 $x \equiv -1 \pmod p$이다.
증명) 합동식의 정의에서 $x^2-1=(x+1)(x-1)$은 $p$의 배수이고, 따라서 $x+1$과 $x-1$ 둘 중 하나는 $p$의 배수여야 합니다.

페르마의 소정리에 의해서 $a^{n-1} = a^{d \cdot 2^r} \equiv 1 \pmod n$입니다. 보조정리에 의해 $a^{d \cdot 2^{r-1}} \equiv 1 \pmod n$이거나 $a^{d \cdot 2^{r-1}} \equiv -1 \pmod n$이 되겠죠. 후자가 성립하면 2가 성립하므로 증명 끝입니다. 전자가 성립할 경우 다시 보조정리에 의해 $a^{d \cdot 2^{r-2}} \equiv 1 \pmod n$이거나 $a^{d \cdot 2^{r-2}} \equiv -1 \pmod n$입니다. 또 전자가 성립하면 계속 보조정리를 써서 $a^{d}$까지 갑니다. 끝까지 갔는데도 $a^d \equiv 1 \pmod n$이면 이번엔 1번이 성립합니다. 결국 둘 중 하나는 성립하게 되어 있습니다.

역으로 어떤 홀수 $n$이 $n$보다 작은 양의 정수 $a$에 대해 다음을 만족하면 $n$은 무조건 합성수입니다.

  1. $a^d \not\equiv 1 \pmod n$
  2. $r=0,1,2 \cdots,s-1$ 모두에 대해 $a^{d \cdot 2^r} \not\equiv -1 \pmod n$

따라서 적당히 $a$를 하나 고르고, 위 조건이 성립하면 $n$이 합성수라고 답하고, 성립하지 않으면 $n$은 아마도 소수(probable prime)라고 답합니다. 합성수를 소수라고 잘못 판별할 확률을 줄이고 싶다면, 최대한 많이 $a$의 값을 바꿔가며 시험해보면 됩니다. 아, 물론 $a=1$이면 1번 조건이 항상 거짓이므로 아무런 도움이 안 됩니다…

$n=221$로 예를 들어보겠습니다. $n-1 = 55 \cdot 2^2$이므로 $d=55$, $s=2$입니다. 첫 번째로 $a=174$를 선택해봅시다.

  1. $a^d = 174^{55} \equiv 47 \not\equiv 1 \pmod {221}$
  2. $a^{d \cdot 2^0} = 174^{55} \equiv 47 \not\equiv -1 \pmod {221}$,
    $a^{d \cdot 2^1} = 174^{110} \equiv -1 \pmod {221}$

2번이 성립하지 않으므로 221은 아마도 소수입니다. 두 번째로 $a=137$을 선택해봅시다.

  1. $a^d = 137^{55} \equiv 188 \not\equiv 1 \pmod {221}$
  2. $a^{d \cdot 2^0} = 137^{55} \equiv 188 \not\equiv -1 \pmod {221}$,
    $a^{d \cdot 2^1} = 137^{110} \equiv 205 \not\equiv -1 \pmod {221}$

1번과 2번 모두 성립하므로 221은 소수가 아닙니다.

결국 $a$를 많이 넣어보지 않으면 높은 확률로 틀린다는 건데, $n$이 작으면 $a$를 이 정도만 넣어봐도 충분하다고 계산해놓은 게 있습니다. $n$이 $2^{32}$보다 작은 합성수이면(unsigned int에 저장 가능한 정수라면) $a$에 2, 7, 61만 넣어봐도 소수로 잘못 판별되는 일이 없다고 합니다. 또, $n$이 $2^{64}$보다 작은 합성수이면(unsigned long long) $a$에 2, 325, 9375, 28178, 450775, 9780504, 1795265022만 넣어봐도 충분합니다.

$2^{32}$ 미만

using ull = unsigned long long;

vector<ull> alist = {2, 7, 61};

// calculate x^y % m
ull powmod(ull x, ull y, ull m) {
    x %= m;
    ull r = 1ULL;
    while (y > 0) {
        if (y % 2 == 1)
            r = (r * x) % m;
        x = (x * x) % m;
        y /= 2;
    }
    return r;
}

// true for probable prime, false for composite
inline bool miller_rabin(ull n, ull a) {
    ull d = n - 1;
    while (d % 2 == 0) {
        if (powmod(a, d, n) == n-1)
            return true;
        d /= 2;
    }
    ull tmp = powmod(a, d, n);
    return tmp == n-1 || tmp == 1;
}

bool is_prime(ull n) {
    if (n <= 1)
        return false;
    if (n <= 10000ULL) {
        for (ull i = 2; i*i <= n; i++)
            if (n % i == 0)
                return false;
        return true;
    }
    for (ull a : alist)
        if (!miller_rabin(n, a))
            return false;
    return true;
}

unsigned long long을 쓰는 이유는 unsigned int 두 개의 곱을 구할 때 오버플로를 방지하기 위해서입니다.

$2^{64}$ 미만

using ull = unsigned long long;

vector<ull> alist = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};

// calculate (x + y) % m; overflow-safe
inline ull addmod(ull x, ull y, ull m) {
    x %= m;
    y %= m;
    return (x >= m - y? x - (m - y) : x + y);
}

// calculate (x * y) % m; overlow-safe
inline ull mulmod(ull x, ull y, ull m) {
    x %= m;
    y %= m;
    ull r = 0;
    while (y > 0) {
        if (y % 2 == 1)
            r = addmod(r, x, m);
        x = addmod(x, x, m);
        y /= 2;
    }
    return r;
}

// calculate x^y % m; overflow-safe
ull powmod(ull x, ull y, ull m) {
    x %= m;
    ull r = 1ULL;
    while (y > 0) {
        if (y % 2 == 1)
            r = mulmod(r, x, m);
        x = mulmod(x, x, m);
        y /= 2;
    }
    return r;
}

// true for probable prime, false for composite
inline bool miller_rabin(ull n, ull a) {
    ull d = n - 1;
    while (d % 2 == 0) {
        if (powmod(a, d, n) == n-1)
            return true;
        d /= 2;
    }
    ull tmp = powmod(a, d, n);
    return tmp == n-1 || tmp == 1;
}

bool is_prime(ull n) {
    if (n <= 1)
        return false;
    if (n <= 10000000000ULL) {
        for (ull i = 2; i*i <= n; i++)
            if (n % i == 0)
                return false;
        return true;
    }
    for (ull a : alist)
        if (!miller_rabin(n, a))
            return false;
    return true;
}

unsigned long long보다 큰 자료형이 없기 때문에 addmod, mulmod 함수를 써서 오버플로가 일어나지 않게 합과 곱을 구해줍니다. 덕분에 코드가 조금 길어졌습니다. 만약 비표준 확장인 __int128을 쓸 수 있는 환경이라면 쓰는 편이 좋습니다. mulmod 함수는 오버플로에서 안전한 대신 많이 느리거든요. (시간 복잡도 $O(\log y)$)

두 경우 모두 $n$이 작으면 무식하게 모든 수로 나눠보는 방식이 더 빠르기 때문에 따로 처리합니다.

시간 복잡도

powmod 함수가 mulmod를 쓰지 않으면 시간 복잡도가 $O(\log y)$입니다. miller_rabin함수는  while문이 $O(\log n)$번 돌고 그 안쪽에서 powmod 함수를 호출하므로 시간 복잡도는 $O(\log^2 n)$이 됩니다. 한편 mulmod를 쓰면 powmod의 시간 복잡도가 $O(\log^2 y)$이므로 miller_rabin의 시간 복잡도는 $O(\log^3 n)$입니다. $a$의 개수가 $k$개이면 밀러-라빈 소수 판정법의 총 시간 복잡도는 $O(k\log^2 n)$ 또는  $O(k\log^3 n)$이 되어 평균적으로 직접 나눠보는 $O(\sqrt n)$ 알고리즘보다 훨씬 빠릅니다.

응용

BOJ 5615:: 아파트 임대

https://www.acmicpc.net/problem/5615

넓이를 $k$라고 하면

\[ 2xy + x + y = k \]

\[ 4xy + 2x + 2y + 1 = (2x+1)(2y+1) = 2k + 1 \]

따라서 $2k+1$이 소수인지 판별하면 됩니다.