Semua Pengetahuan tentang Algoritme Inverse Square Root Cepat

(github.com/francisrstokes)

5 poin oleh GN⁺ 2024-06-03 | 2 komentar | Bagikan ke WhatsApp

fast inverse square root yang terkenal lewat Quake 3 adalah solusi performa pada masanya untuk mengaproksimasi 1 / sqrt(x) dengan reinterpretasi bit float dan koreksi Newton-Raphson
Intinya, pola bit bilangan float 32-bit IEEE-754 dapat diperlakukan seperti aproksimasi log2(x) yang diskalakan dan digeser
0x5f3759df - (i >> 1) adalah bentuk yang memindahkan log2(x^-0.5) = -0.5 * log2(x) ke dalam shift dan pengurangan integer, dan konstanta ajaibnya berasal dari 3/2 * 2^23 * (127 - σ)
Setelah itu, koreksi Newton-Raphson diterapkan sekali dengan y = y * (1.5 - 0.5x * y * y), dan iterasi kedua pada kode Quake diberi komentar
Pada 1999, inverse square root dibutuhkan ratusan hingga ribuan kali per detik untuk pencahayaan dan normalisasi vektor 3D, tetapi pada perangkat keras modern kegunaan praktis trik ini berkurang berkat unit floating-point khusus

Apa yang dilakukan kode Quake

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

Fungsi ini menghitung nilai aproksimasi inverse square root 1 / sqrt(number) untuk number
Bagian paling terkenalnya adalah manipulasi bit yang menafsirkan nilai float sebagai long, lalu menjalankan 0x5f3759df - (i >> 1)
Saat Quake 3 dirilis pada 1999, inverse square root adalah operasi yang lambat dan mahal, tetapi dibutuhkan ratusan hingga ribuan kali per detik dalam persamaan pencahayaan dan perhitungan vektor 3D yang perlu dinormalisasi
Pada perangkat keras modern, perhitungan seperti ini sering tidak dikerjakan di CPU, atau jika dijalankan di CPU pun tetap cepat berkat perangkat keras floating-point khusus yang jauh lebih maju

Representasi float 32-bit IEEE-754

Float 32-bit terdiri dari tiga bagian
- Sign: 1 bit, menunjukkan apakah nilainya positif atau negatif
- Exponent: 8 bit, menentukan rentang tempat nilai berada
- Mantissa: 23 bit, menunjukkan posisi secara linear di dalam rentang tersebut
Nilai normal diinterpretasikan dalam bentuk berikut

N = (-1)^S * 2^(E - 127) * (1 + M / 2^23)

B = 127 adalah nilai bias yang dipakai untuk biased exponent, dan eksponen sebenarnya adalah e = E - B
Mantissa tidak ditulis hanya sebagai m, melainkan dalam bentuk 1 + m
- Jika m = 0, nilainya adalah 2^e
- Saat m mendekati 1, nilainya merepresentasikan hingga tepat sebelum rentang eksponen berikutnya, yaitu 2^(e+1)
Jika semua bit eksponen bernilai 0, itu adalah bilangan sub-normal, dan rumusnya berubah

N = (-1)^S * 2^-126 * m

Sub-normal diperlukan untuk merepresentasikan 0 dan bilangan sangat kecil yang sangat dekat dengan 0
Jika semua bit eksponen bernilai 1, nilainya diperlakukan sebagai nilai khusus
- Jika E = 255 dan M = 0, hasilnya Infinity atau -Infinity
- Jika M != 0, hasilnya NaN

Relasi logaritmik saat bit float dilihat sebagai integer

Jika representasi internal float dilihat sebagai integer 32-bit, bentuknya dapat ditulis sebagai berikut

I_x = 2^31 S + 2^23 E + M

Karena inverse square root bekerja pada input positif, kita dapat mengambil S = 0 sehingga rumusnya menjadi lebih sederhana

L = 2^23
I_x = L E + M

