[C++ 알고리즘] 밀러-라빈 소수 판별법 (Miller-Rabin Primality Test)

이 포스트에서는 알고리즘의 일종인 밀러-라빈 소수 판별법의 원리와 예제 풀이에 대해 다룹니다.

 

밀러-라빈 소수 판별법은 어떤 자연수 N이 소수인지를 확률적으로 판단하는 알고리즘입니다.

몇 가지 경우에 대해서만 검사를 거치므로 다른 소수 판별 알고리즘에 비해 훨씬 빠르게 작동합니다.

다만 int나 unsigned long long과 같이 큰 범위 내에서도 특정 값에 대해서만 검사를 하면 충분히 모두 검사가 된다는 증명이 되었기 때문에, 우리는 그 값을 가져다가 사용하면 주어진 범위 내에서는 완전한 소수 판별이 가능합니다.

 

 

 

** 정정 : unsigned long long에 대해서 a = 37 이하의 소수(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)에 대해서 검사하는 것 역시 충분하며, 이 방법이 외우기도 더 편합니다.

 

밀러-라빈 소수 판별법에 대한 대략적인 증명과 식은 위와 같이 정리를 해놓았습니다. (다만 자세한 부분까지는 증명을 모두 적으면 내용이 아주 길어지므로, 필요에 따라서는 식만 정리된 부분도 있습니다.)

전반적으로는 페르마의 소정리(Lemma) 2개를 활용하여 증명되어 있습니다.

 

어디까지나 알고리즘이기 때문에 코드로 직접 구현하는 것이 중요하고, 따라서 이에 맞는 예제를 풀이해보도록 하겠습니다.

풀이할 문제는 프로그래밍 문제 사이트 백준 Online Judge(BOJ)의 5615번 : '아파트 임대'이고, 문제 난이도는 Solved.ac 기준 Platinum I에 해당합니다.

 

 

5651번 : 아파트 임대

5615번: 아파트 임대

 

단순히 보았을 때에는 소수 판별과 관련이 없어 보이지만, 대응되는 x, y가 1개만 되도록 식을 유지하면서 면적에 해당하는 식 2xy+x+y를 적절히 변형하여 소인수분해를 해주면 해당 식이 소수인지를 통해 x, y가 존재하는지를 확인할 수 있습니다.

 

 

 

따라서 식을 위와 같이 변형하고 2A+1에 대해 밀러-라빈 소수판별법을 사용하여 소수라면 있을 수 없는 면적에 해당합니다.

 

 

#include <bits/stdc++.h>
#define ull unsigned long long
using namespace std;

ull Base[3] = {2, 7, 61};

ull Power(ull x, ull y, ull mod) { // ret = (x^y)%mod
    ull ret = 1;
    x %= mod;
    while(y) {
        if(y%2 == 1) ret = (ret*x)%mod;
        y /= 2;
        x = (x*x)%mod;
    }
    return ret;
}

bool isPrime(ull n, ull a) {
    if(a%n == 0) return true;
    ull k = n-1;
    while(1) {
        ull temp = Power(a, k, n);
        if(temp == n-1) return true;
        if(k%2 == 1) return (temp == 1 || temp == n-1);
        k /= 2;
    }
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL); cout.tie(NULL);

    int N, compositeN = 0;
    cin >> N;
    for(int i=0; i<N; i++) {
        ull A;
        cin >> A;
        for(int j=0; j<3; j++)
            if(!isPrime(A*2+1, Base[j])) {
                compositeN++;
                break;
            }
    }
    cout << N - compositeN;
}

 

isPrime 함수가 바로 밀러-라빈 알고리즘이 사용된 소수판별 함수입니다.

각 Base(= a)에 대해 위에 첨부한 정리가 만족하는지를 확인하고, 하나라도 만족하지 않는다면 합성수에 해당합니다.

 

k가 홀수가 될 때까지 계속 2로 나누면 결국 k = d×2^r로 분해하는 꼴이 되고, 이 때 모든 r에 대해 mod n이 -1 즉 n-1이 되는지를 확인합니다.

마지막에 k가 홀수가 되었을 때 k = d가 되며 이 때 mod N 값이 1인지를 검사합니다.

쉽게 말해서 그냥 위의 정리에 있는 식을 그대로 코드로 구현한 것입니다.

 

다만 거듭제곱을 수행하는 데 있어서 시간복잡도를 O(log k)로 줄이기 위해 빠른 거듭제곱 알고리즘을 사용해주어야 합니다.

간략히만 서술하자면, 지수 k를 이진수로 변환하듯이 하여 2의 거듭제곱꼴로 제곱을 수행할 경우에는 제곱한 값을 다시 제곱하는 원리를 이용합니다.

