Hızlı ters karekök algoritması hakkında bilinmesi gereken her şey

(github.com/francisrstokes)

5 puan yazan GN⁺ 2024-06-03 | 2 yorum | WhatsApp'ta paylaş

Quake 3 ile ünlenen fast inverse square root, 1 / sqrt(x) değerini float bitlerini yeniden yorumlama ve Newton-Raphson düzeltmesiyle hızlıca yaklaşık hesaplayan, döneminin performans çözümüydü
Temel nokta, IEEE-754 32 bit float’ın tamsayı bit deseninin ölçeklenmiş ve kaydırılmış bir log2(x) yaklaşımı gibi ele alınabilmesidir
0x5f3759df - (i >> 1), log2(x^-0.5) = -0.5 * log2(x) ifadesinin tamsayı kaydırma ve çıkarma işlemine aktarılmış hâlidir; sihirli sabit 3/2 * 2^23 * (127 - σ) ifadesinden gelir
Ardından y = y * (1.5 - 0.5x * y * y) bir kez uygulanarak Newton-Raphson düzeltmesi yapılır; Quake kodundaki ikinci yineleme yorum satırına alınmıştır
1999’da aydınlatma ve 3B vektör normalizasyonu için ters karekök saniyede yüzlerce ila binlerce kez gerekliydi; modern donanımlarda ise özel kayan nokta işleme sayesinde aynı numaranın pratik değeri azalmıştır

Quake kodu ne yapıyor

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

Bu fonksiyon, number için ters karekök 1 / sqrt(number) değerinin yaklaşık sonucunu hesaplar
En ünlü kısmı, float değeri long gibi yorumladıktan sonra 0x5f3759df - (i >> 1) yapan bit manipülasyonudur
Quake 3’ün çıktığı 1999’da ters karekök yavaş ve pahalı bir işlemdi; aydınlatma denklemleri ve normalizasyon gerektiren 3B vektör hesaplarında saniyede yüzlerce ila binlerce kez gerekiyordu
Modern donanımlarda bu tür hesaplamalar ya CPU’da yapılmaz ya da CPU’da çalıştırılsa bile gelişmiş özel kayan nokta donanımı sayesinde hızlıdır

IEEE-754 32 bit float gösterimi

32 bit float üç bölümden oluşur
- Sign: 1 bit; pozitif/negatif olup olmadığını gösterir
- Exponent: 8 bit; değerin ait olduğu aralığı belirler
- Mantissa: 23 bit; bu aralık içindeki konumu doğrusal olarak gösterir
Normal bir değer şu biçimde yorumlanır

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

B = 127, biased exponent için kullanılan bias değeridir; gerçek üs e = E - B olur
Mantissa basitçe m ile çarpılmaz, 1 + m biçiminde kullanılır
- m = 0 ise 2^e
- m = 1 değerine yaklaştığında, bir sonraki üs aralığı olan 2^(e+1) değerinin hemen öncesine kadar ifade eder
Üsün tüm bitleri 0 ise sayı sub-normal’dır ve formül değişir

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

Sub-normal, 0’ı ve 0’a çok yakın küçük sayıları temsil etmek için gereklidir
Üsün tüm bitleri 1 ise özel değer olarak işlenir
- E = 255, M = 0 ise Infinity veya -Infinity
- M != 0 ise NaN

Float bitlerine tamsayı olarak bakınca ortaya çıkan logaritmik ilişki

Float’ın iç gösterimine 32 bit tamsayı gibi bakarsak şu formülle ifade edilebilir

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

Ters karekök pozitif girdileri hedeflediği için S = 0 alınırsa formül sadeleşir

L = 2^23
I_x = L E + M

