원주율(pi)을 계산하자
어려운 문제가 가끔 날 유혹해, 그래도 난 좋다네.

차 례

  1. 도입 배경
  2. 개념 파악
  3. 수학적 분석
  4. 알고리즘 연구
  5. 원주율 계산 식
  6. 프로그래밍

1. 도입 배경

 대학에 입학(91년도)하고 학과 컴퓨터실에서 286을 마음대로 사용할 수 있었다. 이 때까지 GW-BASIC 프로그램을 주로 사용하던 나는 실행파일을 만들 수 있는 Quick Basic을 듣게 되어 공부를 했다. 하드가 없던 286(XT)에서 도스 디스켓과 퀵베이직 디스켓을 들고 다니며 이것 저것 프로그램을 짜던 나는 베이직 명령어와 짧은 수학적 지식으로 테트리스 프로그램을 짤 수 있게 되었다. 그 외에도 퍼즐형 게임을 만들며 흐뭇해 하곤 했다. 하지만 학기초부터 궁금했던 원주율(pi, 이하 pi라고 한다.)의 계산은 여전히 머릿속을 떠나지 않는 스스로의 숙제였다.

주위의 선배들이나 교수님 중에서도 완전한 해답을 주시는 분은 없었다. 수학 문제 중에는 풀리는 문제와 미해결 문제가 있는데, pi를 소수 100자리까지 구하는 프로그램을 짜는 것이 나에게 미해결 문제가 될까, 몇 달 동안 고민을 했다. 결국 어설프지만 소수 1500자리까지 계산하는 프로그램이 성공했을 때 무척 기뻤다.

이후에 군 전역하고 복학하여 TC++을 공부하며 C 언어로도 프로그램을 다시 짰는데 해를 구하는 속도가 훨씬 빨라졌다. 알고리즘에서 약간 개선된 점도 있고, 컴퓨터도 486이었기 때문이다. 하지만 여전히 깔끔하게 정리되지 않고 그냥 주먹구구식이 된 것같아 이번에 그 알고리즘을 정리도 할 겸, 여러분들의 도움도 받을 겸 해서 공개하게 되었다. 필자가 한 방법에 문제가 있거나 더 좋은 알고리즘이 있으면 반드시 게시판으로 연락을 주면 고맙겠다.

2. 개념 파악

원주율은 지름의 길이에 대한 원주의 길이의 비율이다.
원주는 원의 둘레의 길이이다. 원이란 한 지점으로부터 일정한 거리(반지름)에 있는 점들의 집합이다. 점은 위치만 있고 크기가 없으며, 선은 길이만 있고 두께는 없다. 원주율은 지름의 길이에 관계없이 일정하다(구부러지지 않은 평평한 2차원 평면 위에서만 다룬다).

원주율은 실수이다. 원주율은 유리수가 아니다. 원주율은 무리수이다. 원주율은 무한 개의 유리수의 합으로 나타낼 수 있다. 유리수는 분수(정수/정수 꼴) 형태로 나타낼 수 있으며, 무리수는 단일 분수로 나타낼 수 없다. 소수는 유한소수와 무한소수가 있으며 모든 유한소수는 무한소수로도 표현할 수 있다. 무한소수 중 순환하지 않는 것은 분수 형태로 나타낼 수 없으며 이들을 무리수라 한다.

원주율은 일정한 값을 가진다.
원주율은 3보다 크고 4보다 작다. 이것은 고대에 발견되었는데, 원에 외접하는 정사각형을 그림으로써 4보다 작음을, 원에 내접하는 정육각형을 그림으로써 3보다 큼을 쉽게 증명할 수 있다.
일정한 지름의 길이를 가지는 원(혹은 원통)을 굴려 그 길이를 재 봄으로써 원주율의 근사치를 측정할 수 있지만 이는 정확한 값이 될 수 없다. 원에 내.외접하는 n각형(n이 클수록 오차가 작아짐)을 계산(피타고라스의 정리를 이용)함으로써 근사치를 구할 수도 있지만 제곱근 역시 무리수라서 나중엔 어려움을 느낀다.삼각함수나 기타 다른 방법으로 계산하려 할 때도 결국엔 어려운 무리수들의 합으로 표현되는 경우가 많다.

