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/y2 – x, sehingga nantinya nilai yang kita cari adalah akar positif dari f(x). Dengan metode Newton, jika kita punya pendekatan awal yn, maka kitadapat menghitung pendekatan berikutnya sebagai yn+1 = yn – f(yn)/f'(yn). Dengan f(y) di atas, akan diperoleh yn+1 = ½yn (3 – xyn2), atau dalam kode C di atas adalah x = x*(1.5f-xhalf*x*x), dengan x nilai pendekatan awal y0 kita. Baris i = 0x5f3759df-(i>>1) menghitung nilai awal y0 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.

3 Comments

  1. Menggugah ingin saya u/ belajar. Ga ngerti !! tapi menarik. Fiuh tapi umur sudah segini. terlambat kah?
    Sepertinya bukan John Carmack author fungsi di atas.

    http://www.beyond3d.com/content/articles/8/

    Regards,

  2. soidera

    buset barus sadar kalo calculuss cuma C itu buruk >_< , pusing bgt. tapi emang keren hacker2 sejati

  3. hi, thanks for the story – very inspirationalVery useful site. I’ll be recomend this blog for my friends.?

Leave a Reply

%d bloggers like this: