7장 분할 정복

재귀호출을 이용하여 \(n\)까지의 합을 계산하는 것보다 더 빠른 합 함수를 알아본다. n까지의 합을 구하는 함수 \(fastSum()\)은 다음과 같이 표현할 수 있다.

\begin{align} fastSum(n)&=1+2+\cdots+n\\&=\left(1+2+\cdots+\frac{n}{2}\right)+\left(\left(\frac{n}{2}+1\right)+\cdots+n\right)\\\end{align}

합을 두 부분으로 분리하면, 첫 번째 부분은 \(fastSum\left(\frac{n}{2}\right)\)으로 나타낼 수 있지만, 두 번째 부분은 그럴 수 없다. 다음은 두 번째 부분 식을 어떻게 정리하는지 보여준다.

\begin{align} \left(\frac{n}{2}+1\right)+\cdots+n&=\left(\frac{n}{2}+1\right)+\left(\frac{n}{2}+2\right)+\cdots+\left(\frac{n}{2}+\frac{n}{2}\right)\\&=\frac{n}{2}\times\frac{n}{2}+\left(1+2+3+\cdots+\frac{n}{2}\right)\\&=\frac{n}{2}\times\frac{n}{2}+fastSum\left(\frac{n}{2}\right)\\\end{align}

이 결과를 가지고 첫 번째 부분의 결과이었던 \(fastSum\left(\frac{n}{2}\right)\) 와 더하면 다음과 같은 식이 된다.

\[fastSum(n)=2\times fastSum\left(\frac{n}{2}\right)+\frac{n^2}{4}\quad(\text{when}\ n\ \text{is even})\]

다음은 이것을 구현한 코드이다. 만약 \(n\)이 홀수인 경우 코드에서 어떻게 처리하는지 잘 살펴보자.

int fastSum(int n) 
{
	if(n == 1) return 1;
    // n이 홀수라면
	if(n % 2 == 1) return fastSum(n-1) + n;
	return 2 * fastSum(n/2) + (n/2) * (n/2);
}

 

행렬의 거듭제곱 알고리즘

class SquareMatrix{};

SquareMatrix identity(int n);

SquareMatrix pow(const SquareMatrix& A, int m) 
{
	if (m == 0) return identity(A.size());
	if (m % 2 > 0) return pow(A, m - 1) * A;
	SquareMatrix half = pow(A, m / 2);
	return half * half;
}

카라츠바의 빠른 곱셈 알고리즘

만약 \(a\)와 \(b\)가 각각 256자리 수라면 \(a_1\)과 \(b_1\)은 첫 128자리, \(a_0\)과 \(b_0\)은 그다음 128자리를 저장한다고 가정하자. 그러면 다음과 같이 표현할 수 있다.

\[a=a_1\times10^{128}+a_0\]

\[b=b_1\times10^{128}+b_0\]

두 수의 곱은 다음과 같이 표현할 수 있다.

\begin{align}a\times b&=(a_1\times10^{128}+a_0)\times(b_1\times10^{128}+b_0)\\&={\color{Red}a_1\times b_1}\times 10^{256}+({\color{Red}a_1\times b_0}+{\color{Red}a_0\times b_1})\times10^{128}+{\color{Red}a_0\times b_0}\end{align}

단순하게 분배 법칙을 사용하여 전개했다. 위 식은 총 4번의 부분 숫자들의 곱셈으로 표현된다(빨간색으로 표시한 부분). 10의 거듭제곱과의 곱은 시프트 연산으로 구현하므로 곱셈 연산으로 보지 않는다. 그리고 \( z_0, z_1, z_2\)을 다음과 같다고 하자.

\[a\times b=\underbrace{a_1\times b_1}_{z_2}\times10^{256}+\underbrace{(a_1\times b_0+a_0\times b_1)}_{z_1}\times10^{128}+\underbrace{a_0\times b_0}_{z_0}\]

위 식에서 두 번의 곱셈으로 \( z_0, z_2 \)를 각각 구하고 다음 식을 이용한다.

\begin{align}(a_0+a_1)\times(b_0+b_1)&=\underbrace{a_0\times b_0}_{z_0}+\underbrace{a_1\times b_0+a_0\times b_1}_{z_1}+\underbrace{a_1\times b_1}_{z_2}\\&=z_0+z_1+z_2\end{align}

위 식의 결과에서 \(z_0\)과 \(z_2\)를 빼서 \(z_1\)을 구할 수 있다.

\[z_1=(a_0+a_1)\times(b_0+b_1)-z_0-z_2\]