Aynı exponent aralığında mantissa konumu doğrusal olarak gösterir; ancak exponent büyüdükçe aynı sayıdaki mantissa adımı sayı doğrusu üzerinde daha geniş bir aralığı kapsar
- E = 127, yani e = 0, yaklaşık [1, 2) aralığıdır
- E = 128, yani e = 1, yaklaşık [2, 4) aralığıdır
- İki aralıkta da mantissa adımı sayısı aynıdır, ancak ikinci aralık iki kat geniştir
Bu yapı nedeniyle float’ın ham bit desenine tamsayı olarak bakıldığında logaritmik bir ilişki ortaya çıkar

Ham bitler `log2(x)` için bir yaklaşımdır

Float’ın bit desenini tamsayı I_x olarak yorumlarsak, bunu log2(x) için parçalı doğrusal bir yaklaşım gibi görebiliriz
Bu ilişki şu yaklaşık formülle ifade edilir

log2(x) ≈ I_x / L - B

Ham bit tamsayısını mantissa boyutu L = 2^23 değerine bölüp exponent bias’ı B = 127 çıkardığımızda log2(x) değerine yakın bir sonuç elde edilir
Mantissa aralığındaki logaritma, doğrusal yaklaşımla ele alınır

log2(1 + x) ≈ x + σ

σ, yaklaşımı ayarlayan bir tuning parametresidir; x, [0, 1] aralığında exponent bölgesi içindeki konumu gösterir

Ters karekökü logaritma özdeşliğine çevirmek

Amaç şu değeri bulmaktır

y = 1 / sqrt(x)

Bu, üs biçiminde şöyle yazılır

y = x^-0.5

Logaritma özdeşliği uygulanınca ters karekök hesabı şu ilişkiye dönüşür

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

Float bitlerinin log2(x) yaklaşımı gibi davranmasından yararlanarak, x’in tamsayı bit gösterimi I_x üzerinden y’nin tamsayı bit gösterimi I_y doğrudan yaklaşık hesaplanabilir

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

Bu formül Quake kodunun kritik tek satırına dönüşür

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

i >> 1, tamsayı bitlerini 1 bit sağa kaydırarak 1/2 ile çarpma işlevi görür
Baştaki 0x5f3759df sabiti 1.5 * L * (B - σ) ifadesine karşılık gelir

`0x5f3759df` sabitinin gerçek kimliği

σ = 0 alınırsa sabit şöyle hesaplanır

1.5 * 2^23 * 127 = 1598029824

Bu değerin onaltılık gösterimi 0x5f400000 olur
Quake’in gerçek sabiti 0x5f3759df ile arasında 566817 fark vardır
Bu farktan Quake koduna karşılık gelen σ değeri şöyle hesaplanır

σ = 377878 / 2^23
σ = 0.04504656

C koduyla aynı sabit şu şekilde hesaplanabilir

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

Burada double kullanılır; tamsayıya dönüştürme bitlerin yeniden yorumlanması değil, normal casting’dir
Bu σ değeri yaklaşımı optimize etmek için seçilmiştir; ancak gerçek optimum değer değildir ve kimin oluşturduğu da kesin değildir

Basit bir hack olmamasının nedeni

0x5f3759df - (i >> 1), float’ın ham bitlerinin logaritma yaklaşımı olmasından yararlanarak ters karekök başlangıç değeri üreten bir formüldür
Karmaşık bir matematiksel ilişkiye dayanır, ancak yürütme aşamasında yalnızca kaydırma ve çıkarma gibi hızlı işlemler kullanır
O dönemde pahalı işlemleri saniyede binlerce kez yapmak gerektiğinden, bu yöntem donanım kısıtlarına uygun bir mühendislik tasarımıydı
Ancak bu algoritma yalnızca normal float değerlerde çalışır
- Sub-normal değerlerde log2(1 + x) ≈ x + σ yaklaşımının varsayımı geçerli değildir
- Sub-normal’de gerçekte 0 + x değerine yakın bir biçim bulunduğundan yaklaşım bozulur

Newton-Raphson düzeltmesiyle hatayı azaltmak