그 수많은 원주율 계산법 중에서 필자가 계산한 방법을 소개하고자 한다.
이는 필자 혼자만의 머리에서 나온 것이 아니라 여러 수학 관련 서적을 보고 연구한 것임을 밝혀 둔다. 당연한 이야기지만 내장되거나 툴에서 표현하는 pi를 사용하는 것이 아니라 완전히 새로운 pi를 구하는 프로그램을 짜야 한다. 필자의 설명에는 수학적 지식을 요하는 부분도 있을 것이지만 반드시 이해하라고 권하지는 않는다.
컴퓨터 프로그래밍을 해 보지 않은 사람이 이해하기는 당연히 어려울 것이다. 여기서는 프로그램 소스 코드로 설명을 하지 않고 알고리즘으로만 설명을 한다. 따라서 프로그래밍을 많이 해 본 사람도 이해하기 어려울 지 모른다. 이는 필자의 논리력이 미숙한 것이므로 나중에 필요에 따라 C 혹은 Basic 관련 언어로 부분에 대한 프로그램 소스를 추가할 수도 있다. 또한 되도록 쉽게 이해시키는 위주로 알고리즘을 소개할 것이므로 이보다 더 빠른 계산을 위한 알고리즘은 얼마든지 있음을 생각할 수 있다.
프로그램을 알고리즘의 조합으로 표현할 수 있다. 그럼 실제로 원주율을 어떻게 계산하는지 수학적, 알고리즘적으로 알아보자.

3. 수학적 분석

arctan 함수를 이용한 pi 계산

탄젠트 함수(tan x)에 대한 역함수는 아크탄젠트 함수(arctan x)이다.
예를 들어 tan x = y ( -pi/2 < x < pi/2, x는 라디안(radian) 단위 )의 관계가 있을 때 arctan(y) = x 가 된다.
Maclaurin 급수에 의해 arctan 식은 다음과 같다(증명은 생략한다).

식.1 (식. 1)

이 때 x값이 0에 가까와질 때 이 식에 의해 매우 빨리 수렴함을 알 수 있다.
또 우리는 pi/4 = arctan(1) 임을 알고 있다. 따라서 아래식
식.2 (식. 2)

으로도 pi를 계산할 수 있음을 알 수 있다. 하지만 위 식은 수렴성이 매우 낮다. 아래의 식은 위 식과는 비교도 안 될 정도로 빠름을 알 수 있다.
식.3 (식. 3)

arctan(1) = arctan(1/2) + arctan(1/3) 임은 쉽게 증명될 수 있지만 여기선 생략한다. 결국 pi = 4 arctan(1) = 4 arctan(1/2) + 4 arctan(1/3) 와 같은 식으로 구할 수 있으며, m arctan(1/n) 을 계산하는 함수를 만들었을 때, n의 값이 클 때 빨리 수렴한다. 필자가 자주 사용했던 식은 다음과 같다.
식.4 (식. 4)

위 식에 대한 증명은 간단하지만 여기선 생략하며 더 빠른 수렴을 보이는 다른 식도 많이 있지만 여기선 생략한다.

4. 알고리즘 연구

여러자리 수의 표현

의문 제기 - 100자리의 정수를 나타내려면? 100자리의 정수의 덧셈이나 뺄셈은 어떻게 할까? 소수 100자리까지 정확하게 수를 나타내려면? 등등의 문제를 해결하기 위해서는 그냥 integer, double 형의 단일 변수로는 유효자리수의 한계를 느낀다. 이를 해결하기 위해서 어떤 방법이 있을까?

방법 소개- 100자리의 수를 표현하는 경우를 예로 들 때, 배열(혹은 포인터 등)에 다음과 같이 수를 넣는 방법을 생각해 보자. a(0) = 1(10의 0제곱)의 자리의 수, a(1) = 10의 자리의 수, a(2) = 100의 자리의 수, ..., a(99) = 10의 99제곱의 자리의 수. 이와 같이 하면 100자리의 수에 대한 덧셈이나 뺄셈도 쉽게 계산하고 표현할 수 있다. 계산을 더 빠르게 하기 위해서는 배열 a()에 각 두자리 수(혹은 그 이상 허용하는 범위까지)를 넣는 방법을 생각할 수 있다. 여기선 기본적으로 한자리 수만 넣는 것으로 설명하려 한다. 프로그래밍 시에는 필요에 따라 허용하는 범위까지 확장시킬 수 있음을 알자.

