제곱근 역수와 마법의 수 0x5f3759df

해커뉴스 스레드[1]에 재미있는 글[2]이 소개되어 있다. 요즘 포스팅 거리가 해커뉴스밖에 없는 듯-_- ㅋㅋㅋㅋ

3d 그래픽과 관련된 프로그래밍을 할 때, 공간내에서 빛의 반사 경로를 계산하려면 곡면 위 각 점에 대해 법선 단위 벡터(곡면에 수직이고 길이가 1인 벡터)를 계산할 필요가 있다. 각각의 점에 붙은 벡터 (x, y, z)에 대해 단위 벡터는 다음과 같이 계산된다.

\displaystyle \frac{1}{\sqrt{x^2 + y^2 + z^2}}

여기서 x^2 + y^2 +z^2은 비교적 빠르다. 결정적인 부분은 제곱근의 역수인데, 공간상의 수많은 점에 대해 이 계산을 하는 것이 비교적 고사양을 요구하는 녹록치 않은 계산이다. 이 계산의 속도가 빨라지면 게임의 처리 속도가 올라갈 수 있다.

마술과도 같은 다음 코드를 보자. ㅋㅋ

float FastInvSqrt(float x) {
  float xhalf = 0.5f * x;
  int i = *(int*)&x;         // evil floating point bit level hacking
  i = 0x5f3759df - (i >> 1);  // what the fuck?
  x = *(float*)&i;
  x = x*(1.5f-(xhalf*x*x));
  return x;
}

이 코드는 근사적이지만 매우 빠른 속도로 제곱근의 역수를 계산해주는 c코드이다. 이건 존 카맥이 Quake III Arena에 써서 유명해진 코드이지만, 그 계산의 기원은 뿌리가 깊다고 한다. 위키피디아에도 이 함수에 관한 항목이 있으니 Fast inverse square root 항목을 읽어보시라. ㅎ

이게 왜 작동을 할까? 특히 4번째 줄의 괴악한 수 0x5f3759df는 또 뭐란 말인가? 주석대로 진짜 what the fuck?이다-_- 위키피디아에 따르면 퀘이크 소스코드에 있는 주석이라고 한다. ㅋ

어느 사람의 블로그 포스트[2]에 상세한 설명과 더불어 일반화하는 방법까지 있는데, 본 포스트에서는 대략적인 설명과 더불어 그 블로그 포스트[2]에서 설명을 안 하고 넘어간 6번째 라인의 계산에 대해 잠시 설명해볼까 한다.

위 함수의 3번째 줄의 코드는 floating point를 그대로 integer 로 바꾸는 줄이다. 이 자체는 별거 없지만 앞으로 비트단위의 연산을 하게 되리라는 예측이 가능하다. 한 번 따라가보자.

먼저 32비트 single-precision floating-point format에 대해 알 필요가 있는데, 맨 첫 비트가 부호를 나타내고, 다음 8개의 비트는 지수(exponent), 나머지 23개의 비트는 유효숫자(significand)를 의미한다. (로그의 fractional part를 의미하는 mantissa와 약간 다름) 저 블로그[2]의 notation대로 이 값을 e, m이라 하고, 이 값들이 정수화된 값을 E, M이라 하면 다음과 같은 관계가 성립한다.

\displaystyle m = \frac{M}{2^{23}}, \displaystyle e = E- 127 … (1)

최초 주어진 floating point는 (1+m)2^e이고 이것이 3번째 줄에서 integer type으로 바뀌면 M+2^{23}E가 된다. 우리의 계산 목적은 주어진 값 x = (1+m_x )2^{e_x}에 대해 y = x^{-\frac{1}{2}} = (1+m_y )2^{e_y}를 계산하는 것이다. 양변에 로그를 취하면

\displaystyle \log_2 (1+m_y )+e_y = - \frac{1}{2} (\log_2 (1+m_x )+e_x ) … (2)

물론 로그함수도 녹록치 않은 계산이지만, 로그는 쉽게 선형근사 가능하다. 만약 적절한 상수 \sigma를 선택한다면 다음과 같이 근사가능하다.

\displaystyle \log_2 (1+v) \approx v+\sigma

따라서 이 결과를 (2)에 적용하면

\displaystyle m_y +\sigma +e_y \approx - \frac{1}{2} (m_x +\sigma+e_x )

