Akar Kuadrat dan 5f3759df

Saat source code untuk game Quake III dibuka, orang menemukan kode-kode C menarik dari John Carmack. Salah satunya adalah fungsi invers akar kuadrat, yang pada intinya ditulis sebagai berikut:

float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0x5f3759df- (i>>1);
x = *(float*)&i;
x = x*(1.5f-xhalf*x*x);
return x;
}

Menariknya, hack semacam ini menghasilkan kecepatan kalkulasi yang amat cepat, kira-kira empat kali lebih cepat daripada menggunakan (float)(1.0/sqrt(x)), walaupun sqrt dalam hal ini menggunakan instruksi assembly FSQRT.

Yang jelas, kode ini memanfaatkan metode Newton-Raphson. Tetapi pendekatan yang digunakan, dan terutama penggunaan heksadesimal 5f3759df tentu sangat menarik. Googlekan angka itu, dan temukan petualangan menarik mencari asal-usul baik bilangan ajaib itu, maupun hacknya sendiri. D Ebery misalnya, menganggap shift i>>1 mengakibatkan interpolasi linear pada si invers akar kuadrat. Tetapi, pertama kali, kita harus ingat bahwa si penulis kode sedang memainkan bit-bit floating point yang secara standar akan dikodekan sesuai IEEE 754-1985 (yang memisahkan mantissa dengan eksponen); dan si penulis juga meyakinkan diri bahwa kode ini jalan baik untuk little endian maupun big endian.

Sekarang, metode Newton-Raphson dulu, sambil dipandu C Lomont. Kita akan menghitung 1/akar(x). Definisikan dulu f(y)=1/y² – x, sehingga nantinya nilai yang kita cari adalah akar positif dari f(x). Dengan metode Newton, jika kita punya pendekatan awal y_n, maka kitadapat menghitung pendekatan berikutnya sebagai y_n+1 = y_n – f(y_n)/f'(y_n). Dengan f(y) di atas, akan diperoleh y_n+1 = ½y_n (3 – xy_n²), atau dalam kode C di atas adalah x = x*(1.5f-xhalf*x*x), dengan x nilai pendekatan awal y₀ kita. Baris i = 0x5f3759df-(i>>1) menghitung nilai awal y₀ ini, dengan mengalikan eksponen x dengan -½. Kemudian bagian yang menarik pun dimulai. Karena, hey, yang dishift kan bukan hanya eksponen, tetapi keseluruhan bilangan.

Ada beberapa pendekatan untuk bagian ini. Aku baca baik versi C Lomont maupun C McEniry. Lomont membagi ulasan untuk eksponen genap dan ganjil segala. Tapi akhirnya yang mereka dapati adalah sebuah pola berulang, yang kemudian dicari minimasi kesalahannya, sehingga diperoleh sebuah nilai dengan kesalahan paling kecil. Dan kode di atas itulah hasilnya, dengan sebuah kiraan r 0.4327448899 dst. Lomont sendiri mendapati bahwa ia memperoleh bilangan yang memberikan akurasi lebih tinggi, yaitu 5f375a86. McEniry mengkritik bahwa brutal force yang sempat digunakan Lomont bisa justru tidak mendeteksi jurang sempit antar celah serangan brutal itu. Tapi, sebagaimana Lomont, ia juga memberikan saran perbaikan, baik untuk versi float maupun versi double.

float InvSqrt(float x)
{
union {float f; unsigned long ul; } y;
y.f = x;
y.ul = ( 0xBE6EB50C – y.ul ) >> 1;
y.f = 0.5f * y.f * (3.0f – x * y.f * y.f);
return y.f;
}

Untuk versi double, keyword float harus diganti double, unsigned long menjadi unsigned long long, dan si konstanta ajaib jadi 0xBFCDD6A18F6A6F54.