피보나치 수열

Reference

• Algorithm (Sanjoy Dasgupta)

---

피보나치 수열

이번에 아주 간단한 피보나치 수열과 n번째 피보나치 수를 구하는 알고리즘에 대해서 알아보도록 하곘습니다.

 

피보나치 수열은 위 이미지에 나타난 수열이며, 수열 내에서 각 숫자는 이전 두 숫자의 합입니다.

정확하게 이야기하면, 피보나치수 \(F_n\)은 아래의 간단한 규칙(점화식)으로 생성됩니다.

\[F_n = \left\{\begin{matrix}
F_{n-1} + F_{n-2} && \text{if } n > 1 \\ 
1 && \text{if } n = 1 \\ 
0 && \text{if } n = 0
\end{matrix}\right.\]

 

피보나치 수열은 \(2^n\)과 거의 같은 속도로 증가합니다. 예를 들어, \(F_{30}\)은 백만이 넘고, \(F_{100}\)은 21자리의 수입니다. 일반적으로 \(F_n = 2^{0.694n}\) 입니다.

 

그럼 \(F_{30}\)이나 \(F_{100}\)의 정확한 값은 얼마일까요 ?

이 질문에 답을 하기 위해서는 n번째 피보나치 수를 계산하는 알고리즘이 필요합니다.

 

---

An exponential algorithm

가장 먼저 \(F_n\)의 재귀적 정의를 그대로 사용하여 구하는 방법이 있습니다.

의사코드를 사용하여 알고리즘을 표현하면 다음과 같습니다.

C++ 코드로 표현하면 다음과 같습니다.

uint64_t fib1(size_t n)
{
	if (n == 0)
		return 0;
	if (n == 1)
		return 1;

	return fib1(n - 1) + fib1(n - 2);
}

이 알고리즘은 피보나치 수열의 정의를 가지고 그대로 만든 알고리즘이기 때문에 문제가 되는 건 없습니다.

다만, 이 알고리즘은 n의 값에 따라서 어떤 수행 시간을 가질까요?

 

이를 계산하기 위해서 \(T(n)\)을 \(fib1(n)\)을 계산하기 위한 기본적인 컴퓨터 연산의 개수라고 해봅시다.

n이 만약 2보다 작다면 몇 개의 연산만 수행하고 바로 알고리즘은 종료됩니다.

따라서 다음이 성립합니다.

\[T(n) \leq 2 \text{ for } n \leq 1\]

만약 n이 보다 커지면 \(T(n - 1)\)과 \(T(n - 2)\)만큼의 시간이 걸리는 \(fib1\)의 재귀호출을 각각 해야하고, 세 번의 연산(n의 값 체크와 마지막 덧셈)을 더 해야 합니다. 따라서 다음이 성립합니다.

\[T(n) = T(n - 1) + T(n - 2) + 3 \text{ for } n > 1\]

 

이를 \(F_n\)의 점화식과 비교해보면, \(T(n) \geq F_n\)이라는 사실을 알 수 있습니다.

이는 알고리즘의 수행 시간은 피보나치 수열값만큼 빨리 늘어난다는 의미이고, n이 매우 작지 않다면 알고리즘은 쓸모없을 만큼 느리다는 것을 뜻합니다.

즉, 지수적인 수행 시간이 필요한데, 예를 들어, \(F_{200}\)을 계산하기 위해서 위 알고리즘은 \(T(200) \geq F_{200} \geq 2^{138}\) 만큼의 기본 연산을 수행해야 합니다.

 

\(fib1(n)\)의 수행 속도는 \(2^{0.694n) \approx (1.6)^n\)에 비례하므로 \(F_{n+1}\)을 계산하는 데 \(F_n\)보다 1.6배의 시간이 걸립니다.

 

결과적으로 정의에 따라 만든 알고리즘은 올바른 결과를 찾기는 하지만 너무 느립니다.

이보다 빠른 다항 시간을 만족하는 알고리즘을 찾아봅시다.

 

---

A polynomial algorithm

fib1 알고리즘을 살펴보면, 연쇄적으로 호출되는 재귀호출들이 너무 많습니다.

즉, 수많은 계산이 반복되고 있습니다.

 

이것보다 합리적인 방법은 중간 결과들(\(F_0, F_1, \cdots, F_{n-1}\))을 나오는 대로 저장하는 것입니다.

(Dynamic Programming 이라고도 합니다.)

 

uint64_t fib2(size_t n)
{
	if (n == 0)
		return 0;

	std::vector<uint64_t> arr(n + 1, 0);
	arr[0] = 0;
	arr[1] = 1;

	for (size_t i = 2; i <= n; i++)
		arr[i] = arr[i - 1] + arr[i - 2];

	return arr[n];
}

 

fib1에서처럼 \(F_n\)의 정의를 그대로 사용하고 있기 때문에 알고리즘에 문제는 없습니다.

그렇다면 이 알고리즘의 수행 시간은 얼마일까요?

n-1번 반복되는 루프 안에는 컴퓨터 연산이 하나만 수행됩니다. 그러므로 fib2가 수행하는 컴퓨터 연산의 개수는 n에 선형적입니다. 지수적 수행시간의 알고리즘에서 다항 시간 알고리즘으로 엄청난 속도 상향이며, 이제 \(F_{200}\)을 구하는데에는 큰 무리가 없을 것입니다.

---

하지만 조금 더 정확하게 분석해보면, 조금 다르긴 합니다.

이때까지 기본적인 컴퓨터 연산들이 일정 시간이 걸린다고 가정하고 각각의 알고리즘들이 수행하는 기본적인 컴퓨터 연산의 개수를 세었습니다. 이는 매우 유용한 단순화방법입니다. 여기서 말하는 연산은 분기, 메모리에 저장, 숫자 비교 그리고 산술 연산 등 다양항 기본 연산들을 모두 포함합니다. 이러한 연산들을 모두 구분하기보다는 하나로 묶어서 생각하는 것이 편리하긴 합니다.

 

하지만 피보나치 수열 알고리즘을 분석할 때에는 너무 넓은 범위를 기본 연산으로 인정했습니다. 예를 들어 32비트로 표현될 수 있는 작은 숫자의 덧셈은 하나의 컴퓨터 수행 단계로 생각하는 것이 타당할 수 있습니다. 하지만 n번째 피보나치 수는 0.694n비트이기 때문에 n이 커지면 32비트를 넘을 수 있습니다.

따라서 임의의 크기를 가지는 숫자에 대한 산술 연산은 상수 시간에 수행될 수 없고, 앞에서 보았던 수행 시간을 엄밀히 따져보아야 조금 더 정확한 분석이 됩니다.

 

후에 기회가 된다면 자세히 살펴보겠지만, n비트의 두 숫자를 더하는 데 걸리는 시간은 n에 비례합니다.

(덧셈을 한자리씩 하던 시절을 떠올려보면 납득하기 쉽습니다.. !)

 

fib1은 \(F_n\)만큼의 덧셈을 수행하므로 이는 \(nF_n\)개의 기본 연산에 해당하고, 마찬가지로 fib2가 수행하는 기본 연산의 개수는 \(n^2\)에 비례합니다. 이는 아직도 다항 시간이므로 지수적 수행 시간이 걸리는 fib1보다는 훨씬 빠릅니다.

---

행렬을 이용한 연산

피보나치 수열은 아래처럼 행렬식으로 나타낼 수도 있습니다.

\[\begin{align*}
\begin{bmatrix}F_{n}\\ F_{n+1}\end{bmatrix} &= \begin{bmatrix}0 & 1\\ 1 & 1\end{bmatrix}\begin{bmatrix}F_{n-1}\\ F_{n}\end{bmatrix}\\ 
 &= \begin{bmatrix}0 & 1\\ 1 & 1\end{bmatrix}^{n} \begin{bmatrix}F_{0}\\ F_{1}\end{bmatrix}\\ 
 &= \begin{bmatrix}0 & 1\\ 1 & 1\end{bmatrix}^{n} \begin{bmatrix}0\\ 1\end{bmatrix}
\end{align*}\]

 

즉, \[\begin{bmatrix}0 & 1\\ 1 & 1\end{bmatrix}^{n}\]

만 계산하여, 1행 2열의 값만 취하면 \(F_n\)을 구할 수 있습니다.

 

이는 C++로 아래와 같이 fib3 함수로 구현할 수 있습니다. 

template<size_t R, size_t C>
struct Matrix
{
	Matrix() {
		mat[0][0] = mat[1][1] = mat[0][1] = mat[1][0] = 0;
	}
	Matrix operator*(const Matrix& rhs) const {
		Matrix<R, C> ret;
		
		ret.mat[0][0] = this->mat[0][0] * rhs.mat[0][0] + this->mat[0][1] * rhs.mat[1][0];
		ret.mat[0][1] = this->mat[0][0] * rhs.mat[0][1] + this->mat[0][1] * rhs.mat[1][1];
		ret.mat[1][0] = this->mat[1][0] * rhs.mat[0][0] + this->mat[1][1] * rhs.mat[1][0];
		ret.mat[1][1] = this->mat[1][0] * rhs.mat[0][1] + this->mat[1][1] * rhs.mat[1][1];

		return ret;
	}

	uint64_t mat[R][C];
};
using Mat2by2 = Matrix<2, 2>;

Mat2by2 MatExp(Mat2by2 mat, size_t n)
{
	if (n == 1)
		return mat;

	Mat2by2 half = MatExp(mat, n / 2);
	
	return (n % 2) ? (half * half * mat) : (half * half);
}

uint64_t fib3(size_t n)
{
	if (n == 0)
		return 0;

	Mat2by2 base;
	base.mat[0][0] = base.mat[1][1] = 1;
	Mat2by2 ret;
	ret.mat[0][0] = 0;
	ret.mat[0][1] = ret.mat[1][0] = ret.mat[1][1] = 1;

	ret = MatExp(ret, n);

	return ret.mat[0][1];
}

행렬의 곱셈에 여러 번의 연산이 있지만, 행렬 곱셈을 컴퓨터 연산 한 번이라고 가정한다면, fib3이 수행하는 산술 연산의 개수는 fib2의 \(O(n)\)보다 훨씬 적은 \(O(logn)\) 입니다. 이는 행렬의 n승을 구할 때, 이전에 구한 n/2승의 값을 이용하기 때문입니다.

 

하지만, 행렬 곱셈을 컴퓨터 연산 한 번이라고 가정한다고 하더라도, 덧셈이 아닌 곱셈을 포함한다는 것을 생각하면 fib2보다 빠르다고 확신은 할 수 없을 것입니다. 참고로, 일반적으로 큰 수들 간의 곱셈은 덧셈보다 느립니다.

(행렬의 곱은 이보다 더 느리겠죠. 정확한 계산은 책을 더 읽어보고 가능하다면 계산해보도록 하겠습니다.. !)

 

---

공식 이용

마지막으로 다음과 같은 피보나치 수열을 계산하는 공식이 있습니다.

\[F_n = \frac{1}{\sqrt5}(\frac{1+\sqrt5}{2})^n - \frac{1}{\sqrt5}(\frac{1-\sqrt5}{2})^n\]

 

double fib4(size_t n)
{
	double sqrt5 = sqrt(5);
	double ret = (1 / sqrt5) * pow((1 + sqrt5) / 2, n) - (1 / sqrt5) * pow((1 - sqrt5) / 2, n);
	return ret;
}

 

공식을 사용하면, 상수 시간 내에 피보나치 수열을 구할 수가 있습니다. 

다만, 사용되는 수들이 무리수이기 때문에 정확하게 정수로 떨어지지 않으므로 주의해야 합니다.