리먼 소인수분해

소개

이 글에서는 정수 $N$의 소인수분해를 결정적으로 (deterministic) $O(N^{1/3})$ 시간만에 할 수 있는 리먼 소인수분해 (Lehman Factorization) 기법을 소개하고, 유도합니다. 교과서에 소개된 구현은 다음과 같습니다.

// 가장 작은 소인수를 return
int lehmanFactorization(int n)
{
    // n > 21이어야 함
    if (n <= 21) return trialDivision(n);
    // n^(1/3)까지의 수로 나눔
    // cbrt는 세제곱근을 구하는 함수
    for (int i = 2; i < ceil(cbrt(n)); i++)
    {
        if (n % i == 0) return i;
    }
    // 이 글에서 주로 설명할 부분
    for (int k = 1; k <= ceil(cbrt(n)); k++)
    {
        double fourKN = 2 * sqrt(k * n)
        for (int a = ceil(fourKN); a <= floor(fourKN + cbrt(sqrt(n)) / 4 / sqrt(k)); a++)
        {
            if (isSquare(a * a - 4 * k * n)) return gcd(a + sqrt(b), n);
        }
    }
    // 이제 n은 소수
    return n;
}

서론

$\sqrt N$보다 작거나 같은 숫자로 정수 $N$을 나누면 $O(\sqrt N)$ 시간만에 소인수분해 할 수 있는 것은 잘 알려져 있습니다. 이 방법을 trial division이라고 부릅니다.

int trialDivision(int n)
{
    for (int i = 2; i <= sqrt(n); i++)
    {
        if (n % i == 0) return i;
    }
    return n;
}

또한, $9991$과 같은 숫자는 $9991=100^2-3^2=(100-3)\times(100+3)=97\times103$인 것을 관찰하여 쉽게 소인수분해 할 수 있습니다. 이 방법을 페르마의 소인수분해 방법이라고 부릅니다.

페르마의 소인수분해 방법으로 $N=pq$$p\ge q$를 만족하는 두 수 $p$, $q$를 찾으면, $N=\left(\frac{p+q}2\right)^2-\left(\frac{p-q}2\right)^2$가 됩니다. $\frac{p+q}2$가 정수인 것이 편리하기 때문에, $p$$q$의 홀짝이 같으면 좋습니다. 이것은 $N$$2$로 나눠지지 않을 때까지 나눠서 쉽게 달성할 수 있습니다. 그리고, 이제 $p$$q$가 홀수이므로 $\frac{p+q}2$가 최대이려면 $q=3$, $p=N/3$이어야 하며, 이 때 $\frac{p+q}2=\frac{N+9}6$입니다.

int fermatFactorization(int n)
{
    if (n % 2 == 0) return 2;
    for (int i = sqrt(n); i <= (N + 9) / 6; i++)
    {
        if (isSquare(i * i - N)) return gcd(i + sqrt(i * i - N), N);
    }
    return n;
}

이 방법은 최악의 경우 $N/6$회 정도의 loop를 돌 수 있기 때문에, 일반적인 수의 소인수분해 방법으로는 적절하지 않습니다.

두 방법의 융합

Trial division은 작은 인수를 가진 숫자를 빠르게 소인수분해 할 수 있지만, 두 인수가 서로 비슷한 숫자의 소인수분해는 느립니다. 반대로, 페르마의 소인수분해 방법은 인수가 서로 비슷한 숫자의 소인수분해가 빠르고, 작은 인수를 가진 숫자에 대해 가장 느립니다. (for 문의 순회 방향을 바꿔도 이 사실이 변하지 않습니다. 이 사실의 증명은 연습문제로 남깁니다.)

따라서, 작은 인수에 대해 나눠본 뒤 페르마의 소인수분해 방법을 적용하는 것을 생각할 수 있습니다.

int trialFermat(int n)
{
    // M은 N에 따라 바뀔 수 있는 수. 아래에 설명합니다.
    for (int i = 2; i <= M; i++)
    {
        if (n % i == 0) return i;
    }
    for (int i = sqrt(n); i <= f(M); i++)
    {
        if (isSquare(i * i - N)) return gcd(i + sqrt(i * i - N), N);
    }
    return n;
}

여기서, $f(M)$$N=pq$를 만족하는 $p$$q$에 대해 $\frac{p+q}2$의 최대값이며, $p=M$, $q=N / M$일 때를 생각하면 $f(M)=\frac{N+M^2}{2M}$입니다.

첫 번째 loop의 반복 횟수는 $M-1$, 두 번째 loop의 반복 횟수는 $\frac{N+M^2}{2M}-\sqrt N$입니다. isSquare(i * i - N) 함수가 n % i == 0보다 느리겠지만, 이 차이를 무시하면 우리의 방법은 $M - 1+\frac{N+M^2}{2M}-\sqrt N$의 시간이 걸립니다. $M$을 잘 선택하여 이 수치를 최소로 만들 수 있습니다. $$M - 1+\frac{N+M^2}{2M}-\sqrt N=\frac{3M}2+\frac{N}{2M}-\sqrt N-1\ge2\sqrt\frac{3N}{4}-\sqrt N-1=(\sqrt 3-1)\sqrt N-1$$ 부등식은 산술-기하 평균 부등식이며, $M=\sqrt\frac{N}{3}$일 때 등호가 성립합니다. $\sqrt 3-1\approx0.73205$이므로, trial division의 $\sqrt N$ 시간이 약 27% 정도 향상되었습니다. 즉, 이 방법 역시 여전히 $O(\sqrt N)$입니다.

