BigFloat로 Pi를 구해보자-13. 가우스-르장드르 알고리즘

저번 연재에서는 코사인 함수를 이용하여 Pi계산속도를 비약적으로 끌어올렸었죠. 하지만 그게 최종 보스가 아닙니다. 수학자들은 수학자들답게 Pi를 구하는 무시무시한 알고리즘을 개발해냈는데, 이름만 들어도 무서운 가우스-르장드르 알고리즘(Gauss-Legendre algorithm) 입니다.

위키에도 잘 설명되어있지만, 굳이 간단하게 소개하자면 다음과 같아요.

먼저 초기값을 다음과 같이 줍니다.

그리고 아래 점화식에 따라 다음 항들을 차례로 계산해줍니다.

그러면 n이 커짐에 따라 다음값은 Pi에 수렴합니다.

으아... 전혀 이해가 되지 않는 알고리즘입니다. 산술-기하 평균의 극한이 1종 타원 적분과 연관있다는 사실을 이용하는 알고리즘이라는데... (솔직히 잘 모르겠어요 어려워요ㅠ)

중요한건 이 알고리즘은 n을 키워감에 따라 정밀도가 2배씩 증가하는 좋은 성질을 가지고 있다는 것입니다. 즉 n자리의 정밀도를 얻기 위해서는 log n번의 반복만으로 충분하다는 것이지요. 게다가 알고리즘에 이용되는 연산은 사칙연산과 제곱근 밖에 없습니다. 우리는 제곱근을 O(n²)에 계산해 낼수 있으니, 1회 반복에 걸리는 시간을 O(n²)으로 맞출수 있지요. 그 결과 우리는 O(n²log n)시간에 n자리 정밀도의 Pi를 계산해 낼수 있습니다!

(저번 연재에서 cos함수의 시간복잡도가 O(n²log n)보다 컸다는 것을 생각해보면, 가우스-르장드르 알고리즘이 얼마나 빠른것인지 알수 있을겁니다.)

이제 코드를 짜봅시다.

void RR_pi_gl(PRR dest) { PRR a = RR_create(dest->size + 1); PRR a2 = RR_create(dest->size + 1); PRR b = RR_create(dest->size + 1); PRR bt = RR_create(dest->size + 1); PRR t = RR_create(dest->size + 1); //오차를 줄이기 위해 크기를 넉넉하게 잡아줍시다! RR_initu(a, 1, 0); RR_init(bt, .5); RR_sqrt(b, bt); RR_init(t, .25); int c = i_bsu(dest->size) + 1; //몇번 반복해야하는지 구합니다. for(int i = 0; i < c; ++i) { RR_mul(bt, a, b); //bt = a_n * b_n RR_move(a2, a); RR_add(a, b); RR_shr(a, 1); // a_{n+1} = (a_n + b_n)/2 RR_sqrt(b, bt); // b_{n+1} = sqrt(a_n * b_n) RR_sub(a2, a); // a2 = a_n - a_{n+1} RR_mul(bt, a2, a2); RR_shl(bt, i); // 2^n*(a_n - a_{n+1})^2 RR_sub(t, bt); // t_{n+1} = t_n - 2^n*(a_n - a_{n+1})^2 } RR_add(a, b); RR_shr(a, 1); RR_mul(a2, a, a); RR_inv(bt, t); RR_mul(dest, a2, bt); //dest = ((a_n + b_n)/2)^2 / t_n RR_delete(a); RR_delete(a2); RR_delete(b); RR_delete(bt); RR_delete(t); }

중간에 쓰인 i_bsu 함수는 최상위비트가 몇번째 비트에서 나타나는지 알려주는 함수입니다. 내부구현은 지저분하니, math.h의 log함수를 이용해서, (int)(log(dest->size)/M_LN2) 이렇게 작성할 수도 있습니다. 마지막에 +1은 혹시나 오차로 인해 정밀도가 충분히 도달 못할지도 모르니, 넉넉하게 한 번더 반복을 돌려주기 위함입니다.