- Fast inverse square root nổi tiếng nhờ Quake 3 là một giải pháp hiệu năng của thời đó, dùng diễn giải lại bit
float và hiệu chỉnh Newton-Raphson để xấp xỉ nhanh 1 / sqrt(x)
- Điểm cốt lõi là mẫu bit số nguyên của
float IEEE-754 32 bit có thể được xử lý như một giá trị xấp xỉ log2(x) đã được scale và shift
0x5f3759df - (i >> 1) là cách chuyển log2(x^-0.5) = -0.5 * log2(x) sang phép dịch bit số nguyên và phép trừ, còn hằng số ma thuật xuất phát từ 3/2 * 2^23 * (127 - σ)- Sau đó áp dụng hiệu chỉnh Newton-Raphson bằng một lần
y = y * (1.5 - 0.5x * y * y), và vòng lặp thứ hai trong mã Quake đã bị comment out
- Năm 1999, nghịch đảo căn bậc hai cần đến hàng trăm đến hàng nghìn lần mỗi giây cho chiếu sáng và chuẩn hóa vector 3D, nhưng trên phần cứng hiện đại tính thực dụng của cùng mẹo này đã giảm đi nhờ xử lý dấu phẩy động chuyên dụng
Mã Quake thực hiện điều gì
float Q_rsqrt(float number) {
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = *(long*)&y;
i = 0x5f3759df - ( i >> 1 );
y = *(float*)&i;
y = y * ( threehalfs - ( x2 * y * y ) );
return y;
}
- Hàm này tính giá trị xấp xỉ của nghịch đảo căn bậc hai
1 / sqrt(number) cho number
- Phần nổi tiếng nhất là thao tác bit: diễn giải giá trị
float như long, rồi thực hiện 0x5f3759df - (i >> 1)
- Khi Quake 3 ra mắt năm 1999, nghịch đảo căn bậc hai là phép toán chậm và tốn kém, nhưng lại cần hàng trăm đến hàng nghìn lần mỗi giây trong các phép tính vector 3D cho phương trình chiếu sáng và chuẩn hóa
- Trên phần cứng hiện đại, các phép tính kiểu này hoặc không được thực hiện trên CPU, hoặc nếu có thì cũng nhanh nhờ phần cứng dấu phẩy động chuyên dụng đã phát triển hơn nhiều
Biểu diễn float IEEE-754 32 bit
float 32 bit gồm ba phần
- Sign: 1 bit, biểu thị số dương hay âm
- Exponent: 8 bit, xác định khoảng mà giá trị thuộc về
- Mantissa: 23 bit, biểu thị tuyến tính vị trí bên trong khoảng đó
- Một giá trị thông thường được diễn giải theo dạng sau
N = (-1)^S * 2^(E - 127) * (1 + M / 2^23)
B = 127 là giá trị bias dùng cho biased exponent, và số mũ thực là e = E - B- Mantissa không được dùng như phép nhân trực tiếp với
m, mà ở dạng 1 + m
- Nếu
m = 0 thì là 2^e
- Khi
m tiến gần 1, nó biểu diễn đến ngay trước khoảng số mũ tiếp theo là 2^(e+1)
- Nếu toàn bộ bit số mũ đều bằng 0 thì đó là số sub-normal, và công thức sẽ khác
N = (-1)^S * 2^-126 * m
- Sub-normal cần thiết để biểu diễn số 0 và các số rất nhỏ gần 0
- Nếu toàn bộ bit số mũ đều bằng 1 thì được xử lý là giá trị đặc biệt
E = 255, M = 0 thì là Infinity hoặc -InfinityM != 0 thì là NaN
Quan hệ log xuất hiện khi xem bit float như số nguyên
- Nếu xem biểu diễn nội bộ của
float như số nguyên 32 bit, ta có thể viết nó như sau
I_x = 2^31 S + 2^23 E + M
- Vì nghịch đảo căn bậc hai áp dụng cho đầu vào dương, có thể đặt
S = 0 để đơn giản hóa công thức
L = 2^23
I_x = L E + M
- Trong cùng một khoảng exponent, mantissa biểu thị vị trí tuyến tính, nhưng khi exponent tăng thì cùng số lượng bước mantissa sẽ phủ một đoạn rộng hơn trên trục số
E = 127, tức e = 0, tương ứng xấp xỉ khoảng [1, 2)E = 128, tức e = 1, tương ứng xấp xỉ khoảng [2, 4)- Cả hai khoảng đều có cùng số bước mantissa, nhưng khoảng thứ hai rộng gấp đôi
- Vì cấu trúc này, nếu nhìn mẫu bit thô của
float như số nguyên thì sẽ xuất hiện quan hệ mang tính logarit
Bit thô là giá trị xấp xỉ của log2(x)
- Nếu diễn giải mẫu bit của
float thành số nguyên I_x, có thể xem nó như một phép xấp xỉ tuyến tính từng đoạn của log2(x)
- Quan hệ này được biểu diễn bằng công thức xấp xỉ sau
log2(x) ≈ I_x / L - B
- Chia số nguyên bit thô cho kích thước mantissa
L = 2^23, rồi trừ exponent bias B = 127, ta sẽ thu được giá trị gần với log2(x)
- Phần log trong từng đoạn mantissa được xử lý bằng xấp xỉ tuyến tính
log2(1 + x) ≈ x + σ
σ là tham số tuning để điều chỉnh phép xấp xỉ, còn x biểu thị vị trí bên trong khoảng exponent trong miền [0, 1]
Chuyển nghịch đảo căn bậc hai thành đồng nhất thức log
- Mục tiêu là tính giá trị sau
y = 1 / sqrt(x)
- Có thể viết lại dưới dạng lũy thừa như sau
y = x^-0.5
- Áp dụng đồng nhất thức log, phép tính nghịch đảo căn bậc hai trở thành quan hệ sau
log2(1 / sqrt(x)) = log2(x^-0.5) = -0.5 * log2(x)
- Dựa vào việc bit
float hoạt động như một giá trị xấp xỉ của log2(x), ta có thể xấp xỉ trực tiếp biểu diễn bit số nguyên I_y của y từ biểu diễn bit số nguyên I_x của x
I_y ≈ -0.5 I_x + 1.5 L (B - σ)
- Công thức này dẫn thẳng đến dòng mã cốt lõi trong Quake
i = 0x5f3759df - ( i >> 1 );
i >> 1 dịch phải số nguyên 1 bit, đóng vai trò nhân với 1/2- Hằng số phía trước
0x5f3759df tương ứng với 1.5 * L * (B - σ)
Bản chất của hằng số 0x5f3759df
- Nếu đặt
σ = 0, hằng số được tính như sau
1.5 * 2^23 * 127 = 1598029824
- Biểu diễn hệ 16 của giá trị này là
0x5f400000
- Nó khác với hằng số thực tế của Quake là
0x5f3759df đúng 566817
- Từ chênh lệch này, có thể tính ra giá trị
σ tương ứng với mã Quake như sau
σ = 377878 / 2^23
σ = 0.04504656
- Trong C, có thể tính cùng hằng số đó như sau
int32_t compute_magic(void) {
double sigma = 0.0450465;
double expression = 1.5 * pow(2.0, 23.0) * (127.0 - sigma);
int32_t i = expression;
return i;
}
// -> 0x5f3759df
- Ở đây dùng
double, và phép chuyển sang số nguyên là casting thông thường chứ không phải diễn giải lại bit
- Giá trị
σ này được chọn để tối ưu hóa phép xấp xỉ, nhưng không phải giá trị tối ưu thực sự, và cũng không rõ chính xác ai là người tạo ra nó
Vì sao đây không chỉ là một mẹo hack đơn giản
0x5f3759df - (i >> 1) là công thức tạo giá trị khởi tạo cho nghịch đảo căn bậc hai bằng cách tận dụng thực tế rằng bit thô của float là một phép xấp xỉ log- Dù dựa trên quan hệ toán học phức tạp, ở bước thực thi nó chỉ dùng các phép toán nhanh như dịch bit và trừ
- Khi đó, do phải xử lý các phép toán đắt đỏ hàng nghìn lần mỗi giây, cách tiếp cận này trở thành một thiết kế kỹ thuật phù hợp với giới hạn phần cứng thời bấy giờ
- Tuy nhiên, thuật toán này chỉ hoạt động với normal float
- Với giá trị sub-normal, giả định xấp xỉ
log2(1 + x) ≈ x + σ không còn đúng
- Trong trường hợp sub-normal, biểu thức thực tế gần với dạng
0 + x, nên phép xấp xỉ bị phá vỡ
Giảm sai số bằng hiệu chỉnh Newton-Raphson
- Giá trị khởi tạo thu được từ thao tác bit đã khá tốt, nhưng vẫn còn sai số có thể đo được
- Dòng sau giúp cải thiện phép xấp xỉ đáng kể
y = y * ( threehalfs - ( x2 * y * y ) );
- Dòng này là một dạng áp dụng phương pháp Newton-Raphson
- Để đưa bài toán nghịch đảo căn bậc hai về dạng phù hợp với phương pháp Newton, ta chuyển nó thành bài toán tìm nghiệm của hàm sau
f(y) = 1 / y^2 - x = 0
- Phương pháp Newton tạo giá trị xấp xỉ tốt hơn
y_(n+1) từ giá trị hiện tại y_n như sau
y_(n+1) = y_n - f(y_n) / f'(y_n)
- Đạo hàm của
f(y) = y^-2 - x là
f'(y) = -2y^-3 = -2 / y^3
Công thức hiệu chỉnh Newton không cần chia
- Nếu dùng nguyên công thức Newton, sẽ cần nhiều phép chia dấu phẩy động
- Một trong những lý do thuật toán này nhanh là tránh phép chia dấu phẩy động
- Sau khi biến đổi đại số, ta thu được dạng chỉ dùng phép nhân mà không cần chia
y_(n+1) = y_n * (1.5 - 0.5x * y_n^2)
- Trong mã Quake,
x2 = number * 0.5F được tính trước để lấy 0.5x, rồi dùng ở dòng sau
y = y * ( threehalfs - ( x2 * y * y ) );
- Sau một lần lặp này, sai số tuyệt đối tối đa là 0.175%, và trong nhiều trường hợp sai số còn thấp hơn
- Mã gốc có vòng lặp Newton thứ hai, nhưng đã bị comment out
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
Nguồn gốc và các thuật toán liên quan
- Thuật toán này không phải do John Carmack phát minh, và nguồn gốc chính xác của nó vẫn chưa được xác định chắc chắn 100%
- Có liên kết đến bài viết của Beyond3D về chủ đề này: The truth is the exact origin is not 100% certain
- Chris Lomont đã viết một bài báo nhằm tìm giá trị sigma tối ưu trong bước xấp xỉ log: InvSqrt.pdf
- CORDIC là thuật toán tính sine và cosine chỉ bằng phép cộng và dịch bit, không cần dấu phẩy động, và cách làm chi tiết của nó khá khác với fast inverse square root
- Điểm chung của hai thuật toán là đều áp dụng hiệu quả những quan sát toán học vào các ràng buộc phần cứng của thời kỳ đó
2 bình luận
Đoạn code thú vị lâu lâu lại được nhắc đến.. hehe
Ý kiến trên Hacker News
Nếu là máy tính được sản xuất từ sau năm 1999 thì nhìn chung sẽ hỗ trợ tập lệnh SSE, trong đó có
_mm_rsqrt_psđể tính 4 nghịch đảo căn bậc hai cùng lúc nhanh hơn: https://www.intel.com/content/www/us/en/docs/intrinsics-guid...Dù vậy, kỹ thuật được nói tới ở đây vẫn chưa hoàn toàn vô nghĩa. Việc chuyển đổi float/int thì nhanh, nhưng vẫn còn phần cứng không có lệnh
rsqrt,sqrt,pow,log, và các phép toán này có thể được xấp xỉ bằng mẹo nàysqrt(x)bằngx * 1/sqrt(x)sẽ nhanh hơnTrong các tập lệnh GPU, ARM, RISC-V, AVR, PIC, 8051, FPGA v.v., cũng thường có phép nghịch đảo căn bậc hai xấp xỉ được tích hợp sẵn, nhưng rất có thể chúng cũng được triển khai bằng kiểu thuật toán này
Nếu soi kỹ một chút thì phần giải thích kiểu như phép tính này ngày nay không diễn ra trên CPU là không đúng. Nhiều người vẫn hiểu lầm rằng game hay ứng dụng nặng về số thực dấu phẩy động sẽ muốn đẩy mọi phép toán số thực lên GPU
Trên thực tế, chỉ các tác vụ lớn và đồng đều mới đáng để chuyển sang GPU. Nếu là chuẩn hóa vector dùng một lần, như dựng ma trận quay để một vật thể hướng về vật thể khác, thì để lại trên CPU sẽ nhanh hơn. Ngay cả khi bỏ qua thời gian truyền sang GPU, với một phép toán số thực đơn lẻ thì CPU vẫn nhanh hơn, vì GPU thường có xung nhịp thấp hơn và đạt FLOP cao nhờ tính song song
Tôi đã thử viết một bản triển khai MMIX, với giả định giá trị đầu vào ban đầu lớn hơn
2^-1021Nếu quan tâm thì Wikipedia cũng có phần giải thích khá ổn về hàm này và lịch sử của nó: https://en.wikipedia.org/wiki/Fast_inverse_square_root
Tôi có gom vài thứ kiểu này ở đây: https://github.com/ncruces/fastmath/blob/main/fast.go
Cũng có một bài StackOverflow liên quan: https://stackoverflow.com/questions/32042673/optimized-low-a...
fastmathĐến giờ bắt lỗi rồi. Công thức float có lỗi đánh máy, phải là
(-1)^Schứ không phải-1^S. Cái sau lúc nào cũng ra-1Phần giải thích rằng diễn giải mẫu bit thô là một phép xấp xỉ tuyến tính từng đoạn của log cũng không thật chính xác. Những đoạn thẳng giữa các điểm dữ liệu trên đồ thị xanh thực ra không tồn tại, vì không thể có chuyện chỉ một nửa bit được đặt thành 1. Đúng hơn thì nó gần với một phiên bản rời rạc của log, và các điểm dữ liệu thực sự tồn tại — tức các giao điểm giữa đường đỏ và đường xanh — theo nghĩa đen chính là log đã được co giãn và tịnh tiến. Ngoài chỗ đó ra thì bài viết ổn
[010000, 010111]chứa 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75Nhưng phần trị do log cơ số 2 của các số này ngụ ý lần lượt là
.0000000,.0010101,.0101001,.0111010,.1001010,.1011001,.1100111,.1110100, và trừ cái đầu tiên ra thì chúng không trùng với001,010v.v. của float. Float trong đoạn[2,4)có khoảng cách tuyến tính, còn log tương ứng thì không, nên như bài viết nói, có thể xem float là một xấp xỉ tuyến tính từng đoạn của logNếu là toàn bộ đồ thị thì sẽ có
2^32lựa chọn trong mẫu tuyến tính từng đoạn, nhưng hình tác giả vẽ không phải toàn bộ đồ thị kiểu đó. Vì bài gốc nói về số nguyên 32 bit và phép toán float IEEE-754 32 bit, tôi nghĩ việc lược bỏ chữ “rời rạc” trong phần giải thích cũng chấp nhận đượcĐây là một bài viết hay, giải thích nhiều khái niệm thú vị, nhưng phần biến đổi đại số ở một mục thì tệ đến mức đáng kinh ngạc
Sau câu “Có nhiều bước chính xác để đi từ dạng đầu tiên sang dạng này, nhưng tôi đưa hết vào cho đầy đủ”, phần triển khai có quá nhiều bước không cần thiết, lại còn có vài lỗi dấu triệt tiêu lẫn nhau. Đặc biệt, khi đi từ dòng thứ hai sang dòng thứ ba, dấu âm không được phân phối đúng cách. Nếu bắt đầu từ sau dòng thứ hai thì có thể đi ngắn gọn hơn nhiều, từ
y_n+1 = y_n + (1 - x * y_n^2) / y_n^2 * (y_n^3 / 2)đếny_n+1 = y_n (1.5 * y_n - 0.5 * x * y_n * y_n), và các bước trung gian cũng sẽ đúng. Với người hiểu đại số thì đó đều là những bước hiển nhiênMagic number trong đoạn mã nổi tiếng đó không phải là hằng số tối ưu. Dùng hằng số khác có lẽ có thể giảm sai số tương đối thêm khoảng 0.5%
Có thể hồi đó khó tìm ra giá trị tối ưu tuyệt đối, nhưng bây giờ thì tương đối dễ. Tôi cũng từng chui vào cái hang thỏ này và có một notebook Jupyter để tìm magic number tối ưu cho
(1/x^2)và(1/x)Điều thú vị nhất với tôi trong bài này lại là link “How Java's Floating-Point Hurts Everyone Everywhere”: https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf
Tác giả là William Kahan, cũng được biết đến như “Old Man of Floating-Point”: https://news.ycombinator.com/item?id=29042853 - An Interview with the Old Man of Floating-Point (1998)
JAVAhurtthì thấy dàn trang kinh khủng. Cảm giác như họ dùng một gói TeX kéo giãn khoảng cách giữa các từ quá mức, lại còn không đều, hoặc OCR từ tài liệu khác rồi sinh ra các khoảng trắng thừaNgay cả phần dùng font monospace cũng có thêm những khoảng cách lạ. Rất khó tập trung để đọc, và dù tôi biết thực ra không phải vậy, nó gần như tạo cảm giác như đang đọc một bản tuyên ngôn của dân khoa học quái kiệt
Tôi từng xem video này trước đây và nó thực sự rất hay: https://www.youtube.com/watch?v=p8u_k2LIZyo