Về một comment trong source code của Quake III Arena

Trong lập trình đồ họa 3D, chúng ta thường xuyên phải tính toán các vector đơn vị (unit vector) để xử lý ánh sáng và phản chiếu. Việc này đòi hỏi ta phải chuẩn hóa vector, và do đó cần thực hiện các phép tính 1∑i=1NΔxi2\frac{1}{\sqrt{\sum_{i=1}^N \Delta x_i^2}}. Vào những năm 1990, dòng mã huyền thoại này đã xuất hiện trong mã nguồn của trò chơi Quake III, đã loại bỏ những chỉ thị tiền xử lý của ngôn ngữ C, nhưng comment thì được giữ nguyên:

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;
}

Trong một thời gian dài cộng đồng nghi ngờ John Carmack là tác giả, nhưng sau khi có thêm thông tin thì điều này được xem là không chính xác. Theo Wikipedia thì Greg Walsh là người viết thuật toán ban đầu, sau này được đưa vào Quake III (có thể thông qua Brian Hook từ 3dfx). Thuật toán được biết đến đưới cái tên “Fast inverse square root” (Căn bậc hai nghịch đảo nhanh). Nó nhanh một cách đáng kinh ngạc so với các phương pháp thông thường thời bấy giờ.

Thủ thuật này ước lượng căn nghịch đảo 1x\frac{1}{\sqrt{x}} bằng phương pháp xấp xỉ Newton, và bắt đầu với một dự đoán ban đầu cực kỳ chính xác. Các bước như sau:

Bước 1: Ép kiểu đối số xx sang số nguyên (Aliasing)

Máy tính lưu trữ số thực dấu phẩy động (float) theo chuẩn IEEE 754, bao gồm ba phần: Dấu (Sign), Phần mũ (Exponent - E) và Phần lẻ (Mantissa - M).

Về cơ bản, một số thực được biểu diễn gần giống như:

Giaˊ trị≈2E⋅(1+M)\text{Giá trị} \approx 2^E \cdot (1 + M)
Nếu bạn coi các bit của số thực này như một số nguyên (integer), định dạng của nó trông rất giống với hàm logarit. Khi chúng ta thực hiện dòng mã i = * ( long * ) &y;, chúng ta không chuyển đổi giá trị (như từ 1.0 sang 1), mà chúng ta đang diễn giải lại các bit của số thực đó như một số nguyên.

Một cách kỳ diệu, giá trị số nguyên này xấp xỉ với log⁡2(x)\log_2(x).

Bước 2: Tính toán logarit

Sử dụng giá trị xấp xỉ này để tính toán xấp xỉ của:log⁡2(1x)=−12log⁡2(x)\log_2\left(\frac{1}{\sqrt{x}}\right) = -\frac{1}{2}\log_2(x) Dòng mã “hack” chính là đây:

i = 0x5f3759df - ( i >> 1 );

Dòng này thực chất đang thực hiện phép tính logarit mà chúng ta đã nói ở trên:

  1. (i >> 1) tương đương với việc chia logarit cho 2 (tức là nhân với 1/21/2).
  2. Phép trừ lấy từ số 0x5f3759df giúp xử lý phần hằng số trong công thức logarit và sự dịch chuyển của phần mũ trong định dạng float.

Kết quả của dòng này cung cấp cho chúng ta một lời giải “xấp xỉ cực tốt” cho 1x\frac{1}{\sqrt{x}}. Chúng ta đã đi từ một số thực, chuyển sang số nguyên để tính toán logarit nhanh chóng, và giờ cần chuyển ngược lại về số thực.

Bước 3: Ép kiểu ngược lại về số thực (Alias back)

Chuyển đổi kết quả số nguyên thu được trở lại dạng số thực. Đây là cách để thực hiện tính toán xấp xỉ cho hàm mũ cơ số 2, đưa giá trị từ miền logarit về lại miền giá trị thực ban đầu.

