Q_rsqrt: xấp xỉ nghịch đảo căn bậc hai trong game Quake III Arena

Quake III Arena là trò chơi bắn súng góc nhìn thứ nhất do id Software phát hành năm 1999. Năm 2005, id Software công bố mã nguồn game engine của trò chơi (nay có trên GitHub). Trong đó có một hàm viết bằng C để tính rất nhanh \(\frac{1}{\sqrt{x}}\):

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 f**k?
	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;
}

Bài viết này giải thích từng dòng của đoạn code trên: vì sao một phép dịch bit và một hằng số trông rất kỳ lạ lại cho ra gần đúng \(\frac{1}{\sqrt{x}}\).

Vì sao game cần tính nhanh 1/√x

Một thao tác game engine phải làm rất nhiều lần là chuẩn hóa (normalize) vector: đưa một vector về vector đơn vị cùng hướng, tức vector có độ dài bằng \(1\). Mỗi giây, trò chơi phải chuẩn hóa hàng triệu vector để tính ánh sáng, phản xạ và các hiệu ứng 3D.

Ví dụ với vector \(\vec v = (3, 2)\) trên mặt phẳng \(Oxy\): độ dài của nó là \(\sqrt{3^2 + 2^2} = \sqrt{13}\). Muốn có độ dài bằng \(1\), ta nhân vector với \(\frac{1}{\sqrt{13}}\):

\[\frac{1}{\sqrt{3^2 + 2^2}} \cdot (3, 2) \approx (0{,}832;\ 0{,}555)\]

Như vậy mỗi lần chuẩn hóa đều cần tính \(\frac{1}{\sqrt{x}}\), với \(x\) là bình phương độ dài của vector (ở ví dụ trên là \(x = 13\)).

Sao không tính 1 / sqrt(x) như bình thường?

Cách viết thông thường trong C là 1 / sqrt(x), với hàm sqrt trong thư viện math.h. Nhưng vào cuối những năm 1990, máy tính chậm hơn bây giờ rất nhiều, và phép chia / cùng hàm sqrt chậm hơn hẳn các phép +, -, *. Vì vậy các lập trình viên game tìm một cách tính nhanh hơn, chỉ dùng những phép tính rẻ.

Ba phần của hàm Q_rsqrt

Các chú thích trong code không giải thích cách tính, nhưng cho ta thấy hàm gồm 3 phần:

1. evil floating point bit level hacking: đọc các bit của số thực như một số nguyên.
2. what the f**k?: phép trừ từ hằng số 0x5f3759df.
3. 1st iteration: một bước tính lại để kết quả chính xác hơn.

Để hiểu từng phần, trước hết ta cần biết máy tính lưu một số thực như thế nào.

Phần I: IEEE-754, cách máy tính lưu số thực

Dữ liệu trong máy tính được lưu thành các bit, mỗi bit là \(0\) hoặc \(1\). Theo chuẩn IEEE-754, một số thực kiểu float dùng \(32\) bit, chia làm \(3\) phần: bit dấu, phần mũ và phần định trị.

Bit dấu

Bit đầu tiên lưu dấu: S bằng \(0\) là số dương, bằng \(1\) là số âm. Khi chuẩn hóa vector, \(x\) là bình phương độ dài nên luôn dương. Vì vậy từ đây ta luôn coi S bằng \(0\) cho đơn giản.

Phần mũ (exponent)

IEEE-754 lưu số thực theo dạng khoa học (scientific notation) trong hệ nhị phân, tức viết số dưới dạng \(1{,}\ldots \cdot 2^{k}\). Ví dụ:

• Số \(0{,}001001_2\) viết lại thành \(1{,}001_2 \cdot 2^{-3}\).
• Số \(101{,}001_2\) viết lại thành \(1{,}01001_2 \cdot 2^{2}\).

\(8\) bit tiếp theo lưu số mũ \(k\) của \(2\). Nhưng \(8\) bit chỉ biểu diễn được các số từ \(0\) đến \(255\), trong khi số mũ có thể âm như ở ví dụ đầu. Vì vậy IEEE-754 lưu \(E = k + 127\); muốn lấy lại số mũ thật thì lấy \(E\) trừ \(127\):

• \(0{,}001001_2 = 1{,}001_2 \cdot 2^{124 - 127}\), nên phần mũ là \(E = 124\), viết bằng \(8\) bit là 01111100.
• \(101{,}001_2 = 1{,}01001_2 \cdot 2^{129 - 127}\), nên phần mũ là \(E = 129\), viết bằng \(8\) bit là 10000001.

