폴라드-로 소인수분해 알고리즘

Motivation

소인수분해는 매우 어려운 과정이다. 실제로 소인수분해는 (입력 비트 수에 대해) 다항 시간에 결정론적으로 풀 수 있는 방법이 알려져 있지 않다.

그러나, 약간의 랜덤성 (확률에 대한 의존) 을 허용한다면 상당히 clever한 알고리즘을 통해 빠르게 소인수분해를 할 수 있는데, 이 방법이 바로 Pollard’s Rho($\rho$) 알고리즘이다.

생일 문제 (Birthday Problem)

고등학교 확/통 교과서에도 실려있는 유명한 문제인데, 다음 문제에 답해 보자.

• 23명의 사람이 한 방에 모여 있다. 이 중, 적어도 한 쌍 이상이 생일이 겹칠 확률은 얼마인가?

계산을 직접 해 본다면, 별로 직관적이지 못한 결론을 얻는다. 대략 50%가 넘어가는데, 365개의 생일 중 23명이 골랐을 뿐인데 이렇게 높은 확률이라는 것이 비직관적이기 때문에 이 결과를 Birthday Paradox라고도 부른다.

이 문제의 핵심은, 마구 랜덤하게 뽑으면 생각보다 많이 겹친다 라는 정보이다. 이를 이용하여, 어떻게 소인수분해를 할 수 있는지 알아보고자 한다.

구체적으로, $\sqrt{n}$ 정도의 사람이 필요하다는 사실을 기억하자. 즉, $n$개 중 $k$개를 뽑아서, 겹치는 쌍을 만들고자 한다면, $k=O(\sqrt{n})$ 정도 뽑으면 겹치는 쌍을 기대할 수 있다. 증명은 생략.

가정

폴라드-로 알고리즘이 잘 작동하기 위해서는, 큰 수 $m$을 소인수분해하되, $m$의 가장 작은 소인수 $p$가 작아야 한다. 예를 들어, RSA에서 쓰는 좋은 N - 즉, 큰 소수 2개를 곱한 수에서는 그 효율이 현저히 낮다. 100자리 소수를 2개 곱해서 얻은 200자리 합성수보다는, 20자리 소수 10개를 곱해서 얻은 수를 잘 소인수분해하는 알고리즘임을 의미한다.

Algorithm

$k$개의 정수 $u_1,u_2,\dots u_k$ 를 랜덤하게 선택하자. 이때, 생일 문제에 의해, $u_i\equiv u_j\mathrm{mod}p$ 인 $i,j$를 얻을 가능성이 상당히 높다. 그런데, 이는 다시 말하면 $u_i-u_j$가 (일반성을 잃지 않고 양수라고 하자) $m$과 $p$를 공약수로 가질 가능성이 상당히 높다는 뜻이다. 우리는 두 수의 gcd를 입력 비트 수에 대한 다항 시간에 찾는 유클리드 알고리즘을 잘 알고 있으므로, $m$의 어떤 소인수를 빠르게 찾을 수 있다는 의미가 된다. 다시 말해, $(\genfrac{}{}{0ex}{}{k}{2})$개의 $u_i-u_j$들과 $m$의 gcd를 확인함으로서 $m$의 어떤 약수를 얻을 확률이 상당히 높다는 뜻이다.

그러나, $O(k^2)$ 개의 조합을 모두 검토하기는 너무 느리다. 폴라드-로 아이디어는 이 과정을 줄이는 것인데, $u_i$를 랜덤하게 뽑는 대신 $u_{i+1}=u_i^2+1$ 을 사용하는 것이다. 이 다항식 $f(x)=x^2+1$ 은 사실 어떤 다항식이든 크게 상관은 없으나, mod $p$에 대한 랜덤성을 해치지 않아야 하고, 계산이 너무 오래 걸리지 않아야 한다. 일반적으로 저렇게 생긴 다항식이 잘 작동함이 알려져 있으나 $x^2+x+1$ 같은걸 쓴다고 큰 문제는 없고… 다만 1차식을 쓰면 안 된다는 사실도 잘 알려져 있다.