Bit manipülasyonuyla elde edilen başlangıç değeri oldukça iyidir, ancak ölçülebilir bir hata kalır
Şu satır yaklaşımı büyük ölçüde iyileştirir

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

Bu satır Newton-Raphson method uygulanmış hâlidir
Ters karekök problemini Newton method’a uygun hâle getirmek için şu fonksiyonun kökünü bulma problemine dönüştürür

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

Newton method, mevcut yaklaşık değer y_n üzerinden daha iyi bir yaklaşık değer y_(n+1) değerini şöyle üretir

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

f(y) = y^-2 - x fonksiyonunun türevi şöyledir

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

Bölmesiz Newton düzeltme formülü

Newton formülünü olduğu gibi kullanırsak birden fazla kayan nokta bölmesi gerekir
Bu algoritmanın hızlı olmasının nedenlerinden biri kayan nokta bölmesinden kaçınmasıdır
Cebirsel olarak düzenlenince bölme olmadan, yalnızca çarpma kullanan bir biçime dönüşür

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

Quake kodunda x2 = number * 0.5F ile 0.5x önceden hesaplanır ve şu satırda kullanılır

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

Bu tek yinelemeden sonra en büyük mutlak hata %0,175 olur; birçok durumda hata bundan daha düşüktür
Orijinal kodda ikinci bir Newton yinelemesi vardır ancak yorum satırına alınmıştır

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

Kökeni ve ilgili algoritmalar

Bu algoritmayı John Carmack icat etmemiştir; kesin kökeni %100 belli değildir
İlgili olarak Beyond3D’deki şu yazıya bağlantı verilmiştir: The truth is the exact origin is not 100% certain
Chris Lomont, logaritma yaklaşımı aşamasında optimal sigma değerini bulmaya çalışan bir makale yazmıştır: InvSqrt.pdf
CORDIC, kayan nokta olmadan yalnızca toplama ve bit kaydırma ile sine ve cosine hesaplayan bir algoritmadır; fast inverse square root ile ayrıntıdaki yöntemi oldukça farklıdır
İki algoritmanın ortak yanı, matematiksel gözlemleri dönemin donanım kısıtlarına uygun şekilde verimli olarak uygulamalarıdır

2 yorum

joyfui 2024-06-03

İnsan unuturken bir anda yeniden ortaya çıkan o ilginç kod.. heh

GN⁺ 2024-06-03

Hacker News yorumları

1999’dan sonra üretilmiş bir bilgisayar büyük olasılıkla SSE komut kümesini destekler; bunun içinde tek seferde 4 ters karekökü daha hızlı hesaplayan _mm_rsqrt_ps de var: https://www.intel.com/content/www/us/en/docs/intrinsics-guid...
Yine de burada ele alınan teknik hâlâ tamamen anlamsız değil. float/int dönüşümü hızlı olsa da rsqrt, sqrt, pow, log komutları olmayan donanımlar hâlâ var; bu tür işlemler bu numarayla yaklaşık olarak hesaplanabiliyor
- SSE’nin kayan nokta reciprocal komutu Intel ve AMD’de biraz farklı sonuçlar verebilir; PC’ler arasında deterministik sonuç bekliyorsanız başınızı ağrıtabilir: https://robert.ocallahan.org/2021/09/rr-trace-portability-di...
- İlginçtir, SSE’de normal karekök komutu da var ama ters karekökten çok daha yavaş; hassasiyet kaybını göze alabiliyorsanız sqrt(x) değerini x * 1/sqrt(x) olarak hesaplamak daha hızlı
- Aslında bilgisayarların ezici çoğunluğu SSE bir yana, i386/amd64 komut kümesinin kendisini bile desteklemiyor; bunu emülasyon olmadan destekleyenlerin oranı da giderek azalıyor
  GPU komut kümeleri, ARM, RISC-V, AVR, PIC, 8051, FPGA vb. üzerinde yaklaşık ters karekök işlemi çoğu zaman yerleşik olabiliyor; ama muhtemelen bunun da bu tür algoritmalarla uygulanmış olması olası
