큰 자리수의 곱셈 알고리즘인 카라츠바 알고리즘에 대해서 핵심만 정리해본다.

약간의 수학

구하고자 하는 곱셈이 $x * y$ 인 경우를 생각하자. 여기서 $x$, $y$는 자연수이다. 편의상 두 수의 자리수가 같다고 하자. (만약 같지 않은 경우라면, 작은 녀석을 큰 녀석의 크기로 맞추고 padding을 추가하면 된다) 그러면, $x$, $y$는 아래처럼 다시 쓸 수 있다.

• $x = x_1B^m + x_0$
• $y = y_1B^m + y_0$

수식으로 보면 복잡해 보이는데, 이는 단순히 $123456 = 123 * 10^3 + 456$ 처럼 특정 수를 high, low로 나눈 것 뿐이다. high에 해당하는 것이 $x_1$, $y_1$이고, low에 해당하는 것이 $x_0$, $y_0$ 이다. 이를 formal하게 나타낸 것이 위 수식이고, $B$는 진법을 나타내고, $m$은 큰 수의 총 자릿수가 $l$이라고 할 때 절반으로 나눈 값이다. ($m = l / 2$, 절반으로 나누면 좋지 않겠는가?) 그러면 구하고자 했던 곱셈은 아래처럼 다시 쓸 수 있다.

이것을 이대로 계산하면, 곱셈의 구간을 절반으로 줄인 대신, 곱셈이 횟수가 4번으로 늘어났기 때문에 시간복잡도 계산을 해보면 $O(N^2)$이 되기 때문에 Naive하게 구했던 경우랑 다를게 없다. 그런데, 만약에 위 수식을 계산하는데에 곱셈의 횟수가 3번만 필요하다면 어떨까? 자리수가 작은 경우에는 위 방식이 상수항을 늘리기 때문에 크게 의미가 없지만, 자리수가 클 경우($N$이 클 경우)에는 충분히 의미가 있다. 카라츠바 알고리즘의 핵심은 $B^m$의 계수에 붙어있는 녀석을 아래처럼 변형하는데에 있다.

• $x_1y_0 + x_0y_1 = (x_1 + x_0)(y_1 + y_0) - x_1y_1 - x_0y_0$

곱셈을 줄인다더니 두번에서 오히려 세번으로 늘리고 있네 라고 할 수 있다. 그런데 식을 자세히 보자. $x_1y_1$과 $x_0y_0$은 어차피 계산해야할 값이었다. 즉, 추가로 계산해야할 곱셈은 $(x_1 + x_0)(y_1 + y_0)$ 하나다.

실제 구현

위 알고리즘은 간단하게 재귀로 아래처럼 구현이 가능하다. $B^m$ 부분은 특정 자리수 부분에다가 더하는 것이므로, 배열에 잘 접근하는 것으로 처리가 가능하다. 또한, 카라츠바는 상수항이 큰 알고리즘이기 때문에, 실제로 써먹을만한 구현은 naive한 곱셈 구현을 섞어서 한다.

// naive 한 곱셈 구현, 시간복잡도 O(N^2)
void naive_mult(ull r[], ull x[], ull y[], int n) {
    for(int i = 0; i < n; ++i) {
        ull cur = x[i];
        ull* rp = &r[i]; ull* yp = &y[0];
        for(int j = 0; j < n; ++j, ++rp, ++yp) {
            *rp += *yp * cur;
        } 
    }
}
void karatsuba(ull r[], ull x[], ull y[], int n) {
    if (n < 50) {
        naive_mult(r, x, y, n); return; // 일정 자리수 이하로 내려가면 karatsuba는 비효율적이다. naive한 곱셈을 쓰자
    }
    int m = n >> 1; // N는 가능한한 N / 2 딱 떨어지게되는 가장 좋은 수를 선택하자.
    ull xs[SZ], ys[SZ], tmp[SZ];
    for(int i = 0; i < m; ++i) xs[i] = x[i] + x[i + m], ys[i] = y[i] + y[i + m];
    for(int i = 0; i < n; ++i) tmp[i] = 0;

    // 4개의 곱셈 -> 3개의 곱셈으로 줄어듬
    karatsuba(tmp, xs, ys, m);
    karatsuba(r, x, y, m);
    karatsuba(r + n, x + m, y + m, m); // B^2m 구간의 결과는 여기로 들어간다.

    for(int i = 0; i < n; ++i) tmp[i] -= r[i] + r[i + n]; // (x0 + x1)(y0 + y1) - x0y0 - x1y1
    for(int i = 0; i < n; ++i) r[i + m] += tmp[i]; // * B^m
}

Time Complexity

이 알고리즘의 구현을 살펴보면, 분할정복을 하는 경우의 시간복잡도 계산을 이용하여 시간복잡도를 구할 수 있다. 그럴듯하게 알고리즘을 어떻게 분할하면서 풀고 있느냐 하면, 구간이 $N / 2$로 줄어들 때마다, $3$개의 하위 문제가 생기고 있다. 또한, 이러한 하위문제를 합치는데 걸리는 시간이 $f(N) = 2N$이다. 이 경우는, 인 경우에 속하기 때문에, 시간복잡도는 이 문서에서 본 것처럼 계산할 수 있다. 결과만 보면,

이 된다.