Di dalam rentang exponent yang sama, mantissa menunjukkan posisi secara linear, tetapi ketika exponent membesar, jumlah langkah mantissa yang sama mencakup segmen garis bilangan yang lebih lebar
- E = 127, yaitu e = 0, kira-kira mencakup rentang [1, 2)
- E = 128, yaitu e = 1, kira-kira mencakup rentang [2, 4)
- Kedua rentang memiliki jumlah langkah mantissa yang sama, tetapi rentang kedua dua kali lebih lebar
Karena struktur ini, pola bit mentah float bila dilihat sebagai integer menunjukkan relasi logaritmik

Bit mentah adalah aproksimasi `log2(x)`

Jika pola bit float ditafsirkan sebagai integer I_x, ia bisa dipandang sebagai aproksimasi linear per segmen terhadap log2(x)
Relasi ini dapat dinyatakan dengan pendekatan berikut

log2(x) ≈ I_x / L - B

Jika integer bit mentah dibagi dengan ukuran mantissa L = 2^23 lalu dikurangi bias exponent B = 127, hasilnya menjadi nilai yang dekat dengan log2(x)
Logaritma di dalam rentang mantissa diperlakukan sebagai aproksimasi linear

log2(1 + x) ≈ x + σ

σ adalah parameter tuning untuk menyesuaikan aproksimasi, dan x menunjukkan posisi di dalam rentang exponent dengan domain [0, 1]

Mengubah inverse square root menjadi identitas logaritma

Tujuannya adalah menghitung nilai berikut

y = 1 / sqrt(x)

Ini dapat ditulis ulang dalam bentuk eksponen sebagai berikut

y = x^-0.5

Jika identitas logaritma diterapkan, perhitungan inverse square root menjadi relasi berikut

log2(1 / sqrt(x)) = log2(x^-0.5) = -0.5 * log2(x)

Dengan memanfaatkan fakta bahwa bit float berperilaku seperti aproksimasi log2(x), kita bisa langsung mengaproksimasi representasi bit integer I_y dari y berdasarkan representasi bit integer I_x dari x

I_y ≈ -0.5 I_x + 1.5 L (B - σ)

Rumus ini mengarah langsung ke satu baris inti pada kode Quake

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

i >> 1 berperan seperti mengalikan bit integer dengan 1/2 melalui shift kanan 1 bit
Konstanta di depannya, 0x5f3759df, sesuai dengan 1.5 * L * (B - σ)

Identitas konstanta `0x5f3759df`

Jika σ = 0, konstanta tersebut dihitung sebagai berikut

1.5 * 2^23 * 127 = 1598029824

Representasi heksadesimal nilai itu adalah 0x5f400000
Nilai ini berbeda 566817 dari konstanta nyata Quake, 0x5f3759df
Dari selisih itu, nilai σ yang sesuai dengan kode Quake dapat dihitung sebagai berikut

σ = 377878 / 2^23
σ = 0.04504656

Dalam C, konstanta yang sama bisa dihitung seperti ini

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

Di sini digunakan double, dan konversi ke integer adalah casting biasa, bukan reinterpretasi bit
Nilai σ ini dipilih untuk mengoptimalkan aproksimasi, tetapi bukan nilai optimal sebenarnya, dan juga tidak pasti siapa yang pertama membuatnya

Mengapa ini bukan sekadar hack sederhana

0x5f3759df - (i >> 1) adalah rumus untuk membuat nilai awal inverse square root dengan memanfaatkan fakta bahwa bit mentah float adalah aproksimasi logaritma
Meskipun didasarkan pada relasi matematis yang kompleks, pada tahap eksekusi ia hanya memakai operasi cepat seperti shift dan pengurangan
Pada masa itu, operasi mahal harus diproses ribuan kali per detik, sehingga pendekatan ini menjadi rancangan rekayasa yang sesuai dengan keterbatasan perangkat keras saat itu
Namun, algoritme ini hanya bekerja pada normal float
- Untuk nilai sub-normal, asumsi aproksimasi log2(1 + x) ≈ x + σ tidak berlaku
- Pada sub-normal, bentuk yang masuk dalam praktik lebih dekat ke 0 + x, sehingga aproksimasinya gagal