무한소수의 유한자리 표현

예1) 1 / 2 = 0.5 = 0.50000 ... 0 (유한소수)
예2) 28 / 3 = 9.11111 ... = 9.11111 ... 1 (순환마디 한 자리의 무한소수)
예3) 1 / 97 = 0.01030 ... = 0.01030 ... ? (순환마디 96 자리의 무한소수)

'예1'과 같이 유한소수의 경우는 굳이 배열을 사용하지 않고도 간단히 계산될 수 있다. 하지만 많은 자리수를 차지하는 유한소수의 경우 컴퓨터가 허용하는 범위를 벗어나면 정확하게 나타낼 수 없을 것이다. '예2'와 '예3'에서처럼 무한소수는 일일이 배열에 넣어 줄 필요를 느낄 것이다. 물론 순환마디를 알 경우 배열에 넣는 쉬운 방법을 생각할 수 있을 것이다.

어쨌든 무한소수를 소수 몇 자리까지 나타낸다는 것은 결국은 근사값을 나타내는 것이 되며 이 때 그 이하자리는 반올림, 올림, 버림의 세 가지 중 한 가지를 선택할 수 있다. 이 때, 오차의 한계를 알 수 있으며, 알려 줘야 정확성을 알 수 있다. 예를 들어 소수 100자리까지 표현할 경우 오차는 소수100자리의 수 1보다 작을 것이며, 그러한 수 10개를 덧셈할 경우 오차는 100자리의 수 1의 10배보다 작을 것이다.

앞으로 원주율 계산을 위해 사용될 함수에 사용될 수들은 정수부분에는 크게 제한을 받지 않는다. 따라서 어떤 수를 소수 n자리까지 나타내는 배열(a 배열을 예로 들자)을 다음의 예와 같이 표현하기로 하자.

예) 348.2378948789773492 를 소수 10자리까지 a 배열에 넣는다.
a(0) = 348, a(1) = 2, a(2) = 3, ... , a(9) = 8, a(10) = 9.
소수 11자리부터는 버림을 한다.
이 때 오차(여기서 및 앞으로 이 글에서 오차는 오차의 절대값을 말함.)는 10^-10 (10의 -10 거듭제곱. '^'는 앞으로도 거듭제곱을 의미함.)보다 작다.

어떤 수의 정수부분을 a(0)에 넣고, 소수 k번째 수를 a(k)에 넣으면 어떤 수는
a(0).a(1)a(2) ... a(k) ... a(n) 과 같은 형태가 된다. a(0)만 여러 자리수가 될 수 있으며(물론 허용범위(예를 들어 integer 형은 32767)를 넘어서지 않는 경우에 한하여), a(1)에서 a(n)까지는 0~9 사이의 한 자리의 수가 된다.

앞으로는 이렇게 표현되는 수에 대한 연산을 다룰 것이며 수 a(0).a(1)a(2) ... a(k) ... a(n) 를 그냥 a()라고 표현할 것이다. 그러니까 수를 배열로 나타낸다고 생각하면 된다.

유한자리 소수의 덧셈, 뺄셈, 나눗셈

a와 b의 덧셈을 하는 경우를 생각해 보자. 그 더한 결과를 a나 b에 저장시킬 수도 있고, 또 다른 수 c에 저장할 수도 있다. 각기 장단점이 있을 수 있으며 호출하는 방법에 따라 결국 같은 결과를 만들 수 있으므로 여기서는 편의상 a+b=c의 연산을 설명한다(참고로 나중에 설명할 m arctan (1/n) 알고리즘시 편리함을 느낄 것이다).

더하기 알고리즘

  1. c()를 0으로 초기화시킨다.
  2. c(n)에 a(n)과 b(n)의 합을 넣는다(이것을 c(n)=a(n)+b(n)로 표현하자).
  3. c(n)이 9보다 크면 c(n-1)=1, c(n)=c(n)-10을 행한다.
  4. k = n-1 (k에 n-1을 넣는다.)
  5. c(k) = c(k) + a(k) + b(k)
  6. c(k)의 값이 9보다 크면 c(k-1)=1, c(k)=c(k)-10을 행한다.
  7. k값이 n-2, n-3, ..., 2, 1 에 대해 각각 5번과 6번을 반복한다.
  8. c(0) = c(0) + a(0) + b(0)