Phần định trị (mantissa)

\(23\) bit còn lại lưu phần đứng sau dấu phẩy của số \(1{,}\ldots\) trong dạng khoa học. Chữ số \(1\) đứng đầu luôn có nên không cần lưu:

• \(1{,}001_2 \cdot 2^{-3}\) có \(23\) bit cuối là 00100000000000000000000.
• \(1{,}01001_2 \cdot 2^{2}\) có \(23\) bit cuối là 01001000000000000000000.

Hai công thức cần nhớ

Gọi \(E\) là phần mũ và \(M\) là phần định trị (đọc như một số nguyên \(23\) bit) của số thực dương \(x\). Khi đó:

• Giá trị của số thực: \(x = \left(1 + \frac{M}{2^{23}}\right) \cdot 2^{E - 127}\).
• Nếu đọc nguyên \(32\) bit của \(x\) như một số nguyên, ta được \(I_x = E \cdot 2^{23} + M\).

Ví dụ \(x = 0{,}001001_2\) (tức \(0{,}140625\) trong hệ thập phân) có biểu diễn 0 01111100 00100000000000000000000, nên \(I_x = 124 \cdot 2^{23} + 2^{20} = 1\,041\,235\,968\).

Phần II: Đọc bit số thực ra logarit cơ số 2

Ý tưởng của Q_rsqrt không tự nhiên lắm. Để hiểu nó, ta thử tính \(\log_2(x)\):

\[\begin{aligned} \log_2(x) &= \log_2 \left[ \left(1 + \frac{M}{2^{23}} \right) \cdot 2^{E - 127} \right] \\ &= \log_2 \left( 1 + \frac{M}{2^{23}} \right) + E - 127 \end{aligned}\]

Với \(t\) trong đoạn \([0, 1]\), ta có phép xấp xỉ \(\log_2(1 + t) \approx t + \epsilon\), trong đó \(\epsilon\) là một hằng số nhỏ được chọn để sai số trên cả đoạn nhỏ nhất. Hàm Q_rsqrt dùng \(\epsilon \approx 0{,}0450466\). Thay vào, ta được:

\[\begin{aligned} \log_2(x) &\approx \frac{M}{2^{23}} + \epsilon + E - 127 \\ &= \frac{1}{2^{23}} \left( E \cdot 2^{23} + M \right) + \epsilon - 127 \\ &= \frac{1}{2^{23}} \cdot I_x + \epsilon - 127 \end{aligned}\]

Nói cách khác: đọc 32 bit của số thực như một số nguyên là đã gần có \(\log_2(x)\), chỉ cần chia cho \(2^{23}\) và cộng thêm một hằng số.

Từ đây, đặt \(y = \frac{1}{\sqrt{x}}\). Nếu tính được \(I_y\) (số nguyên có được khi đọc 32 bit của \(y\)), ta sẽ có \(y\). Mục tiêu bây giờ là tính \(I_y\).

Phần III: "evil floating point bit level hacking"

Quay lại hàm Q_rsqrt. Các dòng đầu chỉ khai báo biến và gán giá trị. Ta xem hai dòng sau:

y  = number;
i  = * ( long * ) &y;                       // evil floating point bit level hacking

Dòng đầu gán y = number, tức biến y lúc này đang giữ số \(x\) cần tính. Dòng thứ hai cần lấy \(I_x\): đọc nguyên 32 bit của y như một số nguyên. Cách đầu tiên ta nghĩ tới có thể là i = (long)y;.

Nhưng phép ép kiểu này đổi giá trị của số thực sang số nguyên (bỏ phần thập phân), chứ không giữ các bit. Ví dụ số \(0{,}5\) có biểu diễn 0 01111110 00000000000000000000000. Đọc dãy bit này như số nguyên ta được \(1\,056\,964\,608\), còn phép ép kiểu (long)0.5 chỉ cho \(0\). Đây không phải kết quả ta cần.

Để đọc đúng các bit, code gốc "lách luật" bằng con trỏ, qua ba bước:

1. &y: lấy địa chỉ của biến y trong bộ nhớ. Đây là địa chỉ của một số thực.
2. (long *) &y: "đánh lừa" C rằng đây là địa chỉ của một số nguyên kiểu long.
3. * (long *) &y: lấy giá trị tại địa chỉ đó (dereference). C đọc vùng nhớ chứa các bit IEEE-754 của y như một số nguyên.

