정수 범위에서 제곱근의 최적화된 구현 방법은?

C/C++를 포함한 거의 대부분의 언어에서는 부동소수점의 제곱근 함수를 지원하고 있다.

그리고, 정수 범위에서 제곱근을 구할 때는 보통 이 제곱근 함수를 활용하여 간단히 구현한다.

 

그런데, 순수하게 정수 범위에서 동작하는 제곱근 함수가 있다면 조금은 더 빠르게 동작할 것 같다.

파이썬의 경우 아예 이러한 함수를 stdlib에서 지원한다.

 

정수 범위에서 제곱근을 구하는 아이디어에는 몇 가지 방식이 있다.

 

1. 뉴턴의 근사법 활용

?

1 2 3 4 5 6 7 8 9 10 11 12

unsigned ISqrtByNewton32(const unsigned n) {     unsigned x = n;     unsigned y = (x + 1) / 2;       while (y < x) {         x = y;         y = (x + n / x) / 2;     }       return x; }

 

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

def improved_i_sqrt(n):     assert n >= 0     if n == 0:         return 0     i = n.bit_length() >> 1    # i = floor( (1 + floor(log_2(n))) / 2 )     m = 1 << i    # m = 2^i     #     # Fact: (2^(i + 1))^2 > n, so m has at least as many bits     # as the floor of the square root of n.     #     # Proof: (2^(i+1))^2 = 2^(2i + 2) >= 2^(floor(log_2(n)) + 2)     # >= 2^(ceil(log_2(n) + 1) >= 2^(log_2(n) + 1) > 2^(log_2(n)) = n. QED.     #     while (m << i) > n: # (m<<i) = m*(2^i) = m*m         m >>= 1         i -= 1     d = n - (m << i) # d = n-m^2     for k in xrange(i-1, -1, -1):         j = 1 << k         new_diff = d - (((m<<1) | j) << k) # n-(m+2^k)^2 = n-m^2-2*m*2^k-2^(2k)         if new_diff >= 0:             d = new_diff             m |= j     return m

 

3. 이진 탐색 알고리즘 활용

미리 제곱수를 계산해놓은 뒤 이진 탐색 아이디어를 적용하는 방법도 있다. (관련 사이트)

 

4. Intrinsic 함수 활용

SSE가 탑재된 CPU에서는 _mm_sqrt_ss() 함수를 활용해서 더 빠른 연산을 수행할 수 있다.

float에서 처리되지만, 엄청난 성능을 보여준다.

?

1 2 3 4 5 6 7

int ISqrtByIntrinsic32(const int n) noexcept {     register __m128 xmm0 = _mm_setzero_ps();     xmm0 = _mm_cvt_si2ss(xmm0, n);     xmm0 = _mm_sqrt_ss(xmm0);     return _mm_cvtt_ss2si(xmm0); }

 

이 외에도 중학교에선가 수학시간에 배웠던 손으로 계산하는 방법을 적용하는 알고리즘도 있다. (관련 사이트)

 

---

 

여러가지 정수 제곱근 알고리즘을 int16_t 및 int32_t 용으로 구현해서 100만번씩 실행해 성능을 비교해봤다.

 

 

놀랍게도 그냥 sqrt() 함수를 실행시켜 정수형으로 형변환한 것이 거의 가장 빨랐다.

Intrinsic 함수를 사용하는 쪽이 조금 더 빠르긴 했지만, 굳이 별도의 알고리즘을 적용할 필요는 사실상 없는 것 같다.

 

?

1 2 3 4 5 6 7 8 9 10 11 12 13

unsigned short 99969066 times CISqrt::ISqrtBySqrt16(): 0.11 CISqrt::ISqrtByFastSqrt16(): 0.109 CISqrt::ISqrtByNewton16(): 1.391 CISqrt::ISqrtByBinSearch16(): 0.25 CISqrt::ISqrtByIntrinsic32(): 0.109   unsigned 100000000 times CISqrt::ISqrtBySqrt32(): 0.204 CISqrt::ISqrtByNewton32(): 2.015 CISqrt::ISqrtByBinSearch32(): 0.64 CISqrt::ISqrtByIntrinsic32(): 0.109 CISqrt::ISqrtByMultiplyImproved(): 1.719

 

이게 뭔 일일가 싶어 디스어셈블 해보니 이유를 할 수 있었다.

컴파일러가 알아서 AVX 명령어로 컴파일하기 때문^[각주:1]이었다.

?

1 2 3 4 5 6 7 8

00007FF6FA381AA0  mov         ecx,ebx  00007FF6FA381AA2  vxorps      xmm1,xmm1,xmm1  00007FF6FA381AA6  vcvtsi2sd   xmm1,xmm1,rcx  00007FF6FA381AAB  vxorpd      xmm0,xmm0,xmm0  00007FF6FA381AAF  vucomisd    xmm0,xmm1  00007FF6FA381AB3  ja          main+93Bh (07FF6FA381ABBh) 00007FF6FA381AB5  vsqrtsd     xmm0,xmm1,xmm1  00007FF6FA381AB9  jmp         main+944h (07FF6FA381AC4h)

 

한줄 요약

 

복잡한 생각 필요 없고 그냥 sqrt() 사용하자

 

1. 아래 화면은 컴파일러 설정을 AVX2로 한 결과이며, SSE2로 설정하면 거기 맞는 코드가 생성됨 [본문으로]