원리는 간단하다. 끝자리부터 한 자리씩 더해 나가는데 더한 값이 9를 넘을 경우는 앞자리로 1을 올려 주고 그자리는 더한 값에서 10을 빼는 것이다. 이런 식으로 계속, 계산할 때 carry(올려준 값이 있으면 1, 없으면 0)를 더해 준다. c(0)는 정수부분이므로 한자리의 수가 아닐 수도 있지만 c(1)부터 c(n)까지는 모두 0에서 9까지의 한자리수이다.
이렇게 하면 a()와 b()의 합이 c()에 저장된다. 여기서 주의해야할 것은 c(0)의 값이 한계값을 넘지 않아야 한다.
만일 a()와 b()가 각각 근사값이었을 경우 c()의 오차의 한계는 2*10^-n(10의 (-n) 거듭제곱의 두 배)이다.

빼기 알고리즘

  1. c()를 0으로 초기화시킨다.
  2. c(n) = a(n) - b(n)
  3. c(n)이 0보다 작으면 c(n-1)=-1, c(n)=c(n)+10 을 한다.
  4. k = n-1
  5. c(k) = c(k) + a(k) - b(k)
  6. c(k)의 값이 0보다 작으면 c(k-1)=-1, c(k)=c(k)+10을 한다.
  7. k값이 n-2, n-3, ..., 2, 1 에 대해 각각 5번과 6번을 반복한다.
  8. c(0) = c(0) + a(0) - b(0)
원리는 간단하다. 끝자리부터 한 자리씩 빼 나가는데 뺀 값이 음수일 경우는 앞자리를 1 내려 주고 그자리는 뺀 값에서 10을 더한다. 이런 식으로 계속, 계산할 때 carry(내려준 값이 있으면 -1, 없으면 0)를 더해 준다. c(0)는 정수부분이므로 한자리의 수가 아닐 수도 있지만 c(1)부터 c(n)까지는 모두 0에서 9까지의 한자리수이다.
이렇게 하면 a()에서 b()를 뺀 값이 c()에 저장된다. 오차의 한계는 2*10^-n이다(실제로 이보다 더 작게 잡을 수도 있지만 편의상 이렇게 잡는다).

나누기 알고리즘

  1. c(0) = a(0)/p의 정수부분
  2. tmp = a(0) - p * c(0) ; 나머지 = 피제수 - 제수 * 몫
  3. k=1
  4. tmp = 10 * tmp + a(k)
  5. c(k) = tmp/p의 정수부분
  6. tmp = tmp - p * c(k)
  7. k = 2, 3, ..., n-2, n-1 에 대해서도 4번,5번,6번을 반복한다.
  8. tmp = 10 * tmp + a(n)
  9. c(n) = tmp/p의 정수부분
원리는 간단하다. '피제수 = 몫 * 제수 + 나머지'를 반복하는데, 처음의 피제수는 a(0)가 되며, 다음은 그 나머지를 10배한 수에 a(1)을 더한 수, 이렇게 계속적으로 적용시켜 셈하는 것이다. 여기서 1,5,9번에 나온 나눗셈의 정수부분을 구하는 것은 내장함수를 사용할 수도 있고, 새로운 함수를 만들어 사용할 수도 있다. 또 2,6번 등은 나머지 구하는 함수를 사용해도 되며 이를 이용하여 1번과 2번, 5번과 6번 등을 순서를 바꿔 나머지부터 구한 뒤 정수부분은 적당히 구하도록 할 수도 있다. 오차의 한계는 10^-n으로 잡자. 나누기는 여러번 연속으로 해도 오차의 한계는 여전히 10^-n으로 잡을 수 있다.
여기서 주의해야 할 것이 있다. 제수 p의 값이 유효 범위의 10분의 1보다 클 경우 p와 tmp가 같은 형일 경우 버그가 있을 수 있다. 예를 들어 p를 integer 형으로 했을 때, integer 형의 허용한계가 32767일 때, p 값이 3276보다 커서는 안 된다. tmp = 10 * tmp + a(k) 를 하는 과정에서 over flow가 되어 잘못된 결과가 나올 수 있기 때문이다. 그러한 문제가 생길 수 있기 때문에 함수(프로시저)에서 p가 한계값을 넘어설 경우 에러처리를 하기를 권한다.

