고속 역 제곱근

고속 역 제곱근(高速逆-根, 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차원 벡터 formula_1의 크기를 유클리드 노름
으로 정하면 이 벡터와 같은 방향의 단위 벡터는 다음과 같다:
여기서 실수 formula_4에 대한 역 제곱근 formula_5이 등장한다. 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 * ) &amp;y; // evil floating point bit level hacking
 i = 0x5f3759df - ( i » 1 ); // what the fuck?
 y = * ( float * ) &amp;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의 codice_1에 비해 훨씬 느리다. 또한 codice_2을 사용하는 것은 64비트 시스템에 대한 코드의 이식성을 떨어뜨릴 수 있다.
실행 예시.
다음은 formula_6에 대해 formula_7를 근사하는 과정이다:
 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%의 오차를 갖게 된다.
알고리즘 설명.
이 알고리즘은 다음 단계들을 수행해서 formula_8를 계산한다:
부동소수점.
이 알고리즘은 부동소수점의 단정밀도 표현에 의존하기 때문에, 이에 관한 설명이 필요하다. 한 예로, 십진수 실수 -118.625에 대한 부동소수점의 단정밀도 표현은 다음과 같다:
그러므로 float 양수 formula_4를 long 비트로 표현하면 formula_18에서 formula_19와 formula_20을 알 수 있고, 양변에 밑이 2인 로그를 취하면
이다.
로그 근사.
여기서 범위 formula_22를 고려하여
로 근사할 수 있는 상수 formula_24를 사용한다.
이고 극댓값은 formula_26에서 0.0860713이기 때문에, 오차의 최댓값이 가장 작은 상수는 그 절반인 formula_27이다.
한편 양수 formula_4를 비트 formula_29로 취급하면
이다. formula_27를 대입하면
이다. 이것은 다음 코드로 쓰였다:
i = 0x5f3759df - ( i » 1 ); // what the fuck?
매직 넘버 0x5f3759df가 어떻게 유도되었는지는 알려져 있지 않으며, 추론되는 상수는 formula_33이다.
뉴턴의 방법.
함수 formula_34에 대한 formula_35에서의 접선
의 formula_4절편은 formula_38보다 방정식 formula_39의 해에 근접하며, 이 값을 다시 formula_38에 대입하여 방정식의 해를 근사할 수 있다. 위의 단계는 이 방법에 대한 신뢰할 만한 초깃값을 제공하여 함수 formula_41에 대한 formula_42에서의 접선