Bước 4: Tinh chỉnh xấp xỉ

Sau khi chuyển đổi ngược lại y = * ( float * ) &i;, chúng ta có một kết quả khá gần đúng. Để làm nó chính xác hơn, thuật toán sử dụng một vòng lặp của Phương pháp Newton.

Nói qua một chút về phương pháp này.

Giá trị y=1xy = \frac{1}{\sqrt{x}} chính là nghiệm của phương trình 1y2−x=0\frac{1}{y^2} - x = 0. Kết quả xấp xỉ thu được từ các bước trước đó có thể được tinh chỉnh bằng cách sử dụng phương pháp tìm nghiệm của hàm số. Nếu ta có một giá trị xấp xỉ yny_n của yy, thì một giá trị xấp xỉ tốt hơn yn+1y_{n+1} có thể được tính toán bằng công thức yn−f(yn)f′(yn)y_n - \frac{f(y_n)}{f'(y_n)}, trong đó f′(yn)f'(y_n) là đạo hàm của f(y)f(y) tại điểm yny_n.

Áp dụng vào phương trình 1y2−x=0\frac{1}{y^2} - x = 0, phương pháp Newton cho ta:

f(y)=1y2−xf(y) = \frac{1}{y^2} - x

f′(y)=−2y3f'(y) = -\frac{2}{y^3}

yn+1=yn−f(yn)f′(yn)=yn+yn32(1yn2−x)=yn(32−x2yn2)y_{n+1} = y_n - \frac{f(y_n)}{f'(y_n)} = y_n + \frac{y_n^3}{2} \left( \frac{1}{y_n^2} - x \right) = y_n \left( \frac{3}{2} - \frac{x}{2} y_n^2 \right)

Công thức này được viết trong mã nguồn dưới dạng:
y = y * ( threehalfs - ( x2 * y * y ) );

Phép tính này giúp “hiệu chỉnh” sai số, đưa giá trị xấp xỉ ban đầu về rất gần với kết quả chính xác tuyệt đối. Trong Quake III, một lần lặp là đủ để đạt độ chính xác cần thiết cho đồ họa.

Tại sao nó lại quan trọng?

Vào thời điểm đó, bộ xử lý (CPU) tính toán phép chia và căn bậc hai rất chậm. Thuật toán này đã thay thế chúng bằng:

  • Một phép dịch bit (>>).
  • Một phép trừ số nguyên.
  • Một vài phép nhân số thực.

Đây là một sự đánh đổi thông minh giữa độ chính xác cực cao và tốc độ xử lý, cho phép các trò chơi 3D đời đầu chạy mượt mà trên phần cứng hạn chế.

---

Ghi chú và Giải thích thêm

  1. Thủ thuật ép kiểu (Pointer Hacking): Dòng i = * ( long * ) &y; là cách để “đánh lừa” trình biên dịch, lấy nội dung nhị phân của một biến float và gán vào một biến long. Trong C hiện đại, điều này có thể vi phạm quy tắc “strict aliasing”, người ta thường dùng memcpy để đạt được mục đích tương tự mà vẫn đảm bảo tính an toàn.
  2. Độ chính xác: Con số ma thuật 0x5f3759df được cho là của Greg Walsh, mặc dù nó thường gắn liền với John Carmack (người sáng lập id Software). Các nhà toán học sau này đã tìm ra những con số thậm chí còn tối ưu hơn một chút, nhưng sự chênh lệch là không đáng kể.
  3. Tính thời đại: Ngày nay, các CPU hiện đại đã có các tập lệnh chuyên dụng (như RSQRTSS trên x86) để tính căn bậc hai nghịch đảo trực tiếp ở cấp độ phần cứng, nên thuật toán này không còn được sử dụng rộng rãi trong mã nguồn mới, nhưng nó vẫn mãi là một biểu tượng của sự sáng tạo tối ưu hóa trong lập trình.

Xem thêm https://en.wikipedia.org/wiki/Fast_inverse_square_root