이렇게 기존에 총 4번의 곱셈을 써야 할 것을, 3번의 곱셈으로도 구할 수 있게 되었다. 코드로 구현하면 다음과 같다. 정수의 길이가 50자리 미만이면 \(O(n^2)\) 곱셈 알고리즘을 사용하는 부분도 잘 봐두길 바란다.

#include <iostream>
#include <algorithm>
#include <vector>
#include <cmath>
using namespace std;

// num[]의 자릿수 올림을 처리한다.
void normalize(vector<int>& num)
{
	num.push_back(0);
	for (int i = 0; i + 1 < num.size(); ++i) {
		/* 
        	karatsuba의 곱셈 알고리즘에서 필요한 부분
            두 수의 뺄셈 연산이 있어 뺀 결과가 음수가 나올 수 있다.
            음수가 나온다면 앞의 자리 숫자에서 숫자를 빌려와
            해당 자리의 숫자를 다시 양수로 만들어준다.
        */
		if (num[i] < 0) {
			int borrow = (abs(num[i]) + 9) / 10;
			num[i + 1] -= borrow;
			num[i] += borrow * 10;
		}
		// num[i] >= 0이다. 직관적인 부분.
		else {
			num[i + 1] += num[i] / 10;
			num[i] %= 10;
		}
	}
	// 숫자의 자릿수가 1보다 크고 현재 최고자리 숫자가 0이라면
	// (vector에 가장 뒤에있는 숫자가 최고자리 숫자임, 앞에서 0넣음)
	while (num.size() > 1 && num.back() == 0)
		num.pop_back();
}

// O(n^2) 곱셈 알고리즘.일반적이다.
// 두 긴 자연수의 곱을 반환한다.
// 각 배열에는 각 수의 자릿수가 1의 자리에서부터 시작해 저장되어 있다.
// a[0]는 1의 자리 숫자. a[max]는 최고자리 숫자. b도 마찬가지
vector<int> multiply(const vector<int>& a, const vector<int>& b)
{
	vector<int> c(a.size() + b.size() + 1, 0);
	for (int i = 0; i < a.size(); ++i)
		for (int j = 0; j < b.size(); ++j)
			c[i + j] += a[i] * b[j];

	normalize(c);
	return c;
}

// a+=b*(10^k)를 구현한다.
void addTo(vector<int>& a, const vector<int>& b, int k)
{
	// a가 b보다 크기 때문에 b만 모두 순회하여 a에 더한다.
	int i = 0;
	int ten_k = static_cast<int>(pow(10, k));
	for (auto iter = begin(b); iter != end(b); ++iter)
		a[i++] += *iter * ten_k;
    normalize(a);
}

// a-=b를 구현한다. a>=b를 가정한다.
void subFrom(vector<int>& a, const vector<int>& b)
{
	int i = 0;
	for (auto iter = begin(b); iter != end(b); ++iter)
		a[i++] -= *iter;
    normalize(a);
}

// 두 긴 정수의 곱을 반환한다.
vector<int> karatsuba(const vector<int>& a, const vector<int>& b)
{
	// 각 수가 몇자리 수인지. 자릿수 길이
	int an = a.size(), bn = b.size();
	// a가 b보다 짧을 경우 둘을 바꾼다. 따라서 앞이 항상 더 크다.
	if (an < bn) return karatsuba(b, a);
	// 기저 사례: a나 b가 비어 있는 경우, 빈 벡터를 반환한다(0이므로).
	if (an == 0 || bn == 0) return vector<int>();
	// 기저 사례: a가 비교적 짧을 경우 O(n^2) 곱셈으로 변경한다.
	if (an <= 50) return multiply(a, b);

	// 다음은 자릿수의 절반이지, 숫자의 절반이 아니다.
	int half = an / 2;
	// a와 b를 밑에서 half자리와 나머지로 분리한다. 
	// a = a_1 * 10^half + a_0
	// b = b_0 * 10^half + b_0
	vector<int> a0(begin(a), begin(a) + half); // 구간에서 뒤 원소는 포함 안함.
	vector<int> a1(begin(a) + half, end(a));
	vector<int> b0(begin(b), begin(b) + min<int>(b.size(), half)); // b가 half자리보다 더 작은 수일 수도 있음.
	vector<int> b1(begin(b) + min<int>(b.size(), half), end(b));

	vector<int> z2 = karatsuba(a1, b1);
	vector<int> z0 = karatsuba(a0, b0);
	addTo(a0, a1, 0); addTo(b0, b1, 0);
	vector<int> z1 = karatsuba(a0, b0);
	subFrom(z1, z0); subFrom(z1, z2);
	vector<int> ret;
	addTo(ret, z0, 0);
	addTo(ret, z1, half);
	addTo(ret, z2, half + half);
	return ret;
}