Yazıya biraz itiraz edecek olursam, bu tür hesaplamaların günümüz CPU’larında yapılmadığı yönündeki açıklama doğru değil. Oyunların ya da kayan nokta işlemi yoğun uygulamaların tüm kayan nokta işlemlerini GPU’ya devretmek istediği yaygın bir yanılgı
Gerçekte GPU’ya devretmenin mantıklı olduğu işler yalnızca büyük ve tekdüze işlerdir. Bir nesnenin başka bir nesneye bakmasını sağlamak için dönüş matrisi oluşturmak gibi tek seferlik vektör normalizasyonu yapıyorsanız, bunu CPU’da bırakmak daha hızlıdır. GPU’ya aktarım süresini hariç tutsanız bile tekil kayan nokta işlemlerinde CPU daha hızlıdır; çünkü GPU’lar genellikle daha düşük saat hızında çalışır ve yüksek FLOP sayılarına paralellik sayesinde ulaşır
- Burada GPU değil, FPU kastedilmiş gibi görünüyor. Eskiden FPU asenkron hesaplama yapardı; bugün CPU’ya entegre bir parça olarak görülüyor
Bir MMIX uygulaması yazmayı denedim ve özgün giriş değerinin 2^-1021 değerinden büyük olduğu varsayımını koydum
İlginizi çekerse Wikipedia’da da bu fonksiyon ve geçmişi hakkında iyi bir açıklama var: https://en.wikipedia.org/wiki/Fast_inverse_square_root
Bunlardan birkaçını bir araya getirdim: https://github.com/ncruces/fastmath/blob/main/fast.go
İlgili bir StackOverflow yazısı da var: https://stackoverflow.com/questions/32042673/optimized-low-a...
- Eskiden yaptığım 80’lerin sonu tarzı 3D motoru yeniden yazmak için bu tür tekniklerden oluşan bir derleme başlatmayı düşünüyordum; faydalı oldu
- fastmath paketinin benchmarklarını da görmek isterim
Kusur bulma zamanı geldi. float formülünde yazım hatası var; -1^S değil (-1)^S olmalı. İlki her zaman -1 olur
Ham bit desenini yorumlamanın logaritmanın parçalı doğrusal bir yaklaşımı olduğu açıklaması da doğru değil. Mavi grafikteki veri noktaları arasındaki çizgiler gerçekte yok ve bir bitin yalnızca yarısının 1 olarak ayarlanması mümkün değil. Daha çok logaritmanın ayrık sürümüne yakın; gerçekten var olan veri noktaları, yani kırmızı çizgiyle mavi çizginin kesiştiği noktalar, kelimenin tam anlamıyla ölçeklenmiş ve ötelenmiş logaritmaya eşit. Bunun dışında iyi bir yazı
- Pek anlayamadım. Çok küçük, 6 bitlik bir float düşünürsek: 1 bit işaret, 2 bit üs, 3 bit mantis olduğunda [010000, 010111] aralığı 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75 değerlerini içerir
  Ama bu sayıların 2 tabanındaki logaritmasının ima ettiği mantisler sırasıyla .0000000, .0010101, .0101001, .0111010, .1001010, .1011001, .1100111, .1110100; ilki hariç float’taki 001, 010 vb. ile aynı değil. [2,4) aralığındaki float’lar doğrusal aralıklıdır ama karşılık gelen logaritma öyle değildir; bu yüzden yazının söylediği gibi float, logaritmanın parçalı doğrusal yaklaşımı olarak görülebilir
