Fast Inverse Square Root
Sayısal İşaret İşleme (DSP)’de sıklıkla bir vektörü normalize etmek için kullanılan bir yöntemdir. 1/kök(x) ( IEEE 754 standardında oluşturulan ikili tek duyarlıklı 32-bit kayan nokta olan x’in karekökü’nün çarpım tersi) için bulunacak sonuca yakınsamak için kullanılan bu algoritma o zamanki bilgisayarlar için hesaplanması oldukça güçtü.
1999 yılında çıkan, Birincil Şahıs Nişancı(FPS) oyunları arasında gelmiş geçmiş en iyi örneklerinden biri olan QUAKE III Arena’nın baş programcılarından John Carmack, idTech 3 (2002 yılında Carmack tarafından GNU General Public License v2.0 olarak yayınlanmıştır) oyun motorunda, 3 boyutlu grafiklerin yansıma, aydınlatma ve gölgelendirme hesaplamalarında sıklıkla kullanılan bu karekök fonksiyonunu gerçekleştirdi.
Yerleşik sqrt(x) fonksiyonu: 404.029ms (%0 hata ile)
Intel Streaming-SIMD ssqrts: 200.395ms (%0 hata ile)
John Carmack’ın tekniği: 72.682ms (0.0990% hata ile)
QUAKE III Arena kaynak kodlarında geçen fonksiyona bakalım:
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // what the f*?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
return y;
}
C dilinde yazılan, float tipinde bir sayı alarak başlayan fonksiyon, ilk satırlarında bizi oldukça basit karşılıyor.
long i; -> 32 bit sayı
float x2, y; -> 32 bit ondalıklı sayı
const float threehalfs = 1.5F; -> 32 bit ondalıklı sayı
Ancak burada bir sorun var, eğer ki 32 bit ondalıklı sayıları binary olarak göstermek istersek
00000000 00000000.00000000 00000000
olarak düşünebiliriz ancak tam ve ondalıklı sayılarla oynadığımızda
00000000 00000010.00000000 00000000 -> 2
00000000 00000100.00000000 00000000 ->4
0000000000000100.00100000 11000100 ->4.128
görüldüğü gibi ~+-2 milyar’a kadar değer tutabilen 32 bit ondalıklı sayıyı tam ortasından ikiye ayırarak kapasitesini 32 bine kadar indirdik. Bu oldukça büyük bir kayıp.
Ancak aramızda yaşayan Carmack gibi dahilerin düşünebileceği bir başka yolumuz daha var.
Matematik kitaplarında sıkça gördüğümüz üstel gösterimi hatırlıyorsunuzdur
4.3.10¹⁴ -> 430 000 000 000 000
aynı gösterimi ikili sisteme uyguladığımızda
0.0001101010 -> 1101010.10-³
IEEE 754 standardında tanımlanan ikili tek duyarlıklı 32-bit kayan nokta gösterimine ulaşıyoruz.
Sign: 0 ise pozitif, 1 ise negatif sayı olduğunu belirtir.
Exponent: Üstel 10^n gösterimindeki n’in değerini gösterir.
Fraction(Logaritmanın ondalık kısmı): 0'dan sonraki ondalıklı sayıların binary değerini gösterir.
Ancak bu şekilde yalnızca pozitif üstel değerlere ulaşabiliyoruz bu yüzden Exponent aralığını belirli bir hata payı ile yarıya indirmeliyiz [-128,127]
x2 = number * 0.5F;
Matematiksel hesaplamalarda Logaritmanın ondalık kısmını onluk tabanda gösterirken 1'den 10'a kadar yazarız(2.10⁴, 2.7⁴ gibi).
İkili gösterimde 1'den 2'ye kadar göstermek isteriz bunu da Logaritmanın ondalık kısmınının ilk hanesini nokta ile ayırarak gösterebiliriz.
Fraction: 0.000000000000000000000 -> 0
Fraction: 1.000000000000000000000 -> 1
Fraction: 1.111111111111111111111 -> 2
Ancak bu gösterim de verimsiz olacaktır.
İkili gösterimde 0'dan farklı tek bir sayı vardır(1). Bu sebeple bu haneye saklamanın bir anlamı olmadığından noktayı bir basamak sola kaydırarak daha verimli hale getirebiliriz.
Üstel sayılarımızı 1'den 2'ye ve 32 bit sayımızın kapasitesini 0'dan 2³²’ye kadar olan tüm sayılardan 0 ile 1 arasına alarak ikili gösterime uygun hale getirdik.
Logaritmanın ondalık kısmı ve üstel sayıyı ikili olarak düşünelim.
Örnek olarak:
M = Logaritmanın ondalık kısmı -> 01001110010000000000000
E = Üstel Sayı -> 11010001
olsun.
Bu sayıyı ikili olarak göstermek logaritmanın ondalık kısmı ile E’yi birleştirmeliyiz. Bunun için E’yi 2²³ ile çarpmalıyız.
M = Logaritmanın ondalık kısmı -> 01001110010000000000000
E = Üstel Sayı -> 11010001
E*2²³ = 11010001 00000000000000000000000*Anlaması ve anlatması oldukça zor olan kısmı bitirmeden, güncel nümerik analiz bilgimin oldukça üstünde olan bu konuyu anlatırken, sık sık yaptığım hatalar ve eksiklikler olabileceğini belirtmek isterim.
“evil floating point bit level hacking”
long ve float tiplerinin IEEE 754 standardındaki en büyük farkları bit manipülasyonu yapılıp yapılamadığıdır.
long tipindeki bir değişkene bit manipülasyonu yapılabilir örnek olarak long değişkeninin ikili değeri 10000000 olsun
10000000 = 128
eğer sol bit kaydırma uygularsak
100000000 = 256
ikiye katlanacak, sağ bit kaydırma uygularsak
1000000 = 64
yarıya inecektir.
Koda geri dönelim…
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
Burada float tipindeki y değişkenini, bit shift yapabileceğimiz long değerine çevrimini görüyoruz. y değişkenindeki ondalık değerleri ve bitleri kaybetmemeli ve birebir şekilde long tipine dönüştürmeliyiz. Bu sebeple değişkenin değerini dönüştürmek yerine bellek adresini dönüştürmeliyiz.
( long * ) &y;
sırasıyla y değişkeninin bellek float tipindeki adresini alıp long tipindeki bellek adresine yönlendiriyoruz. Adres değişmemesine rağmen C derleyicisi bu adreste bulunan sayıyı artık long tipinde görmekte.
* ( long * ) &y;
Ardından adresdeki sayıyı okuyacak pointer’ı ekliyoruz.
“What The F*”
Şimdi amacımızı kısaca hatırlayalım, fonksiyondaki y değişkenine gelecek saya 1/karekök(y) işlemi uygulamak istiyoruz ancak bu standart yolla yapıldığında oldukça maliyetli bir işlem.
IEEE 754 ile y değişkenini float tipinden long tipine dönüştürdüğümüzde oldukça ilginç bir sonuçla karşılaşıyoruz, y değerinin bir şekilde dönüşümden önceki değerinin logaritmasına eşit olduğunu görüyoruz.
Bu da bize 1/karekök(x) ile çalışmak yerine log(1/karekök(x)) ile çalışma kolaylığı sağlıyor.
log( 1/karekök(y) ) = log(y^-1/2) = -1/2log(y) )
Ancak burada görüldüğü üzere hala -1/2 ile bölme işlemi yapmak zorundayız ancak bunu daha önce bahsettiğim bit kaydırma ile çözebiliriz!
0x5f3759df - ( i >> 1 ); -> sağa bit kaydırma ile -1/2 hesaplama
Güzel! emin adımlarla anlayarak ilerlediğimizi varsayıyorum… Peki ya bu garip 0x5f3759df sayısı da neyin nesi?
Oldukça uzun ve ağır matematiksel ifadeler içeren açıklaması olduğundan özet geçmek gerekirse (detaylar için):
Hatırlarsanız, sayıları ikili tabana çevirirken, Exponent aralığını yarıya indirirken belirli bir hata payı ile çevirdik ve Fraction(Logaritmanın ondalık kısmı) hesaplarken de üstel sayıyı ondalık kısımla birleştirmek için üstel ifadeyi 2²³ ile çarpıp bit kaydırma yapmıştık.
x2 = number * 0.5F; -> Exponent aralığının yarıya indirilmesi
M = Logaritmanın ondalık kısmı -> 01001110010000000000000
E = Üstel Sayı -> 11010001
E*2²³ = 11010001 00000000000000000000000
0x5f3759df değeri ise bu dönüşüm ve bit kaydırma işlemlerinden arta kalan hata payının değeridir.
Karekök değerimizi bu bit manipülasyon ve kaydırma teknikleriyle standart teknikten onlarca kez daha verimli ve hızlı bir şekilde hesaplamayı gerçekleştirdik.
y = * ( float * ) &i;
Bu satır ile de yaptığımız çözümün bitlerini geri alabiliriz.
Ancak bu çözüm hata payı sebebiyle yalnızca bir yaklaşık değer. Bu sorunu da Newton–Raphson Metodu yaklaşımı ile çözebiliriz.
Newton–Raphson Metodu’nu kısaca özetlemek gerekirse, Metot verilen fonksiyonun bir kökünü (diğer bir deyişle f(x)=0 değerini sağlayacak bir x değerini) iteratif bir şekilde yaklaşım yaparak bulur.
y = y * ( threehalfs — ( x2 * y * y ) ); // 1st iteration
Carmack’ın Newton-Raphson Metodu’nun birinci iterasyonunu uyguladığını ve yeterli yaklaşımı yaptığını görüyoruz.
(1/y²) -x = 0 => y = 1/karekök(x)
Ancak ikinci iterasyonu yorum satırı olarak bırakıp algoritmanın daha hızlı çalışmasını sağlıyor.
Yazının sonuna geldik :) buraya kadar okuduğunuz için çok teşekkür ediyorum aşağıda faydalandığım kaynakları bulabilirsiniz.