여기서 (1)을 대입하여 정리하면

\displaystyle \begin{aligned} \frac{M_y}{2^{23}} + \sigma + E_y - 127 & \approx - \frac{1}{2}\left( \frac{M_x}{2^{23}} + \sigma + E_x - 127 \right) \\ \frac{M_y}{2^{23}} + E_y & \approx - \frac{1}{2}\left( \frac{M_x}{2^{23}} + E_x \right) + \frac{3}{2}(127 -\sigma) \\ M_y +2^{23}E_y & \approx -\frac{1}{2}(M_x + 2^{23}E_x) + 3\cdot 2^{22}(127 - \sigma) \end{aligned}

마지막 줄의 좌변이 바로 우리가 원하는 값이 된다. 이 값을 floating point로 변환하는 과정이 코드의 5번째 줄이 된다. 다시 말해 마지막 줄의 식은 최초 주어진 값의 정수화된 값 I_x와 목표로 하는 값 I_y 사이에 다음 관계식이 있음을 의미한다.

\displaystyle I_y \approx 3\cdot 2^{22}(127 - \sigma) - \frac{1}{2}I_x

2로 나누는 작업은 1비트씩 옆으로 움직여서 쉽게 얻을 수 있고, 위 오리지널 소스코드에서는 \sigma의 값을 0.0450465로 잡았다고 한다. 그러므로

\displaystyle 3\cdot 2^{22}(127 - \sigma) = 1597463007.854592

가 되고, 정수 부분이 십육진수로 마법의 수 0x5f3759df가 된다. 바로 소스코드의 4번째 줄에 해당한다. ㅎㅎㅎ

블로그[2]에서는 6번째 줄의 계산은 Newton method라고만 나와있고 별다른 설명이 없는데, 어째서 저런 계산이 되는지에 대해 사족을 덧붙여 보자. ㅋㅋ

뭐 Newton method는 일단 한 번 찾은 근의 근사치의 정밀도를 올리는 방법인 줄은 아실 터. 주어진 값 p에 대해 참값 t = 1/\sqrt{p} 와 근사값 t_1이 있다면 t의 방정식 1/t^2 =p을 푸는 셈이므로 f(t)= 1/t^2 -p라 두면

\displaystyle t_2 = t_1 - \frac{f(t_1)}{f'(t_1)} = t_1 + \frac{1}{2}\left(t_1 - pt_1^3\right) = t_1 \left(\frac{3}{2} - \frac{1}{2}pt_1^2 \right)

가 되는 것이다. 이것이 여섯 번째 줄에 해당한다. 이 근사법을 한 번 더 쓰면 정밀도가 더 올라간다.

이 코드는 위키피디아에 따르면 존 카맥은 Terje Mathisen이라는 어셈블리 프로그래머에게 제안을 받았다고 한다. Mathisen도 완전히 창작한 것은 아닌 듯 하고, 이 코드의 기원은 3D 그래픽의 기원까지 거슬러 올라간다고 하니, 웬만한 코드 구루들은 아마 이미 알고 있었을 것 같다. 저 위의 매직 넘버 0x5f3759df를 어떻게 결정했는지는 불명확한 듯. 64비트용 매직 넘버로 0x5fe6eb50c7b537a9가 있다고 한다. ㅎ

 


곰곰히 생각해 봤는데 애초에 floating point format이 주어진 값을 어느정도 로그화한 형태로 저장하는 방식이라 (지표 + 가수), 이것을 비트채로 다루면 거듭제곱이 어떤 선형식의 계산과 유사해지는 현상을 이용하는 테크닉이 아닌가 싶다.

 


[1] https://news.ycombinator.com/item?id=8519365
[2] 0x5f3759df in Hummus and Magnets

답글 남기기

아래 항목을 채우거나 오른쪽 아이콘 중 하나를 클릭하여 로그 인 하세요:

WordPress.com 로고

WordPress.com의 계정을 사용하여 댓글을 남깁니다. 로그아웃 / 변경 )

Twitter 사진

Twitter의 계정을 사용하여 댓글을 남깁니다. 로그아웃 / 변경 )

Facebook 사진

Facebook의 계정을 사용하여 댓글을 남깁니다. 로그아웃 / 변경 )

Google+ photo

Google+의 계정을 사용하여 댓글을 남깁니다. 로그아웃 / 변경 )

%s에 연결하는 중