BigFloat로 Pi를 구해보자-11. Pi 구하기: Machin 급수

저번 연재에서 아크탄젠트 급수를 이용해서 Pi를 계산해보았는데, 머리를 조금 더 굴리면 훨씬 더 빠른 알고리즘을 만들어낼 수 있습니다.

가장 대표적인 공식은 다음과 같아요.

John Machin이 1706년에 발견한 급수인데, Machin의 이름을 따서 Pi를 계산해내는 일련의 아크탄젠트 공식들을 Machin-like 공식이라고 부릅니다.^[각주:1]

그냥 지나가면 재미없으니, 저 공식이 어떻게 유도되었는지 거슬러 올라가봅시다. (복소수를 이용한 증명을 소개하니, 이해하려면 복소수에 대한 기초적인 지식이 필요해요..)

...(식1)

복소수 지수를 다루는 데 가장 중요한 공식이지요. 여기서 앞의 e^r는 복소수의 절대값(크기)이라 하고, theta는 편각이라고 합니다. (복소평면에 그려보면 x축과 이루는 각이 theta, 길이가 e^r이 나오니까요.)

두 복소수를 곱하면, 복소수의 절대값은 곱해지고, 편각은 더해집니다.

자, 이제 만약에 다음 식을 가지고 생각해볼께요.

아크탄제트 함수를 위처럼 변형해봤습니다. (식1)의 실수부가 분모, 허수부가 분자라는 것을 기억하세요. 다른 말로 쓰면,

는 복소수 가 나타내는 각과 같은 값이라는 것이죠.

이므로, 이 값은 (1 + i)가 나타내는 값과 같습니다.

여기까지 장황하게 설명한 것은 아크탄젠트들의 합차를 복소수의 곱셈으로 표현할수 있음을 설명하기 위해서였습니다. 자, 이제 이 식을 복소수로 바꿔볼까요

->

계산을 해보니 (1+i)의 배수가 나옵니다. 편각을 결정짓는건 실수부와 복소수의 비라는 것을 기억하세요. 편각이 Pi/4와 같다는 알 수 있습니다.

수학자들은 비슷한 기교를 통해 변태..같은 공식들을 많이 만들어냈는데, 다음과 같아요,,,

(핡.. 수학자들은 다들 변태야... 더 끔찍한 식은 여기에서 마음껏 볼수 있습니다.)

모두 다 똑같이 아크탄젠트 급수를 이용한다는 점에서는 변함없지만, 인수가 작아졌기때문에 수렴속도는 훨씬 빨라졌지요.

코드로 구현해봅시다.

먼저 아크탄젠트 함수부터~

void RR_atan(PRR dest, PRR a) { //atan(a) = a - a^3/3 + a^5/5 - a^7/7 ... PRR t = RR_create(dest->size + 1); PRR square = RR_create(dest->size + 1); PRR u = RR_create(dest->size + 1); RRU32 i = 1; RR_move(dest, a); RR_mul(square, a, a); RR_mul(u, a, square); while(1) { RR_invu32(t, 2*i++ + 1); RR_mulsub(dest, t, u); RR_mul(t, u, square); PRR temp = t; t = u; u = temp; RR_invu32(t, 2*i++ + 1); RR_muladd(dest, t, u); RR_mul(t, u, square); temp = t; t = u; u = temp; if(RR_iszero(u)) break; } RR_delete(u); RR_delete(square); RR_delete(t); }

저번 연재에서 사용했던 것과 거의 똑같습니다.

이제 이를 이용해 Pi를 계산해보아요.

void RR_pi_machin1(PRR dest) { // pi/4 = 4atan(1/5) - atan(1/239) PRR t = RR_create(dest->size + 1); PRR u = RR_create(dest->size + 1); RR_invu32(u, 5); RR_atan(dest, u); RR_shl(dest, 4); RR_invu32(u, 239); RR_atan(t, u); RR_shl(t, 2); RR_sub(dest, t); RR_delete(u); RR_delete(t); } void RR_pi_machin2(PRR dest) { /* pi/4 = 12atan(1/49) + 32atan(1/57) - 5atan(1/239) + 12atan(1/110443) */ PRR t = RR_create(dest->size + 1); PRR u = RR_create(dest->size + 1); RRU32 n; RR r; r.size = 1; r.alloced = 1; r.sign = 0; r.frac = &n; RR_initu(dest, 0, 0); RR_invu32(u, 49); RR_atan(t, u); n = 12*4; RR_muladd(dest, t, &r); RR_invu32(u, 57); RR_atan(t, u); n = 32*4; RR_muladd(dest, t, &r); RR_invu32(u, 239); RR_atan(t, u); n = 5*4; RR_mulsub(dest, t, &r); RR_invu32(u, 110443); RR_atan(t, u); n = 12*4; RR_muladd(dest, t, &r); RR_delete(u); RR_delete(t); } void RR_pi_machin3(PRR dest) { /* pi/4 = 183atan(1/239) + 32atan(1/1023) - 68atan(1/5832) + 12atan(1/110443) -12atan(1/4841182) + 100atan(1/6826318)*/ PRR t = RR_create(dest->size + 1); PRR u = RR_create(dest->size + 1); RRU32 n; RR r; r.size = 1; r.alloced = 1; r.sign = 0; r.frac = &n; RR_initu(dest, 0, 0); RR_invu32(u, 239); RR_atan(t, u); n = 183*4; RR_muladd(dest, t, &r); RR_invu32(u, 1023); RR_atan(t, u); n = 32*4; RR_muladd(dest, t, &r); RR_invu32(u, 5832); RR_atan(t, u); n = 68*4; RR_mulsub(dest, t, &r); RR_invu32(u, 110443); RR_atan(t, u); n = 12*4; RR_muladd(dest, t, &r); RR_invu32(u, 4841182); RR_atan(t, u); n = 12*4; RR_mulsub(dest, t, &r); RR_invu32(u, 6826318); RR_atan(t, u); n = 100*4; RR_mulsub(dest, t, &r); RR_delete(u); RR_delete(t); }

속도를 한번 비교해볼까요? 길이를 512(512*32비트 정밀도)로 하고, Pi를 구해봤습니다.

pi_machin3함수가 처음구현했던 pi_atan보다 두 배 이상 빠르네요!

(정밀도를 늘리거나 줄여도, 이 비율은 거의 그대로 유지됩니다. 왜냐면 지금까지 구현한 atan, machin1, machin2, machin3 알고리즘은 모두 같은 시간복잡도를 가지고 있거든요.)

하지만 여전히 느립니다. 더 빠른 알고리즘이 없을까요? 다음 연재에서는 훨씬 더 빠른, 코사인 함수를 이용한 incremental한 pi계산법을 다뤄볼게요~