m arctan(1/n) 알고리즘

이제 위 알고리즘들을 이용해 m arctan(1/n) 을 계산하는 방법을 생각해 보자. 여기서 arctan 급수식을 다시 한 번 적어 보면 다음과 같다.

식.5 (식. 5)

m과 n은 여기서 양의 정수이다. 항 단위로 계산을 해 나갈 것이다. for 루프 등을 이용해 계산한다. 계산을 편리하게 하기 위해선 다음과 같은 세 가지의 수(a, b, s 의 배열)를 이용할 수 있다. 그럼, 그 방법을 소개한다.

처음 루프(i=1)

두번째 루프(i=2)
세번째 루프(i=3)
... ( 중 략 ) ...

k번째 루프(i=k)

수렴 및 오차

이런 식으로 a, b, s에 수를 기억시키면 결국 m arctan 1/n 의 계산 결과는 s()의 값과 같게 된다. 이 때 b()=0.00000...0 이 되는 순간, s() = m arctan 1/n 는 수렴한다. 몇 번째 루프에서 수렴할 것인지는 b에 관한 식으로부터 구할 수 있다. 이 때 루프를 두 번 도는 동안 오차의 한계를 마지막 자리 1로 잡을 수 있다. 그러므로 k번째 루프에서 수렴할 경우 오차의 한계는 {(k+1)/2 의 정수값 * 마지막 자리}로 잡을 수 있으며, 참값은 '근사값 - 오차의 한계' ~ '근사값 + 오차의 한계'가 된다.
주의 - 루프 변수 i가 한계값(예를 들어 max_int)을 넘어서지 않도록 해야 한다. 이것은 소수 몇번째 자리까지 구할 것인지, 그러면 몇번째 루프에서 수렴할 것인지 등을 짐작해서 처리하도록 해야 할 것이다.

이상으로 유한자리의 소수끼리의 덧셈과 뺄셈, 유한자리의 소수를 자연수로 나누는 것, 이를 이용해 m arctan 1/n 을 계산하는 것을 알아 보았다. 또 이에 대한 오차의 한계를 앎으로써 어디까지는 참값과 정확히 일치하는지 알 수 있다. 또 여러 개의 m arctan 1/n 에 대한 덧셈과 뺄셈을 할 수 있으며 이에 대한 오차의 한계는 각각을 더함으로써 구할 수 있다. 이제 실제 pi를 구하는 식들에 대해 알아 보자.

5. 원주율 계산 식

여기까지 성실히 읽으신 분은 아래에 소개되는 식으로 pi를 계산할 수 있으리라 생각된다(물론 프로그래밍 경험이 있는 분이라면). 아래의 식들에 대한 증명은 각자 해 보기 바란다.

pi/4 = arctan(1) 따라서,
pi = 4 arctan(1/1)

위 식은 수렴성이 매우 낮다. 권하고 싶지 않은 식이다. 소수 다섯 자리까지 계산하는 데도 엄청난 루프를 돌려야 할 것이다. 앞으로는 pi/4 에 관한 식만 소개한다. pi 식은 이 식의 양변을 4배하여 생각할 수 있다.

pi/4 = arctan(1/2) + arctan(1/3)

수렴 속도가 처음 식과는 비교할 수 없을 정도로 향상되었다. 하지만 더 좋은 식이 요구된다.

pi/4 = 2 arctan(1/3) + arctan(1/7)

pi/4 = 4 arctan(1/5) - arctan(1/239)

pi/4 = 8 arctan(1/10) - arctan(1/239) - 4 arctan(1/515)

이상으로 몇 가지 식을 소개했다. 이러한 식은 얼마든지 만들 수 있다. 여러분도 다른 좋은 식을 만들어서 더 빨리 수렴하도록 노력할 때 즐거움을 느낄 것이다. 이상으로 짧은 글(?)을 마칠까 한다.

6. 프로그래밍

일단, 아래 링크를 보기 바란다. 풀이는 나중에 시간이 되면...

원주율 계산 java 소스 : http://mathman.kr/bbs/read.php?category=6&num=1
원주율 소수 1만자리까지 계산한 결과 : http://mathman.kr/bbs/read.php?category=6&page=1&num=4

최초 작성일 : 1998.3.