더 큰 수로 만들기

위의 소인수분해 방법의 최악의 경우 중 하나는 $31621=107\times 307$입니다. $M=\sqrt\frac{N}{3}\approx104.64$이므로 첫 번째 loop 전체와 두 번째 loop의 거의 전체를 순회한 뒤 답을 찾을 수 있습니다. 하지만, 여기에 $3$을 곱한 $94863$은 비교적 소인수분해가 간단합니다. $94863=3\times107\times 307=307\times321=314^2-7^2$이므로, 이 수는 페르마의 소인수분해 방법의 최선의 경우에 가깝습니다.

이처럼, $2$를 곱했을 때 소인수분해가 쉬워지는 경우도 있을까요? 두 소인수의 비가 $1:2$에 가까운 $20099=101\times199$를 생각해 봅시다. $2$를 곱하면 $40198=2\times101\times 199=199\times202$인데, $199$는 짝수이고 $202$는 홀수이기 때문에 두 정수의 제곱의 차이로 나타낼 수 없습니다. 대신, $8$을 곱해주면 $160792=398\times404=401^2-3^2$가 되어서 소인수분해가 가능합니다. 짝수일 때만 $4$를 곱하는 것은 코드가 복잡하니까 이제부터 모든 경우에 $4$를 곱하도록 합시다.

하지만, 이렇게 적절한 수를 곱하는 방법을 사용하더라도, 시간 복잡도는 $O(\sqrt N)$입니다. 소인수 두 개의 차이가 $2$배를 넘지 않는 경우, 효율적으로 소인수분해 할 수 없기 때문입니다. 그런데, $6048983=2017\times2999$$24$를 곱해 보면 신기한 일이 일어납니다. $$6048983 \times 24 = 145175592 = 12049^2-53^2=(12049-53) \times(12049+53)$$ 두 소인수의 비율이 $\frac 3 2$에 가까웠기 때문에, $6$(과 기본적으로 곱해야 하는 $4$)를 곱했더니 제곱수에 매우 가까워진 겁니다.

이 현상을 시스템적으로 활용한 것이 리먼 소인수분해입니다.

따라서, 우리는 어떤 실수 $q/p$를 더 간단한 유리수 $\frac{q_0}{p_0}$으로 근사하는 것이 가능한지에 대해 알아봐야 합니다.

실수 근사하기

분모가 $N'$ 이하이면서 $0$에서 $1$ 사이인 유리수들을 크기 순서로 정렬한 것을 페리 수열 (Farey sequence)이라고 부릅니다. 페리 수열의 인접한 두 항이 기약분수 $\frac{q_1}{p_1}, \frac{q_2}{p_2}$일 때, 둘의 차이는 항상 $\frac{1}{p_1p_2}$인 것을 보일 수 있습니다.

또한, $p_1+p_2>N'$ (아니라면, $\frac{q_1+q_2}{p_1+p_2}$을 두 항 사이에 넣을 수 있으므로 페리 수열의 연속한 항이 아님)인 것도 알 수 있습니다.

$0$$1$ 사이의 임의의 실수 $q/p$에 대해, 연속된 두 페리 수열의 항 $\frac{q_1}{p_1}, \frac{q_2}{p_2}$가 존재하여 $\frac{q_1}{p_1} \le \frac{q}{p} \le \frac{q_2}{p_2}$을 만족합니다. 이제 $\frac{q}{p}$$\frac{q_1+q_2}{p_1+p_2}$보다 작으면 $\frac{q_1}{p_1}$을, 아니라면 $\frac{q_2}{p_2}$를 선택하는 근사를 생각해 보면, 이 근사의 오차는 $\frac{q_1}{p_1}$을 선택했을 때 $\frac{1}{p_1N'}$을 넘지 않고, $\frac{q_2}{p_2}$을 선택했을 때는 $\frac{1}{p_2N'}$을 넘지 않습니다.

따라서, $0$$1$ 사이의 임의의 실수 $q/p$$q_0\le p_0 \le N'$을 만족하는 $\frac{q_0}{p_0}$으로 근사할 수 있으며, 그 오차는 항상 $\frac{1}{p_0N'}$ 보다 작습니다.

정리

이제 우리는 임의의 $N'$에 대해 실수 $q/p$$|\frac{q}{p}-\frac{q_0}{p_0}|<\frac{1}{p_0N'}$을 만족하게 근사할 수 있는 것을 알았습니다. 같은 $p_0$, $q_0$를 사용하면 $|p_0q-pq_0|<\frac{p}{N'}$인 것을 알 수 있습니다. 식을 간단히 하기 위해, $N'=N''\sqrt\frac{p}{q}$로 두어 봅시다.$$|p_0q-pq_0|<\frac{p}{N'}=\frac{p}{N''\sqrt{\frac{p}{q}}}=\frac{\sqrt{pq}}{N''}=\frac{\sqrt{N}}{N''}$$ $$p_0q_0\le\frac{q_0}{p_0}p_0^2<\left(\frac{q}{p}+\frac{1}{p_0N'}\right)p_0^2<\frac{qN'^2}{p}+\frac{p_0}{N'}=\frac{qN''^2(p/q)}{p}+\frac{p_0}{N'}=N''^2+\frac{p_0}{N'}\le N''^2+1$$

