밀러-라빈 소수 판별법

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

알고리즘의 아이디어는 페르마의 소정리에서 출발합니다. $n$이 홀수라고 했으니 $n-1$은 짝수이고, 따라서 적당한 홀수 $d$와 정수 $s$에 대해 $n-1=d\times2^s$입니다. 만약 $n$이 정말로 소수라면 $n$보다 작은 임의의 양의 정수 $a$에 대해

\[ a^{n-1} = a^{d\times2^s} \equiv 1 \pmod n \]

이고, 합동식의 정의에 따라

\[ a^{d\times2^s} - 1 = \left(a^{d\times2^{s-1}}+1\right)\left(a^{d\times2^{s-2}}+1\right) \cdots\left(a^d+1\right)\left(a^d-1\right) \]

은 $n$의 배수입니다. $n$이 소수라 가정했으므로 우변의 인수 중 적어도 하나는 $n$의 배수여야 합니다. 따라서, 적당한 정수 $a<n$을 골랐을 때 인수 $s+1$개 중 $n$의 배수가 존재한다면 $n$은 아마도 소수일 것이라 유추할 수 있습니다. 합성수를 소수로 잘못 판단할 확률을 줄이고 싶다면, 최대한 많은 $a$에 대해 시험해보면 됩니다.

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

\begin{align*} 174^{55\times2^1}+1 &= 174^{110}+1 \equiv 0 \pmod{221} \\ 174^{55}+1 &= 174^{55}+1 \equiv 48 \not\equiv 0 \pmod{221} \\ 174^{55}-1 &= 174^{55}-1 \equiv 46 \not\equiv 0 \pmod{221} \end{align*}

인수 중 하나가 221의 배수이므로 221은 아마도 소수입니다. 두 번째로 $a=137$을 선택해봅시다.

\begin{align*} 137^{55\times2^1}+1 &= 137^{110}+1 \equiv 206 \not\equiv 0 \pmod{221} \\ 137^{55}+1 &= 137^{55}+1 \equiv 189 \not\equiv 0 \pmod{221} \\ 137^{55}-1 &= 137^{55}-1 \equiv 187 \not\equiv 0 \pmod{221} \end{align*}

인수 중 221의 배수가 하나도 없으므로 221은 확실히 합성수입니다. 실제로 221=13×17입니다.

다행히도, $n$이 작으면 $a$를 이 정도만 넣어봐도 충분하다고 증명한 결과가 있습니다. $n$이 2³²보다 작으면 $a$에 2, 7, 61만 넣어봐도 잘못 판별하는 일이 없고, 2⁶⁴보다 작으면 2, 325, 9375, 28178, 450775, 9780504, 1795265022만 넣어봐도 충분합니다.

다음은 각각 $n$이 2³² 미만일 때와 2⁶⁴ 미만일 때 소수를 확실히 판별하는 코드입니다. 양쪽 다 $n$이 매우 작을 때에는 무식하게 모든 수로 나눠보는 방법이 더 빠르기 때문에 따로 처리합니다. ($a$가 $n$보다 작아야 한다는 조건도 있고요.)


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;
}


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보다 큰 자료형이 C에 없기 때문에 addmod, mulmod 함수를 써서 안전하게 합과 곱을 구해줍니다.

응용

BOJ 5615:: 아파트 임대

넓이를 $k$라 하면 $k=2xy+x+y$이고

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

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

Algorithm