Mengurangi galat dengan koreksi Newton-Raphson

Nilai awal yang diperoleh dari manipulasi bit sudah cukup baik, tetapi masih menyisakan galat yang terukur
Baris berikut sangat meningkatkan aproksimasi

y = y * ( threehalfs - ( x2 * y * y ) );

Baris ini adalah bentuk penerapan metode Newton-Raphson
Agar masalah inverse square root cocok dengan metode Newton, ia diubah menjadi masalah mencari akar dari fungsi berikut

f(y) = 1 / y^2 - x = 0

Metode Newton menghasilkan aproksimasi yang lebih baik dari nilai saat ini y_n menjadi y_(n+1) sebagai berikut

y_(n+1) = y_n - f(y_n) / f'(y_n)

Turunan dari f(y) = y^-2 - x adalah sebagai berikut

f'(y) = -2y^-3 = -2 / y^3

Rumus koreksi Newton tanpa pembagian

Jika rumus Newton dipakai apa adanya, akan ada beberapa operasi pembagian floating-point
Salah satu alasan algoritme ini cepat adalah menghindari pembagian floating-point
Jika disusun ulang secara aljabar, bentuknya menjadi hanya menggunakan perkalian tanpa pembagian

y_(n+1) = y_n * (1.5 - 0.5x * y_n^2)

Dalam kode Quake, x2 = number * 0.5F menghitung 0.5x lebih dahulu, lalu dipakai pada baris berikut

y = y * ( threehalfs - ( x2 * y * y ) );

Setelah satu kali iterasi ini, galat absolut maksimum adalah 0.175%, dan dalam banyak kasus galatnya lebih kecil dari itu
Kode asli juga memiliki iterasi Newton kedua, tetapi diberi komentar

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

Asal-usul dan algoritme terkait

Algoritme ini bukan ciptaan John Carmack, dan asal-usul pastinya tidak 100% diketahui
Ada tulisan dari Beyond3D yang ditautkan terkait hal ini: The truth is the exact origin is not 100% certain
Chris Lomont menulis makalah yang mencoba mencari nilai sigma optimal pada tahap aproksimasi logaritma: InvSqrt.pdf
CORDIC adalah algoritme untuk menghitung sine dan cosine hanya dengan penjumlahan dan bit shift tanpa floating-point, dan detail caranya sangat berbeda dari fast inverse square root
Kedua algoritme itu memiliki kesamaan dalam hal menerapkan pengamatan matematis secara efisien sesuai keterbatasan perangkat keras pada zamannya

2 komentar

joyfui 2024-06-03

Kode unik yang suka muncul lagi tiap kali hampir terlupakan.. hehe

GN⁺ 2024-06-03

Komentar Hacker News

Jika itu komputer yang dibuat setelah 1999, kemungkinan besar mendukung set instruksi SSE, dan di dalamnya ada _mm_rsqrt_ps yang dapat menghitung 4 inverse square root sekaligus dengan lebih cepat: https://www.intel.com/content/www/us/en/docs/intrinsics-guid...
Meski begitu, teknik yang dibahas di sini masih belum sepenuhnya tidak berguna. Konversi float/int memang cepat, tetapi masih ada hardware yang tidak memiliki instruksi rsqrt, sqrt, pow, atau log, dan operasi-operasi ini bisa didekati dengan trik ini
- Instruksi reciprocal floating-point di SSE dapat memberikan hasil yang sedikit berbeda di Intel dan AMD, jadi jika mengharapkan hasil deterministik antar-PC, ini bisa merepotkan: https://robert.ocallahan.org/2021/09/rr-trace-portability-di...
- Menariknya, SSE juga punya instruksi square root biasa, tetapi jauh lebih lambat daripada inverse square root, jadi jika bisa menerima penurunan presisi, menghitung sqrt(x) sebagai x * 1/sqrt(x) lebih cepat
- Sebenarnya, mayoritas komputer yang sangat besar bahkan tidak mendukung set instruksi i386/amd64, apalagi SSE, dan proporsi yang mendukungnya secara native, bukan lewat emulasi, terus menurun
  Pada GPU ISA, ARM, RISC-V, AVR, PIC, 8051, FPGA, dan sebagainya, sering kali ada operasi approximate inverse square root bawaan, tetapi kemungkinan besar diimplementasikan dengan algoritme semacam ini