따라서, 우리는 $N=pq$, $p\ge q$인 임의의 $N$$p, q$ 그리고 $N''$에 대해, 적절한 $q/p$의 근사 $q_0/p_0$이 존재하여, $p_0q_0\le N''^2+1$$|p_0q-pq_0|<\frac{\sqrt{N}}{N''}$을 만족함을 알았습니다.

또한, $(\sqrt{4Np_0q_0})^2=4pqp_0q_0\le(p_0q+pq_0)^2$이므로, $p_0q+pq_0=\sqrt{4Np_0q_0}+E$라고 놓을 수 있습니다. 이 때,$$4Np_0q_0+4E\sqrt{Np_0q_0}<(p_0q+pq_0)^2=|p_0q-pq_0|^2+4pqp_0q_0<\frac{N}{N''^2}+4Np_0q_0$$이므로 $E<\frac{\sqrt{N}}{4N''^2\sqrt{p_0q_0}}$인 것을 알 수 있습니다.

따라서, 정리하면 $N=pq$, $p\ge q$일 때 적절한 수 $k\le N''^2+1$ ($p_0q_0$)이 존재하여 $4kN$은 제곱수에 가깝습니다. 더 정확히 쓰면, $\sqrt{4kN}$이상 $\sqrt{4kN} + \frac{\sqrt{N}}{4N''^2\sqrt{k}}$미만인 어떤 정수의 제곱과 $4kN$의 차이는 또 다른 제곱수가 됩니다.

여기서 예외가 있는데, $q/p$의 근사 $q_0/p_0$가 사실 $0/1$인 경우입니다. 이 때 $p_0q_0=0$이므로 적절한 수를 곱해 $N$을 제곱수에 가까운 수로 만들 수 없게 됩니다. 이것은 $q/p<\frac{1}{N'+1}<\frac{1}{N'}=\frac{1}{N''\sqrt\frac{p}{q}}$일 때입니다. 이것은 $q/p<\frac{1}{N''^2}$과 동치인 것을 쉽게 알 수 있습니다. 따라서, $\frac{\sqrt{N}}{N''}$보다 작은 숫자들로 나눠 보고 시작하면 이 예외에 대한 처리가 끝납니다.

이 방법의 시간 복잡도는 $$\frac{\sqrt{N}}{N''}+\sum_{k=1}^{N''^2+1}\left({\frac{\sqrt{N}}{4N''^2\sqrt{k}}+1}\right)=\frac{\sqrt{N}}{N''}+N''^2+1+\frac{\sqrt{N}}{4N''^2}\sum_{k=1}^{N''^2+1}\frac{1}{\sqrt{k}}=O\left(\frac{\sqrt{N}}{N''}\right)+O(N''^2)$$입니다. $N''=N^{1/6}$으로 놓을 경우, 앞의 두 값이 같아지면서 시간 복잡도가 최소가 되는 것을 알 수 있고, 리먼 소인수분해의 시간 복잡도는 $O(N^{1/3})$이 됩니다.

지금까지 했던 것들을 코드로 정리하면 맨 앞의 코드가 됩니다. $N \le 21$인 경우, 이렇게 찾아낸 소인수가 $1$이거나 $N$일 수 있지만 $21$보다 큰 수에 대해서는 그럴 가능성이 없습니다.

최적화

$a^2-4kN$$8$로 나눈 나머지에 대해서 규칙이 있기 때문에, 절반 정도의 숫자는 계산하지 않고 건너뛸 수 있습니다. 앞에서 big-O notation으로 대충 찾고 넘어간 시간 복잡도를 더욱 분석하여 시간 복잡도를 최소로 하는 정확한 $N''$을 찾는다면 (상수 $c$에 대해 $cN^{1/6}$꼴입니다) 상수 배 더 빨라질 수 있습니다. 또한, 약수가 많은 $k$의 값부터 먼저 시도하면 인수를 더 빨리 찾을 가능성이 높습니다. 숫자가 64비트 정수 범위를 넘어선다면, $O(N^{1/4})$ 이하의 더 나은 시간 복잡도를 갖는 방법으로 바꾸는 것이 좋습니다.

Reference

Crandall, R. E., & Pomerance, C. (2005). Prime numbers: a computational perspective. 2nd ed. New York, NY: Springer.

내용의 대부분을 참고했습니다.

Hardy, G. H. 1., & Wright, E. M. 1. (1968). An introduction to the theory of numbers. 4th ed. [reprinted (with corrections)] Oxford: Clarendon Press.

페리 수열로 실수를 근사하는 부분을 참고했습니다.


댓글 댓글 쓰기