Geometry & Recognition

정수의 제곱근을 정수 연산만 가지고 계산해야 하는 경우가 많이 있는데 (특히, 임베디드 환경에서 프로그래밍을 하는 경우) 이를 구현하는 방법에 대해서 알아보자. 32비트 정수의 제곱근은 16비트 이하이므로 가장 최상위의 비트 값이 세팅이 된 정수를 가지고 제곱근 조건을 만족시키는지 확인하고, 만족시키면 다음의 하위 비트 값이 세팅이 된 정수 값을 기존의 값에 더해서 다시 제곱근 조건을 만족시키는지 확인하는 과정을  마지막 비트까지 반복하면 된다. 그런데 가능한 가장 최상위의 비트 값을 갖는 정수는 0x8000 (= 1<<15) 이므로 양의 정수 x의 제곱근을 구하는 과정은 아래와 같다.

           unsigned s = 0x8000;  //= (1 << 15);
           unsigned t = 0;
           do {                          //새로 세팅된 값이 제곱근 조건을 만족하는가?
                if ((t + s) * (t + s) <= x ) t = t + s;  //만족하면 하위 비트값을 더함;
                s = s >> 1;                                //다음 하위 비트로 내려감;
           } while (s != 0) ;
           return (t);

이 알고리즘은 16번의 반복으로 항상 끝이 나지만, 중간에 곱셈 연산이 들어가 있다. 이 곱셈 연산을 없애는 방법을 알아보자. 먼저           

 (t + s) * (t + s) = t * t + 2 * t * s + s * s;

이므로 if 문 ( ) 내부 식 양변에서 t^2를 빼도록 하자. 이 경우 t의 값이 갱신이 되면

 x - t^2 => x - (told + sold)^2 => x - told * told - (2 * told * sold + sold * sold) ;

이므로, x 값에서도 역시 2 * t * s + s * s의 값을 제거해야 한다. 따라서

                  unsigned nr = 2 * t * s + s * s;
                  if (nr <= x) {
                         x = x - nr;
                         t = t + s;
                  }

로 고쳐야 한다. 이제, s * s 와 2 * t * s을 별도의 변수를 도입하여 컨트롤하면,

          unsigned m = s * s;    // m = 0x40000000;
          unsigned r = 2 * t * s ;// => 2 * (told + sold) * sold
                                        // = 2 * told * sold + 2sold * sold 
                                        // = r + m + m = nr + m;
          do {
                  nr = m + r ;
                  if (nr <= x) {
                         x = x - nr ;
                         r = nr + m;
                         t = t + s;
                  } ;
                  s = s >> 1;
                  m = m >> 2;
                  r = r >> 1;
          } while (s != 0);

로 사용이 가능하다. 그런데 여기서 r = 2 * t * s인데, 마지막으로 루프를 빠져나올 때 s가 1인 상태에서 s / 2를 하므로 (실수 연산을 고려하면) 이때 r 값은 t 값과 같다. 따라서 t 변수를 사용할 필요가 없고, 마찬가지로 s 변수도 사용할 필요가 없어진다 (m을 컨트롤 변수로 이용). 

//16-step integer_sqrt;
unsigned isqrt( unsigned x ) {
       unsigned m = 0x40000000 ;   //=(1 << 30)
       unsigned r = 0, nr;
       do {
            nr = m + r ;
            if (nr <= x) {
                   x -= nr ;
                   r = nr + m;
            }
            m >>= 2;
            r >>= 1;
       } while (m != 0);
       return (r) ;
};

결과는 항상 <= sqrt(x) 이다. 실제로 빠를까?

한번 더 다듬도록 하자. 나중에 r을 시프트 하지 말고 먼저 하면, r = 2 * t * s => r / 2 = t * s 이므로, if 문의 r은 r =  r + m이 된다. 정리하면

unsigned isqrt( unsigned x ) {
       unsigned m = 0 x 40000000 ;   //= (1 << 30)
       unsigned r = 0, nr;
       do {
            nr = m + r ;
            r >>= 1;
            if (nr <= x) {
                   x -= nr ;
                   r += m;
            }
            m >>= 2;
       } while (m != 0);
       return (r) ;
};

이 된다. 이는 잘 알려진 integer sqrt의 값을 구하는 함수이다.

다시 한번 더 형태를 바꾸면, r = r + m 항을 보면, r은 m을 누적하는 변수이고, m은 하나의 비트만 세팅이 된 변수이므로 덧셈을 비트 or 연산으로 바꿀수 있다.

             r += m   ==> r |= m;

그리고 매크로를 도입하여 루프를 풀어쓰기 하면, 아래의 코드를 얻을 수 있다.

unsigned isqrt( unsigned x ) {
       unsigned m = 0x40000000 ;   //= (1 << 30)
       unsigned r = 0, nr;
#define STEP(k)             \
           nr = m + r;      \
           r <<= 1;         \
           if (nr <= x) {   \
                x -= nr;    \
                r |= m;     \
           }                 
       STEP(15) ;  m>>=2;
       STEP(14) ;  m>>=2;
       STEP(13) ;  m>>=2;
       STEP(12) ;  m>>=2;
       STEP(11) ;  m>>=2;
       STEP(10) ;  m>>=2;
       STEP(9)  ;  m>>=2;
       STEP(8)  ;  m>>=2;
       STEP(7)  ;  m>>=2;
       STEP(6)  ;  m>>=2;      
       STEP(5)  ;  m>>=2;
       STEP(4)  ;  m>>=2;
       STEP(3)  ;  m>>=2;
       STEP(2)  ;  m>>=2;
       STEP(1)  ;  m>>=2;
       STEP(0)  ;
       return (r) ;
};

이 코드 또한 인터넷 상에 ARM 개발자 사이트나 게임 개발자 사이트에 널리 알려진 코드이다.