Sedikit mencari celah pada artikelnya, penjelasan seolah-olah perhitungan seperti ini tidak terjadi di CPU modern itu tidak tepat. Ini terkait salah paham umum bahwa game atau aplikasi dengan banyak operasi floating-point ingin memindahkan semua operasi floating-point ke GPU
Dalam praktiknya, yang layak dipindahkan ke GPU hanya pekerjaan besar dan seragam. Jika Anda melakukan normalisasi vektor sekali pakai, seperti menyusun matriks rotasi agar satu objek menghadap objek lain, lebih cepat membiarkannya tetap di CPU. Bahkan tanpa menghitung waktu transfer ke GPU, operasi floating-point tunggal lebih cepat di CPU, karena GPU biasanya ber-clock lebih rendah dan mencapai FLOP tinggi lewat paralelisme
- Sepertinya yang dimaksud di sini bukan GPU melainkan FPU. Dulu FPU melakukan komputasi asinkron, dan sekarang dianggap sebagai bagian yang terintegrasi ke dalam CPU
Saya sempat menulis implementasi MMIX, dengan asumsi nilai input awal lebih besar dari 2^-1021
Kalau tertarik, Wikipedia juga punya penjelasan yang cukup bagus tentang fungsi ini dan sejarahnya: https://en.wikipedia.org/wiki/Fast_inverse_square_root
Saya mengumpulkan beberapa hal seperti ini: https://github.com/ncruces/fastmath/blob/main/fast.go
Ada juga postingan StackOverflow terkait: https://stackoverflow.com/questions/32042673/optimized-low-a...
- Berguna, karena saya memang sedang berpikir untuk mulai mengumpulkan teknik seperti ini agar bisa memakai lagi mesin 3D gaya akhir 80-an yang pernah saya buat
- Saya juga ingin melihat benchmark dari paket fastmath
Waktunya sedikit mengoreksi. Ada salah ketik pada rumus float, dan seharusnya (-1)^S, bukan -1^S. Yang terakhir itu selalu bernilai -1
Penjelasan bahwa menafsirkan pola bit mentah merupakan aproksimasi linear per-bagian dari log juga tidak akurat. Garis di antara titik data pada grafik biru sebenarnya tidak ada, dan tidak mungkin hanya setengah bit yang disetel ke 1. Ini lebih dekat ke versi diskret dari log, dan titik data yang benar-benar ada—yakni titik pertemuan garis merah dan biru—secara harfiah sama dengan log yang diskalakan dan digeser. Selain itu, artikelnya bagus
- Saya kurang paham. Jika kita membayangkan float 6-bit yang sangat kecil, dengan 1 bit tanda, 2 bit eksponen, dan 3 bit mantissa, maka rentang [010000, 010111] memuat 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75
  Namun mantissa yang diimplikasikan oleh log basis-2 dari angka-angka ini adalah .0000000, .0010101, .0101001, .0111010, .1001010, .1011001, .1100111, .1110100, dan selain yang pertama, nilainya tidak sama dengan 001, 010, dan seterusnya pada float. Karena float dalam rentang [2,4) berjarak linear sedangkan log yang bersesuaian tidak, maka seperti yang dikatakan artikel, float bisa dipandang sebagai aproksimasi linear per-bagian dari log
