고속 역 제곱근

위 사진에 나타난 1인칭 슈팅 게임 OpenArena의 조명 처리와 같은 연산에 입사각과 반사각을 계산할 때 고속 역 제곱근 알고리즘을 사용한다.

고속 역 제곱근(高速逆-根, fast inverse square root)은 때때로 Fast InvSqrt()나 16진수 0x5f3759df라고도 하는, IEEE 754 부동소수점 체계의 32비트 실수에 대한 제곱근의 역수를 계산하기 위한 알고리즘이다. 이 알고리즘은 1990년대 초에 실리콘 그래픽스에서 개발한 것으로 추정되며, 1999년에 퀘이크 3 아레나의 소스 코드에 쓰인 것이 가장 유명하다. 이 알고리즘에 관한 논의가 2000년에 중국어 개발자 포럼 CSDN에서 있었고, 공중에 알려진 것은 2002년 또는 2003년에 유즈넷과 같은 공개 포럼에서였다. 이 알고리즘은 컴퓨터 그래픽스에서 조명과 셰이딩을 위한 입사각과 반사각 계산에 사용되면서 대규모 부동소수점 연산의 계산 비용 문제 해소에 도움이 되었다.

이 알고리즘은 실수를 입력 받아 나중에 사용할 절반의 값을 따로 저장한 다음, 입력 받은 float 실수를 long 비트로 취급한다. 한 비트를 오른쪽 논리 시프트한 결과가 매직 넘버 0x5f3759df에서 감산되며, 더 정확한 근사를 위해 이 근사치를 다시 float 실수로 취급하여 뉴턴의 방법을 한 번 사용한다. 이 알고리즘은 float 실수의 나눗셈을 사용하는 것보다 거의 네 배 더 빠르게 역 제곱근을 계산한다.

존 카맥은 퀘이크 3 아레나에 사용된 소스 코드에 대해 이드 소프트웨어에서 퀘이크의 최적화를 도왔던 탁월한 어셈블리 프로그래머인 Terje Mathisen이 작성했을 것으로 보았다. 그러나 이러한 방법이 하드웨어 개발과 소프트웨어 개발 모두에서 뿌리 깊게 사용되고 있었던 것으로 나타났고, 알려진 최초의 사용으로는 게리 타롤리의 SGI Indigo를 위한 구현이 실리콘 그래픽스와 3dfx 인터랙티브를 거쳤다. Rys Sommefeldt는 Ardent Computer의 Greg Walsh가 MATLAB을 설계한 클리브 몰러와의 협의를 거쳐 이 알고리즘을 만들었다고 결론지었다.

개요

위 사진에 나타난 곡면 상의 표면에 수직인 벡터, 즉 법선 벡터는 조명과 셰이딩을 위한 계산에 광범위하게 사용된다.

역 제곱근은 단위 벡터를 구하는 데 사용되며, 법선 단위 벡터를 사용하여 표면을 향하는 빛의 입사각과 반사각을 결정할 수 있다. 3차원 벡터 ${\boldsymbol {v}}(v_{1},\ v_{2},\ v_{3})$의 크기를 유클리드 노름

$\|{\boldsymbol {v}}\|={\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}}$

으로 정하면 이 벡터와 같은 방향의 단위 벡터는 다음과 같다:

${\frac {\boldsymbol {v}}{\|{\boldsymbol {v}}\|}}\left({\frac {v_{1}}{\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}}},\ {\frac {v_{2}}{\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}}},\ {\frac {v_{3}}{\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}}}\right)$

여기서 실수 $x$에 대한 역 제곱근 ${\frac {1}{\sqrt {x}}}$이 등장한다. 3차원 컴퓨터 그래픽스에서는 이러한 연산을 매초 수백만 번 수행해야 했고, 이것은 특수한 하드웨어가 출현하기 전까지 번거로운 작업이었다. 1990년대 초에는 실수의 처리 능력이 정수 처리에 비해서 뒤처져 있었으나, 고속 역 제곱근 알고리즘은 float 실수의 나눗셈을 우회하여 이러한 문제 해소에 도움이 되었다. 퀘이크 3 아레나에서 그래픽 연산 속도를 개선하기 위해서 사용되었고, 이후 FPGA를 사용하는 일부 특수한 하드웨어 버텍스 셰이더에 구현되기도 했다.

코드

다음 코드는 퀘이크 3 아레나에 사용된 것으로, C 전처리 지시문을 생략했지만 본래의 정확한 주석까지 포함하였다.

float Q_rsqrt( float number )
{
	long i;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y = number;
	i = * ( long * ) &y;                     // evil floating point bit level hacking
	i = 0x5f3759df - ( i >> 1 );             // what the fuck?
	y = * ( float * ) &i;
	y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
//	y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

	return y;
}