이런식으로 결정론적으로 계산하는 것의 장점은, $u_i$ 의 수열이 주기성을 갖게 된다. 우리가 어떤 주기성을 갖는 수열 $u_i$에서 주기를 찾고자 할 때는 Tortoise and Hare 라는 좋은 방법을 쓸 수 있는데, 이 방법은 다음과 같다.

수열의 주기를 $r$이라고 할 때, $s\equiv t\mathrm{mod}r\Rightarrow u_s=u_t$ 가 된다. 따라서, $s$를 적당히 큰 (주기성을 보이는) $r$의 배수로 잡고, $t=2s$로 잡으면 항상 $u_s\equiv u_{2s}$ 이다. 따라서, 우리는 $k$를 하나씩 늘리면서 $u_{2k}-u_k$ 만 확인하여도 반드시 주기를 놓치지 않음이 보장되며, 이때 주기성을 가지는 첫번째 $r$의 배수가 대략 $O(\sqrt{p})$ 정도 스케일이기 때문에 $k=O(\sqrt{p})$ 정도 확인하면 그 안에 gcd$(u_{2k}-u_k,m)$ 가 1이 아닌 $k$를 기대할 수 있다.

예시

Niven 의 정수론 책에 수록된 예시는 다음과 같다.

$m=36,287$, $f(x)=x^2+1\mathrm{mod}m$, $u_0=1$일 때, $u$의 수열을 계산하면…
$1,2,5,26,677,22886,2439,33941,24380,3341,22173,25652,26685,29425,22806$… 이다. 이때, 각 $k$에 대해 gcd$(u_{2k}-u_k,m)$ 를 계산하면, $u_{14}-u_7$에서 gcd$(u_{14}-u_7,m)=131$ 이다. 131은 $m$의 약수이다.

구현

위 알고리즘을 그대로 코드로 옮기면 된다.

• 단, $u$ 수열을 미리 구해놓고 Tortoise-Hare를 돌리는 것은 어디까지 구해야 할지 모르는 상황에서는 별로 적절하지 않다. $x^2+1$ 은 구하기 쉬운 다항식이므로, 매 스텝마다 $y$는 두 스텝씩, $x$는 한 스텝씩 나간다고 생각하면서 진행시키자.
• $u_0$는 랜덤하게 뽑았는데, 별로 상관은 없다.
• is_prime 부분은 일반적인 소수 판정 함수를 쓰면 된다. 보통은 밀러-라빈 판정법을 많이 쓴다 (그냥 판정하면 느리니까).

inline int64_t mulmod(int64_t x, int64_t y, int64_t m)
{
    return ((__int128)x * y % m);
}
int32_t PollardRho(int32_t n)
{
    if (n==1) return n;
    if (n % 2 == 0) return 2;

    int32_t x = (rand()%(n-2))+1;
    int32_t y = x;
    int32_t c = 1;
    int32_t d = 1;
    while (1)
    {
        x = (mulmod(x, x, n) + c + n)%n;

        y = (mulmod(y, y, n) + c + n)%n;
        y = (mulmod(y, y, n) + c + n)%n;

        d = __gcd(abs(x-y), n);
        if (d == 1) continue;
        if (is_prime(d))
            return d;
        else return PollardRho(d);
    }
}

시간 복잡도

앞서 설명한 바와 같이, Birthday Paradox에 의해 실제로 구해야 하는 $u$의 값은 $O(\sqrt{p})$개이다. $n$에 대해서는 $p\le \sqrt{n}$ 임을 가정하면 (소수가 아님은 밀러라빈 등으로 확인하자) $O(n^{1/4})$ 알고리즘이라고 할 수 있다. 이는 다항식 $x^2+1$ 이 진짜 랜덤한 수들을 준다는 가정 하에 이루어진 계산인데, 실제로는 당연히 뭔가 더 복잡하고 끔찍한 분석이 필요하다. 이쪽은 잘 모르기도 하니 일단은 넘어가기로 하자.

당연히, 매우 큰 수를 다룰 때는 곱셈이나 모듈러 등이 유의미하게 오래 걸린다. 이 경우 적절하게 시간 복잡도에도 이런 항들을 곱해야 할 것이다. 우선은 64비트 (곱셈이 $O(1)$에 수행) 가정 하에서의 복잡도가 위와 같다고 알아두면 된다.