- Bukan aproksimasi linear per-bagian yang kontinu, melainkan aproksimasi linear per-bagian yang diskret. Memang benar garis biru itu tidak kontinu, tetapi interpretasinya keliru. Grafik biru terdiri dari 256 titik terpisah yang tersebar merata di sepanjang sumbu x, bukan hanya beberapa titik perpotongan
  Jika itu seluruh grafik, akan ada 2^32 kemungkinan di dalam pola per-bagian linear tersebut, tetapi yang digambar tulisan aslinya bukan seluruh grafik seperti itu. Karena tulisan tersebut membahas operasi integer 32-bit dan float IEEE-754 32-bit, menurut saya menghilangkan kata “diskret” dalam penjelasannya masih bisa diterima
Tulisan ini bagus dan menjelaskan banyak konsep menarik, tetapi penurunan aljabar di salah satu bagiannya sangat buruk
Setelah “Ada banyak langkah tepat untuk berpindah dari bentuk pertama ke bentuk ini, tetapi saya sertakan semuanya demi kelengkapan”, penurunannya memuat banyak langkah yang tidak perlu dan juga beberapa kesalahan tanda yang saling meniadakan. Khususnya, saat berpindah dari baris kedua ke baris ketiga, tanda negatif tidak didistribusikan dengan benar. Setelah baris kedua, sebenarnya bisa mulai dari y_n+1 = y_n + (1 - x * y_n^2) / y_n^2 * (y_n^3 / 2) lalu langsung ke y_n+1 = y_n (1.5 * y_n - 0.5 * x * y_n * y_n) dengan jauh lebih singkat, dan langkah-langkah di tengah pun benar. Menurut saya itu semua hanyalah langkah yang jelas bagi orang yang paham aljabar
Magic number pada potongan kode terkenal itu bukan konstanta optimal. Dengan konstanta lain, galat relatifnya mungkin bisa dikurangi sekitar 0,5% lagi
Dulu mungkin sulit mencari nilai paling optimal, tetapi sekarang relatif mudah. Saya sendiri sempat terjun ke lubang kelinci ini dan punya notebook Jupyter untuk mencari magic number optimal bagi (1/x^2) dan (1/x)
- Tautan ke paper yang membahas pertanyaan itu ada di bagian bawah tulisan
Hal paling menarik dari tulisan ini bagi saya justru tautan “How Java's Floating-Point Hurts Everyone Everywhere”: https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf
Penulisnya adalah William Kahan, yang juga dikenal sebagai “Old Man of Floating-Point”: https://news.ycombinator.com/item?id=29042853 - An Interview with the Old Man of Floating-Point (1998)
- Tidak terkait langsung dengan topiknya, tetapi saya mulai membaca PDF JAVAhurt dan tata letaknya mengerikan. Rasanya seperti memakai paket TeX yang memberi jarak antar-kata berlebihan dan tidak merata, dan seolah-olah itu hasil OCR dari dokumen lain yang lalu menyisakan spasi ekstra
  Bahkan di bagian font monospace pun ada jarak tambahan yang aneh. Sangat sulit fokus saat membacanya, dan meskipun jelas bukan begitu, kesannya hampir seperti manifesto ilmuwan eksentrik
Saya sangat menyukai video ini yang pernah saya tonton dulu: https://www.youtube.com/watch?v=p8u_k2LIZyo

Semua Pengetahuan tentang Algoritme Inverse Square Root Cepat

Apa yang dilakukan kode Quake

Representasi float 32-bit IEEE-754

Relasi logaritmik saat bit float dilihat sebagai integer

Bit mentah adalah aproksimasi log2(x)

Mengubah inverse square root menjadi identitas logaritma

Identitas konstanta 0x5f3759df

Mengapa ini bukan sekadar hack sederhana

Mengurangi galat dengan koreksi Newton-Raphson

Rumus koreksi Newton tanpa pembagian

Asal-usul dan algoritme terkait

Bacaan terkait

2 komentar

Komentar Hacker News

Bit mentah adalah aproksimasi `log2(x)`

Identitas konstanta `0x5f3759df`