Cách này cho đúng \(I_x\) khi long có \(32\) bit, bằng float, như trên các trình biên dịch Quake III dùng năm 1999. Trên Linux 64 bit, long có \(64\) bit; phần VII nói cách viết đúng ngày nay.

Phần IV: Hằng số ma thuật 0x5f3759df

Ta cần tính \(y = \frac{1}{\sqrt{x}}\). Lấy logarit cơ số \(2\) hai vế:

\[\log_2(y) = \log_2\left(x^{-1/2}\right) = -\tfrac{1}{2} \log_2(x)\]

Thay công thức ở phần II vào cả hai vế:

\[\begin{gathered} \tfrac{1}{2^{23}} I_y + \epsilon - 127 \approx -\tfrac{1}{2} \left( \tfrac{1}{2^{23}} I_x + \epsilon - 127 \right) \\ \Rightarrow\ I_y \approx \tfrac{3}{2} \cdot 2^{23} \, (127 - \epsilon) - \tfrac{1}{2} I_x \end{gathered}\]

Với \(\epsilon \approx 0{,}0450466\), số \(\frac{3}{2} \cdot 2^{23} \, (127 - \epsilon) \approx 1\,597\,463\,007\), viết ở hệ \(16\) chính là 0x5f3759df. Còn \(\frac{1}{2} I_x\) tính được bằng phép dịch phải i >> 1. Phép dịch bỏ phần lẻ (7 >> 1 = 3, 6 >> 1 = 3), nhưng sai số này rất nhỏ so với \(I_x\) nên chấp nhận được.

Vậy hai dòng sau tính \(I_y\) rồi đọc ngược 32 bit đó thành số thực, cũng bằng mẹo con trỏ như phần III:

i  = 0x5f3759df - ( i >> 1 );               // what the f**k?
y  = * ( float * ) &i;

Đến đây ta đã có \(y \approx \frac{1}{\sqrt{x}}\), với sai số tương đối lớn nhất khoảng 3,4%.

Phần V: Phương pháp Newton tìm nghiệm