1990년대 초에 역 제곱근을 계산하는 통상적인 방법은 첫 번째 근사치를 순람표로 정하는 것이었다. 이 코드는 순람표를 사용하는 것보다 빠르다는 것을 증명하였고, 다른 일반적인 알고리즘보다 네 배 정도 더 빨랐다. 이 코드에 사용되는 상수 0x5f3759df는 그 의미를 즉시 파악하기 어렵기 때문에 이 알고리즘의 매직 넘버라고 불린다. 이 알고리즘은 뉴턴의 방법을 사용하여 비교적 정확한 결과를 만들어 내지만, 소수점의 손실 때문에 부정확하고 1999년부터 도입된 x86 SSE의 rsqrtss에 비해 훨씬 느리다. 또한 long을 사용하는 것은 64비트 시스템에 대한 코드의 이식성을 떨어뜨릴 수 있다.

실행 예시

다음은 $x=0.15625$에 대해 $x^{-1/2}\approx 2.52982$를 근사하는 과정이다:

0_01111100_01000000000000000000000  = 1.25 * 2^-3, Bit pattern of both x and i
0_00111110_00100000000000000000000  = 1.125 * 2^-65, Shift right one position: (i >> 1)
0_10111110_01101110101100111011111  = 1.432430... * 2^+63, The magic number 0x5f3759df
0_10000000_01001110101100111011111  = 1.307430... * 2^+1, The result of 0x5f3759df - (i >> 1)

이 결과는 2.61486으로 약 3.4%의 오차를 갖는다. 뉴턴의 방법을 한 번 사용하면 2.52549로 약 0.17%의 오차를 갖게 된다.

알고리즘 설명

이 알고리즘은 다음 단계들을 수행해서 $x^{-1/2}$를 계산한다:

1. float 양수 $x$를 long 비트로 취급하면 $\log _{2}(x)$로 다룰 수 있다.
2. 매직 넘버를 사용하여 $-\log _{2}(x)/2$의 근사치를 계산한다.
3. long 비트 $\log _{2}(x^{-1/2})$를 다시 float 양수로 취급하면 $x^{-1/2}$로 다룰 수 있다.
4. 뉴턴의 방법을 한 번 사용하여 근사치를 수정한다.

부동소수점

 이 부분의 본문은 IEEE 754입니다.

이 알고리즘은 부동소수점의 단정밀도 표현에 의존하기 때문에, 이에 관한 설명이 필요하다. 한 예로, 십진수 실수 -118.625에 대한 부동소수점의 단정밀도 표현은 다음과 같다:

1. 음수이므로 sign은 1이 된다.
2. $118.625_{10}=1110110.101_{2}$이므로 $118.625_{10}=(2_{10})^{110_{2}}\times 1.110110101_{2}$와 같이 소수점 앞을 1로 만든다.
3. exponent는 음수 지수를 고려해서 지수에 127을 더한 $110_{2}+1111111_{2}=10000101_{2}$가 되고, fraction은 소수점 뒤의 23자리 11011010100000000000000가 된다.
4. 정리하면 다음과 같다:

그러므로 float 양수 $x$를 long 비트로 표현하면 $x=2^{\mathrm {E} }(1+\mathrm {M} )$에서 $\mathrm {E}$와 $\mathrm {M}$을 알 수 있고, 양변에 밑이 2인 로그를 취하면

$\log _{2}(x)=\mathrm {E} +\log _{2}(1+\mathrm {M} )$

이다.

로그 근사

여기서 범위 $0\leq \mathrm {M} <1$를 고려하여

$\log _{2}(1+\mathrm {M} )\approx \sigma +\mathrm {M}$

로 근사할 수 있는 상수 $\sigma$를 사용한다.

$(\log _{2}(1+\mathrm {M} )-\mathrm {M} )'={\frac {1}{\mathrm {M} \ln 2+\ln 2}}-1$

이고 극댓값은 $\mathrm {M} ={\frac {1}{\ln 2}}-1$에서 0.0860713이기 때문에, 오차의 최댓값이 가장 작은 상수는 그 절반인 $\sigma =0.0430357$이다.

float 양수에 대한 로그 근사(파란색)와 로그 함수(회색).

한편 양수 $x$를 비트 $n_{x}$로 취급하면

${\begin{aligned}n_{x}&=2^{23}(\mathrm {E} +127)+2^{23}\times \mathrm {M} \\&=2^{23}(\mathrm {E} +\sigma +\mathrm {M} )+2^{23}(127-\sigma )\\&\approx 2^{23}\log _{2}(x)+2^{23}(127-\sigma )\end{aligned}}$

이다. $\sigma =0.0430357$를 대입하면

${\begin{aligned}n_{x^{-1/2}}&\approx 2^{23}\log _{2}(x^{-1/2})+2^{23}(127-\sigma )\\&=3/2\times 2^{23}(127-\sigma )-n_{x}/2\\&\approx \mathrm {0x5f37bcb6} -n_{x}/2\end{aligned}}$

