카라추바 알고리즘 - 알고리즘카

1960년, 안드레이 콜모고로프는 소련 정보과학의 다른 선구자들과 함께 한 학술 세미나에 모여 “$n^2$ 가설"을 제시했습니다. 두 $n$자리 수를 $O(n^2)$보다 빠르게 곱하는 것은 불가능하다는 것입니다. 이는 적어도 4천 년 전 수메르인들이 고안했고 당시까지 누구도 능가하지 못한 “세로셈” 곱셈이 두 수를 곱하는 점근적으로 최적인 알고리즘이라는 것을 의미합니다.

일주일 후, 23세의 대학원생 아나톨리 카라추바는 실행 시간이 $O(n^{\log_2 3})$인 곱셈 방법을 제안하여 이 가설을 반박하였습니다. 이 글에서는 바로 이 방법을 살펴보겠습니다.

#알고리즘

알고리즘의 기본 아이디어는 매우 간단합니다. 길이가 $n$인 두 수의 곱셈을 약간의 대수적 트릭을 이용해 길이가 $\frac{n}{2}$인 수의 곱셈 세 번으로 환원합니다. 재귀의 각 단계마다 원소의 수는 늘어나지만, 가장 아래 단계에서는 원소가 모두 합해 $O(n^{\log_2 3})$개뿐이며, 이것이 바로 이처럼 특이한 점근적 복잡도가 나오는 이유입니다.

우리가 다룰 알고리즘에서는 수 대신 다항식을 곱합니다. 이는 동등한 문제입니다. $x$ 대신 진법의 밑을 대입하면, 수의 각 자릿수를 차례로 계수로 삼을 수 있습니다.

$$ \begin{aligned} A(x) &= a_0 + a_1\cdot x + a_2 \cdot x^2 + \dots + a_n \cdot x^n \\ &= a_0 + a_1\cdot 10 + a_2 \cdot 10^2 + \dots + a_n \cdot 10^n \end{aligned} $$ 그러면 길이가 $n = 2k$로 같은 두 다항식 $a(x)$와 $b(x)$가 있고, 이 둘을 곱하려 한다고 합시다. 이들의 계수를 똑같은 크기의 두 부분으로 나누어 다음과 같이 나타냅시다. $$ a(x) = a_1 (x) + x^k a_2(x) \\b(x) = b_1 (x) + x^k b_2(x) $$ 이제 곱 다항식 $p_1$과 $p_2$를 재귀적으로 계산합시다: $$ p_1(x) = a_1(x) \cdot b_1(x) \\ p_2(x) = a_2(x) \cdot b_2(x) $$ 그리고 다항식 $t$도 계산합니다: $$ t(x) = ( a_1(x) + a_2(x) ) \cdot (b_1(x) + b_2(x)) $$ 원래 다항식들의 곱셈 결과(크기 $2n$의 다항식)는 이제 다음 공식으로 계산할 수 있습니다. 주의하십시오, 대수학이 나옵니다: $$ c(x) = a(x) \cdot b(x) = p_1(x) + x^k \cdot (t(x) - p_1(x) - p_2(x)) + x^{2k} \cdot p_2(x) $$

필요한 대입을 수행해 보기만 하면 공식이 올바른지 확인할 수 있습니다.

#분석

필요한 연산의 수를 세어 보면, 크기가 $n$인 두 다항식을 곱하기 위해서는 크기가 $\frac{n}{2}$인 세 개의 곱 $p_1$, $p_2$, $t$를 계산하고, 상수 번의 덧셈, 뺄셈, 시프트($x^k$를 곱하는 연산)를 수행해야 하며, 이들은 모두 합쳐 $O(n)$에 수행할 수 있다는 것을 알 수 있습니다.

마스터 정리에 따르면, 이 경우 전체 알고리즘의 점근적 복잡도는 $\Theta (n^{\log_2 3}) \approx \Theta (n^{1.58})$입니다. 문제가 크기가 $b = 2$배 작은 $a = 3$개의 부분으로 나뉘고, 결합은 $O(n)$에 이루어지기 때문입니다.

#구현

간단히 하기 위해 $n$이 2의 거듭제곱이라고 가정하겠습니다. 그렇지 않은 경우에는 상황에 따라 다음 두 가지 임시방편 중 하나로 이를 해결할 수 있습니다:

다항식의 계수를 가장 가까운 2의 거듭제곱까지 0으로 채울 수 있습니다. 최악의 경우 이 방법은 $2^{1.58} \approx 3$배 더 오래 걸립니다.
다항식에서 마지막 계수를 «떼어 내어» 크기 $(2k + 1)$의 문제를 크기 $2k$의 문제와 상수 번의 덧셈으로 환원할 수 있습니다.

구현하기가 더 간단하므로 첫 번째 방법을 사용하겠습니다.

효율적인 구현을 위한 주요 고려 사항은 다음과 같습니다.

불필요한 메모리 할당을 최대한 줄여야 하며, 이를 위해 기존 배열을 재사용해야 합니다.
모든 산술 연산은 컴파일러가 벡터화할 수 있도록 배열에 대한 단순한 선형 순회로 구현해야 합니다.
if (n == 1) c[0] = a[0] * b[0] 형태의 기저 사례를 사용하는 대신, 문제 크기가 어느 정도 이하가 되면 더 효율적인 제곱 시간의 단순 곱셈을 사용하는 것이 합리적입니다.

먼저 핵심 재귀 프로시저의 코드를 제시하고, 그다음에 이 코드가 왜 동작하는지 설명하겠습니다.

void karatsuba(int *a, int *b, int *c, int n) {
    if (n <= 64) {
        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                c[i + j] += a[i] * b[j];
    } else {
        int k = n / 2;
        int l[k], r[k], t[n] = {0};
        for (int i = 0; i < k; i++) {
            l[i] = a[i] + a[k + i];
            r[i] = b[i] + b[k + i];
        }
        karatsuba(l, r, t, k); // считает t
        karatsuba(a, b, c, k); // считает p1
        karatsuba(a + k, b + k, c + n, k); // считает p2
        int *t1 = t, *t2 = t + k;
        int *s1 = c, *s2 = c + k, *s3 = c + 2 * k, *s4 = c + 3 * k;
        for (int i = 0; i < k; i++) {
            int c1 = s2[i] + t1[i] - s1[i] - s3[i];
            int c2 = s3[i] + t2[i] - s2[i] - s4[i];
            c[k + i] = c1;
            c[n + i] = c2;
        }
    }
}

세 번의 재귀 호출이 끝나면 배열 $c$는 $p_1$과 $p_2$를 이어 붙인 것입니다.

그 후, 다항식 $c$ 자체를 계산하려면 $c$를 머릿속으로 네 개의 같은 부분으로 나누고, 다항식 $t$를 두 절반 $t_1$과 $t_2$로 나눈 다음, 식을 보면서 각 부분이 어떻게 바뀌는지 생각해 보는 것이 가장 간단합니다.

$s_1$: 변하지 않습니다 — 이는 $p_1$의 앞쪽 절반입니다.
$s_2$: $s_2 + t_1 - s_1 - s_3$로 표현됩니다. 즉, $t - p_1 - p_2$의 「첫 번째」 절반만큼 변합니다.
$s_3$: $s_3 + t_2 - s_2 - s_4$로 표현됩니다. 즉, $t - p_1 - p_2$의 「두 번째」 절반만큼 변합니다.
$s_4$: 변하지 않습니다 — 이것은 $p_2$의 후반부입니다.

벡터화 때문에 가능한 한 「가벼운」 데이터 타입을 사용하고, 가능하다면 AVX를 활성화하여 컴파일하는 것이 중요합니다:

#pragma GCC optimize("O3")
#pragma GCC target("avx2")

이 구현은 충분히 효율적이어서 크기가 $4 \cdot 10^5$인 두 다항식을 1초 안에 곱할 수 있습니다.

#아이디어의 일반화

비슷한 방법을 행렬 곱셈에도 적용할 수 있으며, 이를 슈트라센 알고리즘이라고 합니다. 이 알고리즘에서는 두 행렬을 $8 = 4 + 4$개의 부분으로 나누어 블록 단위로 곱하는데, 복잡한 대수 계산을 통해 8번의 곱셈 중 한 번을 없앨 수 있으므로 $O(n^{\log_2 7}) \approx O(n^{2.81})$의 점근적 복잡도를 얻습니다.

슈트라센 알고리즘과 카라츠바 알고리즘 모두 객체를 더 많은 부분으로 나누면 더 나은 점근적 복잡도를 달성할 수 있습니다. 그러나 실제로는 이러한 방법이 사용되지 않는데, 이런 알고리즘의 점근적 복잡도에는 비실용적으로 큰 상수가 숨어 있기 때문입니다.