그리고 거듭제곱을 수행하는 과정에서 unsigned long long 범위조차도 벗어나 overflow가 발생할 여지가 있고 어차피 mod 값을 계산할 것이므로 mod 값으로 modulo 연산을 가능한 부분마다 계속 수행해줍니다.

 

 

 

풀이를 제출하면 위와 같이 정답 처리를 받을 수 있습니다.

 

---

+ 2^32 ~ 2^64 - 1인 큰 N에 대해서

 

해결책 1 : 곱셈, 덧셈에 대한 오버플로우 처리

소수인지 판별하려는 N의 크기가 2^32 이상인 경우에는 Power 함수에서 두 수를 곱하거나, 심지어는 더하기만 해도 오버플로우가 발생할 수 있습니다.

따라서 이 경우 오버플로우가 발생하지 않도록 modulo 연산을 더욱 추가하여 함수를 재조정하여 작성해 줄 필요가 있습니다.

 

#include <bits/stdc++.h>
#define ull unsigned long long
using namespace std;

ull Base[] = {2,325, 9375, 28178, 450775, 9780504, 1795265022};

ull Add(ull x, ull y, ull mod) { // ret = (x+y)%mod
    x %= mod, y %= mod;
    if(x+y >= mod) return x+y-mod;
    else return x+y;
}

ull Mul(ull x, ull y, ull mod) { // ret = (x*y)%mod
    x %= mod, y %= mod;
    ull ret = 0;
    while(y > 0) {
        if(y%2 == 1) ret = Add(ret, x, mod);
        x = Add(x, x, mod);
        y /= 2;
    }
    return ret;
}

ull Power(ull x, ull y, ull mod) { // ret = (x^y)%mod
    x %= mod;
    ull ret = 1;
    while(y > 0) {
        if(y%2 == 1) ret = Mul(ret, x, mod);
        x = Mul(x, x, mod);
        y /= 2;
    }
    return ret;
}

bool isPrime(ull n, ull a) {
    if(a%n == 0) return true;
    ull k = n-1;
    while(1) {
        ull temp = Power(a, k, n);
        if(temp == n-1) return true;
        if(k%2 == 1) return (temp == 1 || temp == n-1);
        k /= 2;
    }
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL); cout.tie(NULL);

    int N, compositeN = 0;
    cin >> N;
    for(int i=0; i<N; i++) {
        ull A;
        cin >> A;
        for(int j=0; j<4; j++)
            if(!isPrime(A*2+1, Base[j])) {
                compositeN++;
                break;
            }
    }
    cout << N - compositeN;
}

 

그런데 이 코드는 2^64 미만의 범위를 다루는 문제임과 동시에 시간 제한이 아주 넉넉할 때만 통과가 됩니다.

코드를 몇 번 고쳐보면서 시간을 재봤는데 Mul 함수를 호출하는 과정에서 시간이 많이 소요됩니다.

그리고 위의 코드의 경우에는 추가로 작은 n에 대해서는 밀러-라빈 알고리즘을 돌리는 것보다 그냥 이중 for문으로 소수를 판정하는 것이 더 빠릅니다.

 

 

해결책 2 : 128bit 자료형 사용

#include <bits/stdc++.h>
#define int128 __int128
using namespace std;

int128 Base[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};

int128 Power(int128 x, int128 y, int128 mod) { // ret = (x^y)%mod
    x %= mod;
    int128 ret = 1;
    while(y > 0) {
        if(y%2 == 1) ret = (ret*x)%mod;
        x = (x*x)%mod;
        y /= 2;
    }
    return ret;
}

bool isPrime(int128 n, int128 a) {
    if(a%n == 0) return true;
    int128 k = n-1;
    while(1) {
        int128 temp = Power(a, k, n);
        if(temp == n-1) return true;
        if(k%2 == 1) return (temp == 1 || temp == n-1);
        k /= 2;
    }
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL); cout.tie(NULL);

    int N, compositeN = 0;
    cin >> N;
    for(int i=0; i<N; i++) {
        unsigned long long input;
        cin >> input;
        int128 A = input;
        for(int j=0; j<4; j++)
            if(!isPrime(A*2+1, Base[j])) {
                compositeN++;
                break;
            }
    }
    cout << N - compositeN;
}

 

그냥 위와 같이 128bit의 정수를 지원하는 경우에는 128bit 자료형을 사용하여 풀이하면 됩니다.

 

+ 최소한 위의 문제에서는, 소수를 위의 7개를 모두 검사 대상으로 사용하지 말고, 적당히 3개나 4개 정도만 잡아서 검사해도 반례가 나오지 않습니다.

왜냐하면 unsigned long long 범위 내에서 반례가 없도록 잡은 소수들이기 때문에, int 범위에서는 적당히만 잡아도 웬만한 테스트케이스들에 대해서는 반례없이 통과 됩니다.