BigFloat로 Pi를 구해보자-9. 제곱근 구현하기

최적화하는건 머리가 너무 아프므로, 이번엔 산뜻하게 제곱근 계산하는 함수나 만들어보렵니다. 사칙연산을 다 구현하면 BigFloat 자료형을 다 만들었다고 착각하면 아니되옵니다.ㅋ 이제 시작이지요. 그중에서 제일 기본적인게 제곱근입니다. 계산도 어렵지 않아요.

나눗셈 구현할때 사용했던 Newton-Raphson 방법을 기억하시려는지요? 이를 제곱근 계산에도 이용할 수 있습니다.(바빌론에서 제일 처음 사용되었다고 해서 바빌로니아 방법이라고 부르기도 합니다.)

이 공식 기억나지요? f(x)=0의 근을 근사해주는 Newton-Raphson의 공식이었죠. 우리는 제곱근을 구하고 싶으므로 f(x)를 다음과 같이 놓겠습니다.

(f(x)의 근은 sqrt(S) (혹은 -sqrt(S) )가 되겠지요?)

그럼 아래와 같이 정리가 될겁니다.

나눗셈과 덧셈, 시프트 연산만 가지고 충분히 제곱근을 계산할 수 있겠네요!

Newton-Raphson 공식에서 빠른 수렴을 위해서는 초기값이 좋아야겠지요. 초기값 X0은 어떻게 추정하면 좋을까요?

임의의 양수 S는 다음과 같이 나타낼수 있습니다. 당연히 n은 정수입니다.

정수 n은 비트수를 세는것 만으로 구할수 있습니다. n이 짝수일 경우 m비트 정수는 sqrt(S)와 최대 1.414배밖에 차이나지 않을 겁니다. n이 홀수일 경우 m비트 정수는 sqrt(S)와 최대 2배밖에 차이나지 않습니다.

즉 S의 비트수를 세어 n을 구한뒤, 를 구하면 우리가 구하고자 하는 sqrt(S)와 최악의 경우에도 2배밖에 차이가 나지 않는다는 것이지요. 이정도면 빠르고 괜찮은 근사법입니다.

이제 이 사실들을 바탕으로 코드를 작성해봅시다.

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

void RR_sqrt(PRR dest, PRR a) {     if(a->sign) return; // 음수면 그냥 무시하는 걸로 합시다.     RR_initu(dest, 0, 0);     if(RR_iszero(a)) return; // 0이면 계산할 가치도 없죠     RRU32 w = i_bscs(a->frac, a->size); //1이 나타나는 최상위 비트를 찾아냅니다. //소수점을 31번과 32번 비트 사이로 하기로했었죠?     if(w < 32)     { //32보다 작다는것은 a->frac[0]에 존재한다는 겁니다. //즉 1보다 크거나 같다는 것.         dest->frac[0] = 1 << ((31 - w) / 2);     }     else     { //32이상 이라는 것은 a->frac[0]이 0이라는 것이고 //이는 1보다 작다는 것.         RRU32 shift = (w - 32) / 2;         w = shift >> 5;         shift &= 31;         if(w + 1 >= dest->size) return;         dest->frac[w + 1] = 0x80000000 >> shift;     }       PRR t = RR_create(dest->size + 1); //임시로 사용할 변수 생성       RRU32 dsize = dest->size;     RRU32 asize = a->size;           t->size = 4;     dest->size = 4;     a->size = 4; //처음엔 크기 4로 시작합니다.           RR_inv(t, dest); // t = 1/x     RR_muladd(dest, t, a); // nx = x + a/x     RR_shr(dest, 1); // 2로 나누는것과 같습니다.       RRU32 stop = 123456789; //임의의 값     while(1)     {         RR_inv(t, dest); // t = 1/x         RR_muladd(dest, t, a); // nx = x + a/x         RR_shr(dest, 1); //2로 나누는것과 같습니다.           if(dest->size == dsize)         {             if(stop == dest->frac[dsize - 2]) break; //제일 마지막 바로 전 자리가 변하지 않는다면, //이미 안정한 수준에 도달했다는것이므로 //계산을 중단합니다.             stop = dest->frac[dsize - 2];         }           t->size = MIN(t->size*3/2, t->alloced);         dest->size = MIN(dest->size*3/2, dsize);         a->size = MIN(a->size*3/2, asize); //크기를 1.5배씩 증가시켜요.     }       RR_delete(t); }

중간에 웃지 못할 코드가 있는데요, stop = 123456789; 여기입니다.

Newton-Raphson 방법은 한 번 반복할때마다 정밀도가 대략 2배씩 늘어납니다. 그런데 실제 구현에서는 여러가지 한계에 의해서 실제로 그렇지 않은 경우도 있어서, 우리가 원하는 정밀도까지 계산해내기 위해 몇번 반복해야할지를 가늠하는게 쉽지 않습니다.(ㅠㅠ)

그래서 frac배열의 제일 마지막 바로 전 자리를 지켜보면서, 이 값이 바뀌지 않는 지점에 다다르면 계산을 중단하는 방법을 사용했어요. (제일 마지막 자리는 오차때문에 매번 반복 때마다 계속 값이 바뀌게 됩니다.)

시간복잡도는 Newton-Raphson 방법으로 나눗셈을 구현했을때와 마찬가지 방법으로 유도하면 O(n²)이 나옵니다.