- Sürekli bir parçalı doğrusal yaklaşım değil, ayrık bir parçalı doğrusal yaklaşım. Mavi çizginin sürekli olmadığı doğru, ama yorum yanlış. Mavi grafik yalnızca birkaç kesişim noktasından değil, x ekseni boyunca eşit aralıklarla yerleştirilmiş 256 ayrı noktadan oluşuyor
  Tüm grafik olsaydı parçalı doğrusal desen içinde 2^32 seçenek olurdu; ama asıl yazının çizdiği şey böyle bir tüm grafik değil. Yazı 32 bit tamsayı ve IEEE-754 32 bit float işlemlerini ele aldığı için açıklamada “ayrık” kelimesinin atlanmasını kabul edilebilir buluyorum
İlginç kavramların çoğunu açıklayan iyi bir yazı, ama bir bölümdeki cebirsel açılım şaşırtıcı derecede kötü
“İlk biçimden bu biçime gitmenin kesin adımları çok, ama bütünlük için hepsini ekledim” sonrasındaki açılımda gereksiz adım çok ve birbirini götüren birkaç işaret hatası var. Özellikle ikinci satırdan üçüncü satıra geçerken eksi işareti düzgün dağıtılmamış. İkinci satırdan sonrasında y_n+1 = y_n + (1 - x * y_n^2) / y_n^2 * (y_n^3 / 2) ile başlayıp y_n+1 = y_n (1.5 * y_n - 0.5 * x * y_n * y_n) sonucuna çok daha kısa biçimde gidilebilir; ara adımlar da doğru olur. Cebir bilen biri için bunların yalnızca bariz adımlar olduğunu düşünüyorum
Ünlü kod parçasındaki magic number en iyi sabit değil. Başka bir sabit kullanırsanız göreli hatayı muhtemelen yaklaşık %0,5 daha azaltabilirsiniz
O dönemde mutlak optimum değeri bulmak zor olmuş olabilir, ama bugün nispeten kolay. Ben de bir ara bu tavşan deliğine düşüp (1/x^2) ve (1/x) için en iyi magic number’ı bulan bir Jupyter notebook’a sahip olmuştum
- Bu soruyu inceleyen makalenin bağlantısı yazının alt tarafında var
Bu yazıda bana en ilginç gelen şey “How Java's Floating-Point Hurts Everyone Everywhere” bağlantısıydı: https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf
Yazarı, “Old Man of Floating-Point” olarak da bilinen William Kahan: https://news.ycombinator.com/item?id=29042853 - An Interview with the Old Man of Floating-Point (1998)
- Konudan bağımsız olarak JAVAhurt PDF’sini okumaya başladım ve dizgisi berbat. Kelime aralıklarını aşırı, üstelik düzensiz biçimde açan bir TeX paketi kullanılmış gibi; başka bir belgeden OCR yapılıp fazladan boşluklar eklenmiş hissi veriyor
  Sabit genişlikli yazı tipi bölümlerinde de tuhaf ek aralıklar var. Okumaya odaklanmak gerçekten zordu; gerçekte öyle olmadığını bilsem de neredeyse bilim delisi manifestosu gibi hissettirdi
Daha önce gördüğüm şu video gerçekten çok iyiydi: https://www.youtube.com/watch?v=p8u_k2LIZyo

Hızlı ters karekök algoritması hakkında bilinmesi gereken her şey

Quake kodu ne yapıyor

IEEE-754 32 bit float gösterimi

Float bitlerine tamsayı olarak bakınca ortaya çıkan logaritmik ilişki

Ham bitler log2(x) için bir yaklaşımdır

Ters karekökü logaritma özdeşliğine çevirmek

0x5f3759df sabitinin gerçek kimliği

Basit bir hack olmamasının nedeni

Newton-Raphson düzeltmesiyle hatayı azaltmak

Bölmesiz Newton düzeltme formülü

Kökeni ve ilgili algoritmalar

İlgili okumalar

2 yorum

Hacker News yorumları

Ham bitler `log2(x)` için bir yaklaşımdır

`0x5f3759df` sabitinin gerçek kimliği