이다. 이것은 다음 코드로 쓰였다:

i = 0x5f3759df - ( i >> 1 ); // what the fuck?

매직 넘버 0x5f3759df가 어떻게 유도되었는지는 알려져 있지 않으며, 추론되는 상수는 $\sigma \approx 0.0450466$이다.

뉴턴의 방법

 이 부분의 본문은 뉴턴 방법입니다.

함수 $f(x)$에 대한 $(x_{0},\ f(x_{0}))$에서의 접선

$y=f'(x_{0})(x-x_{0})+f(x_{0})$

의 $x$절편은 $x_{0}$보다 방정식 $f(x)=0$의 해에 근접하며, 이 값을 다시 $x_{0}$에 대입하여 방정식의 해를 근사할 수 있다. 위의 단계는 이 방법에 대한 신뢰할 만한 초깃값을 제공하여 함수 $f(y)={\frac {1}{y^{2}}}-x$에 대한 $(y_{0},\ f(y_{0}))$에서의 접선

$z=-{\frac {2}{{y_{0}}^{3}}}(y-y_{0})+{\frac {1}{{y_{0}}^{2}}}-x$

의 $y$절편

$y_{1}=y_{0}\left({\frac {3}{2}}-{\frac {x}{2}}(y_{0})^{2}\right)$

으로 방정식 $f(x^{-1/2})=0$의 해인 역 제곱근을 근사할 수 있다. 이것은 다음 코드로 쓰였다:

y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration

정확도

libstdc의 역 제곱근과 휴리스틱적인 고속 역 제곱근의 차이.

오른쪽 그래프는 뉴턴의 방법을 한 번 사용하여 개선된 오차를 보여 준다. 0.01에 대해 표준 라이브러리는 10.0을, 위 알고리즘을 적용한 함수는 9.982522를 반환하며, 로그 스케일에서 서로의 차이는 일정한 영역을 벗어나지 않는다. Chris Lomont는 더 정확하게 근사하는 매직 넘버 0x5f375a86을, Matthew Robertson은 배정밀도 표현에서 유효한 매직 넘버 0x5fe6eb50c7b537a9를 찾아냈다.

문서

• Blinn, Jim (July 1997). “Floating Point Tricks”. 《Computer Graphics & Applications, IEEE》 17 (4): 80. doi:10.1109/38.595279. 
• Blinn, Jim (2003). 《Jim Blinn's Corner: Notation, notation notation》. Morgan Kaufmann. ISBN 1-55860-860-5. 
• Eberly, David (2001). 《3D Game Engine Design》. Morgan Kaufmann. ISBN 978-1-55860-593-0. 
• Eberly, David (2015). “Fast and Accurate Inverse Square Root” (PDF). Geometric Tools. 2009년 2월 24일에 원본 문서 (PDF)에서 보존된 문서. 2015년 5월 9일에 확인함. 
• Hennessey, John; Patterson, David A. (1998). 《Computer Organization and Design》 2판. San Francisco, CA: Morgan Kaufmann Publishers. ISBN 978-1-55860-491-9. 
• Kushner, David (August 2002). “The wizardry of Id”. 《IEEE Spectrum》 39 (8): 42–47. doi:10.1109/MSPEC.2002.1021943. 
• Lomont, Chris (February 2003). “Fast Inverse Square Root” (PDF). 2009년 2월 13일에 확인함. 
• McEniry, Charles (August 2007). “The Mathematics Behind the Fast Inverse Square Root Function Code” (PDF). 2015년 5월 11일에 원본 문서 (PDF)에서 보존된 문서. 2009년 2월 13일에 확인함. 
• Middendorf, Lars; Mühlbauer, Felix; Umlauf, George; Bodba, Christophe (2007년 6월 1일). 〈Embedded Vertex Shader in FPGA〉. Rettberg, Achin. 《Embedded System Design: Topics, Techniques and Trends》. IFIP TC10 Working Conference:International Embedded Systems Symposium (IESS). et al. Irvine, California: Springer. ISBN 978-0-387-72257-3. 
• Striegel, Jason (2008년 12월 4일). “Quake's fast inverse square root”. 《Hackszine》. 오라일리 미디어. 2009년 2월 15일에 원본 문서에서 보존된 문서. 2013년 1월 7일에 확인함. 

외부 링크

• A Brief History of InvSqrt by Matthew Robertson
• 0x5f3759df, further investigations into accuracy and generalizability of the algorithm by Christian Plesner Hansen
• Origin of Quake3's Fast InvSqrt()
• Quake III Arena source code
• Margolin, Tomer (2005년 8월 27일). “Magical Square Root Implementation In Quake III”. 《CodeMaestro》. The Coding Experience. 2012년 4월 14일에 원본 문서에서 보존된 문서. 2012년 5월 15일에 확인함.