Để hiểu phần còn lại, ta cần phương pháp Newton (Newton's method), một cách tìm nghiệm gần đúng của phương trình \(f(y) = 0\).

Ý tưởng: bắt đầu từ một giá trị \(y_0\) mà ta nghĩ là gần nghiệm. Kẻ tiếp tuyến của đồ thị \(f\) tại \(y_0\); điểm tiếp tuyến cắt trục hoành là \(y_1\), thường gần nghiệm hơn \(y_0\):

\[y_1 = y_0 - \frac{f(y_0)}{f'(y_0)}\]

Lặp lại công thức này, ta được \(y_2, y_3, \ldots\) ngày càng gần nghiệm. Nếu \(y_0\) đã khá gần nghiệm thì chỉ một bước là đủ tốt. Có thể đọc thêm về phương pháp Newton ở Paul's Online Notes.

Phần VI: Một bước Newton làm kết quả chính xác hơn

Phần còn lại của hàm là hai dòng sau (dòng thứ hai đã bị tắt bằng //):

y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

Tìm \(y\) chính là giải phương trình:

\[y = \frac{1}{\sqrt{x}} \iff y^{-2} - x = 0\]

Đặt \(f(y) = y^{-2} - x\), ta có \(f'(y) = -2y^{-3}\). Áp dụng một bước Newton từ giá trị \(y_0\) có được ở phần IV:

\[\begin{aligned} y_1 &= y_0 - \frac{y_0^{-2} - x}{-2y_0^{-3}} = y_0 + \frac{y_0 - x y_0^{3}}{2} \\ &= y_0 \left( \frac{3}{2} - \frac{x}{2} \cdot y_0^{2} \right) \end{aligned}\]

Với threehalfs \(= 1{,}5\) và x2 \(= \frac{x}{2}\), công thức này chỉ cần phép nhân và phép trừ: y = y * ( threehalfs - ( x2 * y * y ) );. Vì \(y_0\) đã khá gần đúng nên một bước là đủ, và bước thứ hai trong code gốc bị tắt. Sau một bước, sai số tương đối lớn nhất giảm từ khoảng 3,4% xuống còn khoảng 0,175%.

Phần VII: Undefined behaviour và cách viết đúng ngày nay

Mẹo đọc bit bằng con trỏ ở phần III, theo chuẩn C và C++, là undefined behaviour: đọc vùng nhớ của một float qua con trỏ long * vi phạm quy tắc strict aliasing, và trình biên dịch được phép tối ưu theo cách làm hỏng kết quả. Thêm nữa, long không phải lúc nào cũng có \(32\) bit: trên Linux 64 bit nó có \(64\) bit, nên đọc 4 byte của float vào một long 8 byte sẽ ra số sai.

Cách viết an toàn là chép đúng 4 byte bằng std::memcpy vào kiểu std::uint32_t (luôn có đúng 32 bit). Từ C++20 có thể dùng std::bit_cast<std::uint32_t>(number). Trong Python, module struct làm việc tương tự:

#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstring>

static_assert(sizeof(float) == 4, "float phải là số thực 32 bit IEEE-754");

float q_rsqrt(float number) {
    const float x2 = number * 0.5f;
    std::uint32_t i;
    std::memcpy(&i, &number, sizeof i);  // đọc 32 bit của số thực như một số nguyên
    i = 0x5f3759df - (i >> 1);           // hằng số "ma thuật"
    float y;
    std::memcpy(&y, &i, sizeof y);       // đọc ngược 32 bit đó thành số thực
    y = y * (1.5f - x2 * y * y);         // một bước Newton
    return y;
}

int main() {
    const float xs[] = {0.25f, 2.0f, 13.0f, 100.0f, 12345.0f};
    for (float x : xs) {
        float nhanh = q_rsqrt(x), chuan = 1.0f / std::sqrt(x);
        std::printf("x = %8.2f   q_rsqrt = %.6f   1/sqrt = %.6f   sai số = %.4f%%\n",
                    x, nhanh, chuan, std::fabs(nhanh - chuan) / chuan * 100);
    }
}

import math
import struct


def q_rsqrt(number: float) -> float:
    i = struct.unpack("<I", struct.pack("<f", number))[0]  # đọc 32 bit của số thực như một số nguyên
    i = 0x5F3759DF - (i >> 1)                              # hằng số "ma thuật"
    y = struct.unpack("<f", struct.pack("<I", i))[0]       # đọc ngược 32 bit đó thành số thực
    return y * (1.5 - 0.5 * number * y * y)                # một bước Newton


for x in (0.25, 2.0, 13.0, 100.0, 12345.0):
    nhanh, chuan = q_rsqrt(x), 1 / math.sqrt(x)
    print(f"x = {x:8.2f}   q_rsqrt = {nhanh:.6f}   1/sqrt = {chuan:.6f}   sai số = {abs(nhanh - chuan) / chuan * 100:.4f}%")

Chạy chương trình, ta được:

Kết quả luôn nhỏ hơn giá trị đúng một chút, và sai số không vượt quá khoảng 0,175% với mọi số thực dương.

Nhanh đến đâu, ngày nay còn dùng không?

Theo phân tích của Chris Lomont (2003), trên máy tính thời đó Q_rsqrt nhanh hơn cách tính 1 / sqrt(x) khoảng \(4\) lần, với sai số tương đối lớn nhất khoảng 0,175%. Trong game, sai số này không nhìn thấy được trên màn hình, nên Quake III dùng hàm này cho rất nhiều phép chuẩn hóa vector.

Ngày nay CPU có sẵn lệnh tính xấp xỉ \(\frac{1}{\sqrt{x}}\) trong phần cứng, ví dụ lệnh rsqrtss trong bộ lệnh SSE của Intel, vừa nhanh hơn vừa chính xác hơn. Vì vậy Q_rsqrt không còn được dùng để tăng tốc. Giá trị của nó nằm ở bài học: hiểu cách máy tính lưu số thực, và thấy toán học (logarit, phương pháp Newton) biến thành vài dòng code chạy rất nhanh.

Luyện ngay trên Big-O Coding

Thử áp dụng các ý trong bài: tính căn bằng tìm kiếm nhị phân và làm việc với từng bit của số nguyên.

Nguồn tham khảo

• id Software: mã nguồn Quake III Arena (GitHub)
• Beyond3D: Origin of Quake3's Fast InvSqrt()
• Chris Lomont: Fast Inverse Square Root (2003)
• Paul's Online Notes: Newton's Method
• cppreference: Undefined behavior
• cppreference: std::bit_cast
• Nemean: Fast Inverse Square Root — A Quake III Algorithm (YouTube)
• Nick: Clever Code: Fast Inverse Square Root (YouTube)