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와 같다는 알 수 있습니다.

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

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

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

코드로 구현해봅시다.

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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

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를 계산해보아요.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108

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계산법을 다뤄볼게요~