[C++ 알고리즘] 폴라드 로 알고리즘 (Pollard's rho algorithm, 큰 수 소인수분해 알고리즘)

이 포스트에서는 가장 빠른 인수분해 알고리즘인 폴라드 로 알고리즘에 대해 다룹니다.

 

일반적인 소인수분해 알고리즘의 경우에는 O(n^(1/2))이 최선이지만, 폴라드 로 알고리즘으로 소인수분해를 수행할 경우 평균적으로 O(n^(1/4)) 시간에 소인수분해를 수행할 수 있습니다.

 

대신 아마도 소인수분해가 성공적으로 이루어졌을 것이라고 추측이 가능한 것뿐이지, 모든 경우에 대해 완벽히 소인수분해를 수행했다고 하기 위해서는 특정 범위 내에서의 (확인을 통한) 증명이 필요합니다.

 

다만 unsigned long long과 같은 C++에서 사용할 수 있는 가장 큰 범위의 자연수 내에서도 일반적으로는 오류없이 수행이 되기 때문에, 충분히 알고리즘 문제에서 활용이 가능한 알고리즘입니다.

 

 

 

폴라드 로 알고리즘의 아이디어는, f(x) = (x^2 + c) mod n 함수에 대해 f(x)를 반복적으로 수행하다보면 사이클이 형성되는데 이 때 사이클을 이동하며 찾은 두 수 차와 n 사이의 GCD가 1이 아니라면 그 GCD가 n의 약수가 된다는 것입니다.

위의 내용은 이를 수식적으로 정리해둔 것으로, mod p에 대해 같은 modulo 값을 가지는 x_s와 x_t가 존재할 경우 p는 n의 약수임이 증명되어 있고 이를 활용하여 코드로 구현하는 것입니다.

 

그리고 가장 중요한 점이 폴라드 로 알고리즘을 돌리기 위해서는 소수가 아닌 n에 대해서만 돌려야합니다.

즉, 소수인지를 먼저 판별하고 아닐 경우에만 폴라드 로 알고리즘을 사용해야 시간적 손해와 정확도가 최적화됩니다.

따라서 가장 빠르게 소수를 판별할 수 있는 알고리즘인 밀러-라빈 알고리즘을 먼저 사용해주어야 합니다.

 

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

↑ 밀러-라빈 알고리즘에 대한 설명은 위의 링크에 정리되어 있습니다.

 

위의 알고리즘을 코드로 구현하기 위해, 적절한 예제를 풀이해보도록 하겠습니다.

풀이할 예제는 프로그래밍 문제 사이트 백준 Online Judge(BOJ)의 4149번 : '큰 수 소인수분해'이며, 문제 난이도는 Solved.ac 기준 Diamond V에 해당합니다.

 

 

4149번 : 큰 수 소인수분해

4149번: 큰 수 소인수분해

 

아주 간단한 문제인데도 불구하고 Diamond V의 난이도를 가진 문제입니다.

소인수분해를 수행할 수의 입력 범위가 무려 2^62 미만이기 때문에, unsigned long long으로 입력을 받아야만 합니다.

 

 

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

ull 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 (ull)ret;
}

bool checkPrime(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;
    }
}

bool isPrime(ull n) {
    if(n == 2 || n == 3) return true;
    if(n%2 == 0) return false;

    ull a[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    for(int i=0; i<12; i++)
        if(!checkPrime(n, a[i])) {
            return false;
            break;
        }
    return true;
}

ull GCD(ull a, ull b) {
    if(a < b) swap(a, b);
    while(b != 0) {
        ull r = a%b;
        a = b;
        b = r;
    }
    return a;
}

ull findDiv(__int128 n) {
    if(n%2 == 0) return 2;
    if(isPrime(n)) return n;

    __int128 x = rand()%(n-2) + 2, y = x, c = rand()%10 + 1, g = 1;
    while(g == 1) {
        x = (x*x%n + c)%n;
        y = (y*y%n + c)%n;
        y = (y*y%n + c)%n;
        g = GCD(abs(x-y), n);
        if(g == n) return findDiv(n);
    }
    if(isPrime(g)) return g;
    else return findDiv(g);
}

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

    ull N;
    cin >> N;
    vector<ull> ans;
    while(N > 1) {
        ull div = findDiv(N);
        ans.push_back(div);
        N /= div;
    }
    sort(ans.begin(), ans.end());
    for(int i=0; i<ans.size(); i++) cout << ans[i] << "\n";
}

 

폴라드 로 알고리즘의 구현 코드는 위와 같습니다.

 

가장 먼저 주의할 점은 입력 범위가 2^62 미만이고 unsigned long long의 범위는 2^64 미만이기 때문에, 곱셈을 수행하는 과정에서 unsigned long long 범위에서조차도 overflow가 발생하게 됩니다.

따라서 비표준 자료형인 __int128을 사용해주거나, 또는 overflow 방지용 modulo 계산이나 multiply, addition을 따로 구현해주어야 합니다.

다만 백준 Online Judge에서는 __int128 자료형을 지원하기 때문에, 이를 활용하여 문제를 풀이하는 것을 권장합니다.

 

f(x) 식에서의 c 값은 1~10 사이의 정수로 잡아주면 적당하다고 합니다.

x와 y 모두 2~n 사이의 랜덤한 수를 동일하게 잡아주고 x는 1칸씩, y는 2칸씩 이동하며 적절한 x_s와 x_t를 찾도록 구현합니다.

g(x_s - x_t와 n의 최대공약수)가 1이 아닌 수가 나올 때까지 사이클을 순회하고, g 값이 n으로 구해진다면 다른 랜덤 값에 대해 다시 검사를 수행합니다.

1과 n이 아닌 다른 g값이 얻어진다면 이는 n의 약수가 됩니다.

이 g가 소수인지 아닌지를 판별하여 소수라면 벡터에 push_back 해주고, 소수가 아니라면 폴라드 로 알고리즘을 재귀적으로 돌려 다시 소인수분해를 해줍니다.

 

 

 

풀이 코드를 제출하면 위와 같은 결과를 얻을 수 있습니다.

overflow만 주의하면 모든 테스트케이스에 대해 통과할 수 있으니 만약 일부 틀린 케이스가 존재한다면 변수들의 자료형에만 집중해서 검토해보면 좋을 것입니다.