흔히 컴퓨터의 성능을 테스트 해볼 때 벤치마크 프로그램을 사용합니다. 이런 벤치마크 프로그램은 여러 종류가 있으나 슈퍼파이나, wPrime과 같이 수학 상수를 계산하는 프로그램도 있습니다. 수학 상수라 하면 여러 가지가 있겠으나, 보통 계산되는 수는 원주율(π)이나 오일러 상수 e, 그리고 루트2 정도가 있습니다.

이들 상수들은 사실 공학이나 자연계에서 그렇게 긴 자리가 필요한 것은 아니기 때문에, 사실 컴퓨터로 긴 자릿수를 계산하는 것은 뻘짓 입니다. 그러나 이런 긴 자릿수 계산은 컴퓨터의 성능을 테스트하기 가장 쉬운 방법이며 계산하는 방법에 문제가 있는지, 혹은 하드웨어에 어떤 오류가 있는지 발견 할 수 있는 좋은 방법이기도 합니다. 오늘 이 글에서는 수학 상수들이 어떤 원리로 계산되는지 알아보겠습니다. 최대한 어려운 부분은 뺐으나 고등학생 이상이 읽기를 권장합니다.

1. 원주율 계산

원주율은 고대부터 수 천년동안 인간에 의해 계산되어진 중요한 상수중 하나입니다. 사실 여기 기글에 계신 분들은 적어도 초등학교 6학년 이상은 될 터이니 원주율이 어떤 숫자인지는 다들 잘 아실 겁니다. 처음 컴퓨터로 원주율을 계산 한 것은 흔히 최초의 전자식 컴퓨터라 불려지는 ENIAC(사실 ENIAC이 최초의 전자식 컴퓨터는 아닙니다.)으로, 원주율 소수점 이하 2037자리까지 계산 하는데 70시간정도 걸렸다고 합니다. 이후 컴퓨터 성능이 좋아지고, 수학이 발전함에 따라 원주율을 계산하는데 더 효율적인 식들과 방법들이 발견되면서, 계산시간이 극적으로 단축되었습니다.

우선 원주율은 어떤 방법으로 계산할까요? 원을 하나 그려가지고 지름 길이를 구한다음 원의 넓이를 나누는 원주율의 정의를 그대로 이용할까요? 그렇지 않습니다. 현실적으로 완전한 원은 없기 때문에 이런 방법을 이용하면 정확한 원주율 값을 구하기가 어렵습니다. 실제로 고대의 사람들은 이런 방법으로 원주율을 구했지만 아르키메데스가 겨우 소수점 이하 2자리까지 정확히 값을 구했을 뿐입니다. 그러나 17세기 이후 미적분이 발견되고 나서 상황이 바뀌었는데 바로 무한급수를 이용하여 계산하는 것입니다.

무한급수는 말 그대로 무한히 급수가 진행 되는 것입니다. 급수라는 것은 항과 항 사이에 덧셈 관계로 나타내진 수로써 다음과 같은 것들은 모두 무한급수입니다.

원주율에 관한 무한급수 식은 역 삼각함수를 미적분을 이용하면 쉽게 유도 할 수가 있는데 그리하여 발견된 원주율 무한급수 식은 엄청나게 많았습니다. ENIAC에서 원주율을 계산하던 때만 해도 이런 형태의 식을 사용했었습니다. 이들 공식은 항을 많이 더하면 더할수록 원주율 값에 일정한 속도로, 즉 값이 선형으로 수렴합니다. 다음 공식은 ENIAC에서 사용된 마친의 공식(Machin's formula)으로 역 탄젠트 함수를 이용한 급수입니다. 프로그래밍 하기 간단하여 컴퓨터 초창기에 계산용으로 많이 사용되었습니다.

저 식에 대입해서 계산한 다음 4배하면 정말로 원주율 값이 튀어나옵니다. 4arctan(1/5)는 하나의 항을 계산할 때 마다 1.3979자리씩 수렴하고, arctan(1/239)는 하나의 항을 계산할 때 마다 4.7558자리씩 수렴합니다. 현대의 공식에 비하면 그다지 빠른 공식이라고는 할 수 없습니다. 

시간이 흐르고 1970년대 쯤 되어 타원적분을 활용한 원주율 계산 공식이 나왔는데, 이는 슈퍼파이에 사용된 Gauss-Legendre 알고리즘입니다. 실제로 가우스와 르장드르가 만든 공식은 아니고 이 수학자들의 연구 결과를 토대로 만든 공식인데 원리를 보면 다음과 같습니다.

1. 우선 초기값을 다음과 같이 줍니다.

2. 다음과 같은 계산을 반복합니다.

3. 적당히 반복 한 후에 원주율은 다음 식으로 구합니다.

이 공식의 특징은 한번 계산 할 때 마다, 이전에 비해 정확한 자릿수가 2배가 된다는 것입니다. 앞에서 본 무한급수들은 모두 일정한 속도로 가까워 졌지만 이것은 2배씩 가까워지므로 매우 적은 반복 횟수로 긴 자릿수의 원주율을 쉽게 구할 수가 있습니다. 이를 어려운 말로 2차 수렴(2nd-order convergence)이라 하는데, 슈퍼파이가 1M자리를 계산할 때 단 19번만 루프를 도는 것은 바로 그 때문입니다.

하지만 최근에 와서는 이 공식을 잘 이용하지 않습니다. 각 항의 연계가 심해서 멀티스레드화 하기도 쉽지 않고 가장 빠른 공식도 아니기 때문입니다. 현재 원주율을 계산하는데 이용되는 가장 빠른 공식은 러시아 수학자 형제인 Chudnovsky 형제가 1987년도에 만든 공식으로써 대략 한 개의 항을 계산 할 때 마다 14.181647자리를 계산할 수 있습니다.

하지만 이 공식도 앞에 보았던 무한급수들에 비해서 꼴이 좀 복잡해졌다 뿐이지 선형수렴을 하는데 어찌하여 Gauss-Legendre 공식보다 더 빨리 계산될 수 있는 것일까요. 우선 공식 자체도 항 자체가 독립적이기 때문에 멀티스레드화 하기 편하고, 이분 재귀 연산법(Binary splitting method)을 사용하면 계산량 자체도 줄어드는데다 급수 식을 정수의 팩토리얼과 곱셈 형태로 바꿀 수 있으므로 나눗셈 연산을 거의 할 필요가 없는데다가(나눗셈은 곱셈에 비해 매우 느립니다.) 재귀 호출로 작업을 쪼갤 수 있어서 속도가 극적으로 빨라지게 됩니다. 여담이지만 DJ PI 2.0의 멀티스레드 계산은 저 공식을 이용한 것입니다. 이 공식은 최근에 나온 대부분의 원주율 계산 프로그램이나 mathematica와 같은 프로그램에서도 사용됩니다.

2. e의 계산

e는 자연로그에서 밑으로, 그리고 미적분에서 자주 쓰이는 수학 상수입니다. e는 원주율에 비하면 매우 계산하기 쉬운 상수입니다. e는 정의를 이용하는 것 보다는 매클로린급수로 유도한 이 급수를 이용하여 x=1 을 만들어 놓고 계산을 합니다.

이 식은 발견 된지는 꽤 오래됐지만, 수렴하는 속도가 매우 빨라서 현재도 잘 사용되는 공식입니다. 하지만 이 공식도 그냥 계산하는 것 보다는 앞에서 언급한 이분 재귀 연산법을 이용하면 더 빠르게 계산할 수 있습니다. e는 비교적 쉬운 계산 공식과 빠른 수렴속도 때문에 사람들이 그다지 많이 도전하는 것 같지는 않습니다.

3. 제곱근의 계산

제곱근의 경우에는 앞에서 본 것과는 좀 다릅니다. 그래도 원주율에 비해서 수렴 속도도 빠르고 그 공식자체는 매우 간단합니다. 제곱근은 그 값을 계산할 만한 적당한 급수 식(제곱근의 급수 식은 수렴이 너무 느립니다.)이 없기 때문에, 빠른 계산을 위해서 주로 뉴턴-랩슨 법을 이용합니다. wPrime 역시 이 방법을 사용합니다.

뉴턴-랩슨 법은 쉽게 말하자면 방정식의 근을 미분을 이용해 구하는 것입니다. 미분을 하여 구한 접선이 점점 근에 가까워지는 원리를 이용하는데, 예를 들어서 루트2를 구하기 위해서는 f(x)=x^2-2=0의 근을 뉴턴-랩슨 법을 이용하여 구하면 됩니다. 뉴턴-랩슨 법은 제곱근 뿐 아니라 고차방정식의 해를 구하는데도 잘 사용됩니다.

다음과 같은 함수가 있을 경우, 접선이 x축과 만나는 xn+1에 수선을 긋고 또 접선을 긋고 하는 식으로 함수의 근에 접근하는데

이것을 공식으로 풀이해 보자면, 다음과 같습니다.

뭐 f(xn)자리에 루트2를 계산하고 싶다면 x^2-2를 집어넣으면 됩니다. 귀찮으니까 미리 정리해놓은 점화식을 봅시다.

다음 식은 a를 대입하면 루트a를 계산하는 점화식입니다.

뉴턴-랩슨 법의 특징이라 한다면 점화식 하나로 값을 구할 수 있으므로, 계산 속도가 빠르다는 것입니다. 그리고 Gauss-Legendre 공식처럼 2차 수렴의 특징 또한 갖습니다. 따라서 항을 하나 계산할 때마다 정확한 자리는 두배씩 늘어납니다.

그럼 실제로 계산이 되는지 한번 비교 해 봅시다. a에 2를 집어넣고 x0 = 1을 넣고 계산 해 보면 (초기값은 아무것이나 해도 상관없습니다.)

x0 = 1

x1 = 1.5

x2 = 1.4166666666666666666666666666667

x3 = 1.4142156862745098039215686274509

x4 = 1.4142135623746899106262955788902

x5 = 1.4142135623730950488016896235025

x6 = 1.4142135623730950488016887242097

실제 루트 2 = 1.4142135623730950488016887242097....

반복을 6번 밖에 안했는데 실값에 무려 31자리나 가까워졌습니다.  

4. 부동 소수점 수의 한계와 계산을 위한 다른 방법

컴퓨터가 이들 상수를 계산하기 위해서는 적법한 데이터 형이 필요합니다. 흔히 컴퓨터에서 사용하는 부동 소수점 수 형식은 IEEE754로 double형의 경우 그 사이즈가 8바이트입니다. 하지만 이 Double형 하나로 100만자리씩 값을 구할수는 없습니다.

컴퓨터는 소수점 이하 수도 이진법으로 표현합니다. 정수 부분은 사실 어떤 진법으로 표현하더라도 정밀도의 문제가 없으나, 소수점 이하의 수를 표현하게 될 경우 진법간의 차이가 생깁니다. 예를 들자면 십진 소수 0.1은 이진수에서는 무한소수가 되어버립니다.

double형의 경우 유효숫자부가 53bit로, 결론적으로 정확하게 표현될 수 있는 자리는 1.11022e-16까지이므로 대략 소수점 이하 15자리까지 입니다. 이정도 정밀도로는 100만자리는 꿈도 꿀 수 없습니다.

따라서 새로운 자료구조를 정의 해야 합니다. 문제는 덧셈 뺄셈은 쉬울것 같지만 곱셈 나눗셈은 그렇지 않다는 것입니다.

곱셈 나눗셈은 초등학교에서 배운 방법대로 정의할 경우 시간복잡도가 커서, 계산속도에 지대한 영향을 끼치게 됩니다. 따라서 곱셈은 고속 푸리에변환을 이용하여 계산하고, 나눗셈은 뉴턴-랩슨법을 이용하여 역수를 구하는 방법을 이용하면 속도를 빠르게 할 수 있습니다.

문제는 고속 푸리에변환을 하는 프로그램을 구현하는게 공대 수학을 배우지 않은 일반인에게는 만만치 않습니다. 따라서 잘 짜여진 라이브러리를 이용하는데, GMP(Gnu Multi-Precision Library)가 그 대표적인 라이브러리 중 하나입니다. 이 라이브러리는 mathematica와 같은 프로그램에서도 활용 되었으며 또 여담이지만 DJ PI v2.0에서도 사용되었습니다.