빠른 역제곱근 알고리즘은 양의 실수 x에 대해 1/√x를 빠르게 근사하는 방법입니다. 이름에 ‘역’이 들어가지만 1/x를 구하는 알고리즘은 아닙니다. 32비트 부동소수점 수의 비트 표현을 정수처럼 읽고, 이 값을 절반으로 줄여 이른바 마법 상수 0x5f3759df에서 뺀 다음, 뉴턴–랩슨 보정식을 한 번 적용합니다.
이 코드는 id Software가 공개한 Quake III Arena 소스의 Q_rsqrt 함수로 널리 알려졌습니다. 당시 3차원 그래픽에서는 벡터를 정규화하는 계산이 매우 자주 필요했습니다. 벡터 길이 √(x²+y²+z²)를 직접 구한 뒤 나누는 대신, 길이 제곱의 역제곱근을 구해 각 성분에 곱하면 같은 정규화 결과를 얻습니다.
구하는 값은 y≈1/√x이며, 뉴턴 보정은 y←y(1.5−0.5xy²)입니다.
역제곱근은 어디에 쓰이나
벡터 v=(a,b,c)의 길이는 ‖v‖=√(a²+b²+c²)입니다. 단위 벡터를 만들려면 v/‖v‖를 계산합니다. s=a²+b²+c²라고 두면 각 성분에 1/√s를 곱할 수 있습니다. 조명 방향, 표면 법선, 물리 시뮬레이션처럼 방향만 필요할 때 이 계산이 반복됩니다.
| 계산 | 식 | 역할 |
|---|---|---|
| 길이 제곱 | s=a²+b²+c² | 제곱근 전의 값 |
| 역제곱근 | y=1/√s | 길이의 역수 |
| 정규화 | (ay,by,cy) | 길이 1인 방향 벡터 |
현대 프로세서와 컴파일러에서는 제곱근과 나눗셈 명령, SIMD 명령, 하드웨어 근삿값이 발전했습니다. 따라서 이 고전 코드를 그대로 쓰는 것이 언제나 더 빠르다고 볼 수 없습니다. 성능은 하드웨어, 컴파일러, 정확도 요구, 배치 크기에 따라 측정해야 합니다. 이 글의 핵심은 역사적으로 유명한 근사법의 수학 구조입니다.
IEEE 754 단정도 부동소수점의 구조
정규화된 양의 32비트 단정도 값은 부호 1비트, 지수 8비트, 분수 23비트로 저장됩니다. 값은 대략 (1+분수)×2^(지수−127) 형태입니다. 부호가 양수라면 지수 비트가 커질수록 값의 로그도 거의 선형으로 커집니다. 빠른 역제곱근은 이 ‘비트 패턴이 로그와 비슷하게 움직이는’ 성질을 이용합니다.
x=2ᵉm으로 쓰면 1≤m<2이고, log₂x=e+log₂m입니다. 역제곱근의 로그는 log₂(x⁻¹ᐟ²)=−½log₂x입니다. 즉 로그 영역에서는 입력의 로그를 절반으로 줄이고 부호를 바꾸는 일차 변환입니다. 부동소수점 비트를 정수로 해석한 값은 지수와 가수를 한 줄에 놓은 근삿값처럼 작동하므로 오른쪽으로 한 비트 이동한 i>>1이 대략 절반에 해당합니다.
마법 상수에서 비트 절반을 빼는 이유
i = 0x5f3759df − (i >> 1)
입력 x의 비트 패턴을 정수 i로 읽고 i를 한 비트 오른쪽으로 이동하면 정수값이 대략 절반이 됩니다. 이를 적절한 상수에서 빼면 로그 관점에서 −½log₂x에 맞는 지수와 가수의 근삿값을 얻습니다. 결과 비트를 다시 float로 읽으면 1/√x에 가까운 초깃값 y₀가 됩니다.
상수 0x5f3759df는 단순히 IEEE 754의 바이어스만 옮긴 값은 아닙니다. 가수 구간에서 log₂(1+f)를 선형으로 근사할 때 생기는 오차까지 고려해 경험적·수학적으로 조정된 상수입니다. 다른 오차 기준을 최소화하면 인접한 다른 상수가 제안될 수 있습니다. 그러므로 이 상수가 유일하게 가능한 값이라는 설명은 정확하지 않습니다.
뉴턴–랩슨 보정식 유도
비트 조작으로 만든 값은 빠르지만 아직 거친 근사입니다. 목표 y=1/√x는 1/y²=x를 만족합니다. 뉴턴법에 편한 함수 f(y)=1/y²−x를 잡을 수 있습니다. 도함수는 f′(y)=−2/y³입니다.
yₙ₊₁=yₙ−f(yₙ)/f′(yₙ)=yₙ(3−xyₙ²)/2=yₙ(1.5−0.5xyₙ²)
코드에서는 x2=0.5x를 미리 계산하고 y=y(1.5−x2·y²)로 적습니다. 초깃값이 충분히 가까우면 뉴턴법은 오차를 빠르게 줄입니다. Quake III 소스에는 이 보정을 한 번 수행하고 두 번째 반복은 주석 처리한 형태가 있습니다. 한 번은 속도와 정확도의 절충이었으며, 더 높은 정확도가 필요하면 추가 반복을 고려할 수 있습니다.
수치 예시: x=4
정확한 값은 1/√4=0.5입니다. 비트 단계에서 만든 초깃값이 예를 들어 0.48이라고 가정해 보겠습니다. 뉴턴 보정에 넣으면 y₁=0.48×(1.5−0.5×4×0.48²)입니다. 0.48²=0.2304이므로 괄호는 1.0392이고, y₁=0.498816입니다. 한 번의 보정으로 0.5에 훨씬 가까워집니다.
이 예의 0.48은 원리를 보여 주기 위한 가정값입니다. 실제 비트 단계의 값은 입력과 구현, 부동소수점 표현에 따라 정해집니다. 알고리즘을 평가할 때는 상대오차 |근삿값−정확값|/|정확값|을 여러 입력 구간에서 측정해야 합니다.
뉴턴 보정이 오차를 줄이는 방식
근삿값을 정확값에 대한 비율로 y=(1+δ)/√x라고 놓겠습니다. 보정식에 대입하면 새 비율은 (1+δ)(1.5−0.5(1+δ)²)입니다. 전개하면 1−1.5δ²−0.5δ³이 됩니다. 일차 오차 항 δ가 사라지고 주된 오차가 δ²에 비례합니다. 초깃값이 가까운 영역에서는 오차가 대략 제곱되어 빠르게 작아지는 이유입니다.
뉴턴법이 언제나 모든 초깃값에서 안전한 것은 아닙니다. 하지만 비트 근사가 양의 정상 범위 입력에서 목표값 가까이에 오도록 설계되어 있기 때문에 한 번의 반복만으로도 실용적인 근사가 됩니다. 입력 영역과 요구 정확도를 벗어나면 별도 처리가 필요합니다.
고전 C 코드의 동작 순서
- 양의 float 입력 x의 절반 x2=x×0.5를 저장합니다.
- x의 32비트 패턴을 정수 i로 해석합니다.
- i=0x5f3759df−(i>>1)을 계산합니다.
- 바뀐 비트를 다시 float y로 해석합니다.
- y=y(1.5−x2y²)로 뉴턴 보정을 한 번 합니다.
- y를 1/√x의 근삿값으로 반환합니다.
역사적 코드는 포인터 형 변환으로 같은 메모리를 long과 float로 읽었습니다. 현대 C와 C++에서 이런 별칭 위반 방식은 컴파일러 최적화 아래 정의되지 않은 동작을 일으킬 수 있습니다. 재현하려면 비트 복사나 표준이 허용하는 비트 캐스트를 사용하고, float가 실제로 IEEE 754 binary32인지 확인해야 합니다. 또한 C의 long 크기는 플랫폼마다 다를 수 있으므로 정확히 32비트인 정수형을 써야 합니다.
0, 음수, 무한대와 비정규수는 따로 생각해야 한다
실수 범위에서 1/√x는 x>0일 때 양의 실수입니다. x=0이면 값은 발산하며 IEEE 754 연산에서는 보통 양의 무한대와 관련된 처리가 필요합니다. x<0이면 실수 결과가 없어 NaN 처리가 필요합니다. 고전 코드는 이런 입력을 완전하게 처리하는 범용 수학 함수가 아닙니다.
비정규수, 무한대, NaN도 정상적인 양의 정규수와 비트 구조가 다르게 작동합니다. 정확한 라이브러리 함수를 대체하려면 예외 입력, 반올림, 플랫폼의 부동소수점 환경을 모두 다뤄야 합니다. 게임 엔진 내부에서 제한된 정상 입력을 빠르게 처리하던 함수와 표준 수학 라이브러리의 계약은 범위가 다릅니다.
빠른 역제곱근과 정확한 역제곱근의 차이
| 구분 | 빠른 고전 근사 | 표준 수학 계산 |
|---|---|---|
| 목표 | 낮은 비용의 충분한 근사 | 규정된 정확도와 예외 처리 |
| 핵심 | 비트 초깃값 + 뉴턴 보정 | 하드웨어·라이브러리 구현 |
| 입력 | 보통 양의 정상 float 가정 | 0, 음수, 무한대, NaN 고려 |
| 이식성 | 표현과 언어 규칙 확인 필요 | 표준 인터페이스 사용 |
시각 효과나 방향 벡터에서 작은 오차가 허용될 수 있지만, 과학 계산이나 누적 오차에 민감한 계산에서는 요구 정확도를 먼저 정해야 합니다. ‘빠르다’는 명칭도 역사적 맥락입니다. 현재 환경에서는 컴파일러가 제공하는 내장 함수와 하드웨어 명령을 벤치마크해 비교해야 합니다.
역제곱 법칙과 역제곱근은 다르다
이름이 비슷해 혼동하기 쉽습니다. 역제곱근은 x⁻¹ᐟ²=1/√x이고, 역제곱은 r⁻²=1/r²입니다. 중력의 크기가 거리 제곱에 반비례하는 뉴턴의 만유인력과 거리 역제곱 법칙에서 쓰는 1/r²과 이번 알고리즘의 1/√x는 다른 함수입니다.
다만 벡터 길이 제곱 x=r²를 입력한다면 1/√x=1/r이 됩니다. 벡터 정규화에서는 바로 이 관계를 사용합니다. 변수에 무엇이 들어 있는지를 확인하지 않고 지수만 보면 공식을 잘못 해석할 수 있습니다.
자주 하는 오해
- 빠른 역제곱근이 1/x를 계산한다고 설명합니다. 실제 목표는 1/√x입니다.
- 0x5f3759df가 어떤 입력에도 정확한 답을 만든다고 생각합니다. 이 단계는 초깃값 근사입니다.
- 뉴턴 보정 한 번이면 올바르게 반올림된 표준 함수 결과와 항상 같다고 가정합니다.
- 역사적 포인터 캐스팅 코드를 현대 C/C++에 그대로 복사해도 안전하다고 생각합니다.
- 현재의 모든 CPU에서도 표준 sqrt와 나눗셈보다 반드시 빠르다고 단정합니다.
빠른 역제곱근 FAQ
마법 상수는 왜 0x5f3759df인가요?
IEEE 754의 지수 바이어스와 가수의 로그 근사를 한 번에 보정해 1/√x에 가까운 비트 패턴을 만들도록 선택된 상수입니다. 오차 기준에 따라 근처의 다른 상수도 제안될 수 있으므로 유일한 수는 아닙니다.
뉴턴–랩슨법을 두 번 쓰면 더 정확한가요?
초깃값이 수렴 구간에 있으면 일반적으로 추가 반복이 오차를 더 줄입니다. 대신 곱셈과 덧셈 비용이 늘어납니다. 필요한 정확도와 실행 환경에 따라 반복 횟수를 정해야 합니다.
오늘날에도 이 코드를 사용해야 하나요?
자동으로 그렇지는 않습니다. 플랫폼의 역제곱근 명령, 컴파일러 내장 함수, 표준 라이브러리와 정확도·속도를 실제로 비교해야 합니다. 고전 코드는 근사 계산과 부동소수점 표현을 이해하는 좋은 사례이지만 범용 대체품은 아닙니다.
핵심 정리
빠른 역제곱근은 양의 float x에서 1/√x를 구하는 근사 알고리즘입니다. 입력 비트를 정수처럼 읽어 0x5f3759df−(i>>1)을 계산하면 로그 영역의 −½ 변환에 가까운 초깃값이 생깁니다. 이어 y←y(1.5−0.5xy²)라는 뉴턴 보정을 적용해 오차의 일차 항을 없앱니다.
역사적 구현의 속도는 당시 하드웨어 맥락에서 이해해야 합니다. 현대 코드에서는 비트 재해석의 언어 규칙과 예외 입력을 지키고, 표준 함수 및 하드웨어 명령과 벤치마크해야 합니다. 무엇보다 이 알고리즘은 1/x가 아닌 1/√x를 계산한다는 점이 출발점입니다.