Kirjautuminen

Haku

Tehtävät

Keskustelu: Koodit: C: Neliöjuuriapproksimaatio IEEE-liukuluvuilla

Sivun loppuun

Metabolix [22.06.2006 03:23:06]

#

Neliöjuuri on perin hidas laskutoimitus, ainakin vanhemmilla prosessoreilla, mutta peliohjelmoinnin ja erityisesti 3D-matematiikan puolella se on hyvin merkittävässä osassa. Tässä on kaksi funktiota, jotka ovat laitteisto- ja pahimmillaan myös alustariippuvaisia, mutta laskevat neliöjuuren suunnilleen oikein huomattavasti alkuperäistä tapaa nopeammin.

Perustana kummallekin funktiolle on liukuluvun esitysmuoto muistissa. Tyypillisesti luvun alimmissa biteissä sijaitsee itse luku siten, että se on yksi tai suurempi mutta alle kaksi. Ylempänä on eksponentti ja ylimmässä bitissä vielä etumerkki. Luvut siis esitetään tällaisessa muodossa:
etumerkki * luku * 2eksponentti

Jos "luku" on yksi, neliöjuuri onnistuu tietenkin jakamalla eksponentti kahdella.
1*26 = 64
1*26/2 = 1*23 = 8
Vaan entäpä, jos eksponentti onkin pariton? Silloin sen puolitus pyöristyy alempaan, mutta samalla sen alin bitti siirtyykin osaksi lukua, jolloin luku kasvaa puolikkaalla.
1*21 = 2
1,5*20 = 1,5
Tässä kohti tuleekin maksimaalinen virhe. 2:n neliöjuuri on suunnilleen 1,414, eli virhe on noin 6%. Näiden ääritapausten välillä virhettä tulee aina jokunen prosentti. Tämä tarkkuus riittää moneen tarkoitukseen jo vallan mainiosti.

Parannuksena ja hidastuksena tähän tulosta voi tarkentaa pienellä ylimääräisellä laskemalla. Kahden neliöjuureksi tuli 1,5. Siispä 2 / 1,5 on selvästikin 1,5 - vaan eipä olekaan. Siitä tulee pyöreästi 1,333. Neliöjuuren oikea arvo löytyy siis näiden väliltä, joten vastaukseksi kelvatkoon niiden keskiarvo, tässä tapauksessa noin 1,4167. Yhä lähemmäksi oikeaa arvoa pääsisi jatkamalla samalla tavalla vielä toisenkin kerran. Ensimmäisenkin jälkeen virhe on kuitenkin jo niin paljon pienempi (0,17%), että jatkaminen ei useimmiten kannata.

Itse algoritmin tarkempaan toimintaan tutustuakseen voi tutustua ensin IEEE:n liukulukustandardiin. Näin se toimii teknisistä syistä. Lukuja käsitellään kokonaislukuina long ja long long, koska ne sattuvat olemaan kooltaan yhtä isoja kuin liukuluvut mutta vain kokonaisluvuille voi tehdä tarvittavia bittioperaatioita. Etenkin vanhoissa koneissa ero voi olla merkittävä, nykyään liukulukuoperaatiotkin ovat jo melko nopeita.

Muoks. Kommenttikeskustelussa selvisi monta asiaa ja testejä tuli tehtyä runsaasti, joten muokkasin vinkkiä hieman tehokkaampaan suuntaan. Suosittelen kaikille kommenttienkin selaamista. Inline-Assemblyllä tehty versio pesee kaikki muut miten tahansa päin, joten sitä suosittelen, jos vain suinkin mahdollista.

// Maksimaalinen virhe:
// 1,5 / sqrt(2) - 1
// 6,1%
#define pika_sqrtf(vastaus, luku) \
    *(unsigned long *)&(vastaus) = ((*(unsigned long *)&(luku)) >> 1) + (127UL << 22)

#define pika_sqrt(vastaus, luku) \
    *(unsigned long long *)&(vastaus) = ((*(unsigned long long *)&(luku)) >> 1) + (1023ULL << 51)
// Muistathan, että funktiot ovat hitaampia...
inline float f_pika_sqrtf(float val)
{
	pika_sqrtf(val, val);
	return val;
}

inline double f_pika_sqrt(double val)
{
	pika_sqrt(val, val);
	return val;
}
// Maksimaalinen virhe:
// ((1,5 / sqrt(2)) + (sqrt(2) / 1,5)) / 2 - 1
// 17 * sqrt(2) / 24 - 1
// 0,17%
inline float tarkempi_sqrtf(float luku)
{
	float juuri;
	pika_sqrtf(juuri, luku);
	return ((luku / juuri) + juuri) / 2.0f;
}
inline double tarkempi_sqrt(double luku)
{
	double juuri;
	pika_sqrt(juuri, luku);
	return ((luku / juuri) + juuri) / 2.0;
}

Deewiant [09.07.2006 20:18:18]

#

Omalla koneellani todellakin standardi sqrt oli noin kaksi kertaa nopeampi, sekä tietysti myös tarkempi, kuin kumpikin noista. 3DNow!:sta, SSE:stä, ja SSE2:sta kun tosiaan löytyy näitä neliöjuuri-instruktioita, joita prossu voi suorittaa ihan suoraan.

Merri [10.07.2006 15:59:58]

#

Koodia voisi optimoida ainakin sillä, että laittaa (127 << 22) sijaan sen lopputuloksen, eli 532676608 tai heksana 1FC00000.

Metabolix [10.07.2006 16:01:26]

#

Jaa, luuletko, ettei kääntäjä optimoi sitä noin automaattisesti? Varmaankin tuosta tulisi aivan yhtä nopea, vaikka siitä tekisi C++-tyyliin template-funktion, jossa iffillä tarkistetaan, onko parametri float vai double.

Edit: No niin, testattu on molemmat asiat, g++:n assembly-tuloste on täsmälleen sama ja eroa syntyy templateillakin vain funktion nimeen.

Merri [11.07.2006 06:36:43]

#

Mutta minulla onkin käytössä vuonna 1996 valmistunut kääntäjä, jolloin tämmöisillä asioilla on oikeasti väliä, etenkin kun kielestä ei edes löydy bitshiftiä. Eli toisin sanoen kokeilin vastaavaa VB6:lla. Funktiokutsuna sain oman kikkareeni toimimaan nopeammin kuin VB:n natiivi Sqr, yli puolitoistakertaisesti (noin 1.7). Oli laiska olo, joten mittasin vain epätarkan GetTickCountin avulla kauanko kestää tehdä tietty määrä laskutoimituksia.

Tulokset muuttuivat huimasti muuttamalla funktiokutsun inlinekoodiksi: inlinenä koodi oli neljä kertaa nopeampi kuin VB:n natiivi funktio. Joten kysymys kuuluu: miten paljon eroa nykyisillä C-kääntäjillä tulee sillä, että toteuttaakin laskun inlinenä eikä erillisenä funktiokutsuna?


Tässä vielä iki-ihana VB-koodi kansan nähtäväksi:

Option Explicit

Private Declare Sub CopyMemory Lib "kernel32" Alias "RtlMoveMemory" (ByRef lpvDest As Any, ByRef lpvSrc As Any, ByVal cbLen As Long)
Private Declare Function VarPtrArray Lib "msvbvm60.dll" Alias "VarPtr" (Var() As Any) As Long

Private SA1D(5) As Long
Private lngSQR() As Long
Private lngSQRPtr As Long

Public Sub CleanSQR()
    CopyMemory ByVal VarPtrArray(lngSQR), 0&, 4&
End Sub
Public Sub InitSQR()
    SA1D(0) = 1
    SA1D(1) = 4
    SA1D(4) = 1
    CopyMemory ByVal VarPtrArray(lngSQR), VarPtr(SA1D(0)), 4&
End Sub
Public Function FastSqr(ByVal Value As Single) As Single
    If lngSQRPtr <> 0 Then
        SA1D(3) = lngSQRPtr
        lngSQR(0) = (((lngSQR(0) And &H7FFFFFFF) \ 2) + &H1FC00000) Or (lngSQR(0) And &H80000000)
        FastSqr = Value
    Else
        lngSQRPtr = VarPtr(Value)
        SA1D(3) = lngSQRPtr
        lngSQR(0) = (((lngSQR(0) And &H7FFFFFFF) \ 2) + &H1FC00000) Or (lngSQR(0) And &H80000000)
        FastSqr = Value
    End If
End Function

Eli ohjelman alussa täytyy kutsua InitSQR ja lopussa CleanSQR.

Metabolix [11.07.2006 10:06:15]

#

Juu, se on sitten asia erikseen. Turbo C -kääntäjällä (MS-DOSin aikaa) käänsin tuon (ilman inline-määrettä, kun se ei moista vielä tue), ja sekin osasi itsekseen laskea noin yksinkertaisen toimituksen valmiiksi. Väkisin vain tahtoo tulla aina ylimääräisiä nollanlisäyksiä, 32-bittisillä kääntäjillä double-versioon ja 16-bittisellä jo float-versioonkin. Lisäksi jos funktio tietää käsittelevänsä liukulukuja, assembly-tulosteessa näkyy väkisinkin yksittäinen FLD-komento, jota ei varsinaisesti kaivata. Tämän sai pois tekemällä funktiosta version, joka ottaa parametrinaan osoittimen lukuun ja joka siis muuttaa itse lukua. Jos tuosta haluaa inline-funktioita tukemattomalla C:llä inline-version, muoks: lisäsin sen vinkkiin. Vinkin oli tarkoitus selittää, kuinka tuo tehdään. Jokainen voi siitä itse optimoida omaan tyyliinsä sopivan ratkaisun.

Niin, se nolla lisätään, kun pitkää lukua käsitellään kahdessa rekisterissä ja lisättävän toinen puoli sattuu olemaan nolla (0x1FC0, 0x0000):

	add	ax,0
	adc	dx,8128

Ja Merrille vielä: ylintä bittiä ei tarvitse ottaa mukaan. Bittisiirrossa se muuttuu nollaksi C-versiossa, ja muutenkin, jos se olisi ykkönen, kyseessä olisi negatiivinen luku, jonka neliöjuuresta tulee jotakin erikoista, floatin (Singlen) tapauksessa näin TC:llä käännetystä versiosta:

Juurrettava           tulos
x <= -4.0             -0.0
-4.0 < x <= -2.0      +NaN
-2.0 < x < -1.0       Floating point error: Domain.
x = -1.0              +INF
-1.0 < x <= -0.0      Isoja lukuja

Or (lngSQR(0) And &H80000000) on siis nähdäkseni turha.

Merri [11.07.2006 23:40:30]

#

Jep, tarkistelin tuon myöhemmin itsekin ja korjailin pois. Lopputuloksena kaksinkertainen nopeus Sqr:ään verrattuna, jos ei laita virheentarkistusta negatiivisen luvun varalle. Jos tarkistuksen laittaa, niin sitten lennetään taas kovin kovin hitaaksi.

Metabolix [12.07.2006 14:26:41]

#

Tein nyt tuosta testiversion Delphillä. Double-version tulokset näyttävät tältä:

2434 ms Sqrt
2760 ms Juuri
2351 ms AsmJuuri
 253 ms Juuri inlinenä
  86 ms AsmJuuri inlinenä

Nämäkin ovat siis hyvin luotettavia timeGetTime-aikoja, mutta kuten tästä selvästi näkee, inlinenä suoritetut pätkät ovat ylivoimaisia. Kaikista ajoista on vähennetty tyhjän for-silmukan aika, joten ajat ovat lähempänä itse operaatioiden suhteellisia aikoja. Assemblyllä toteutettua versiota vielä tehostin hieman sillä, että bittisiirto suoritetaan vain luvun merkittävämmälle puoliskolle. Lopputuloksena funktio näyttää tältä:

function AsmJuuri(Val: Double): Double; overload;
begin
  asm
    mov eax, dword ptr Val
    mov dword ptr Result, eax
    mov eax, dword ptr Val + 4
    shr eax, 1
    add ax, $1FF8
    mov dword ptr Result + 4, eax
  end;
end;

{Singlelle vastaava, vaikken tätä käyttänytkään...}
function AsmJuuri(Val: Single): Single; overload;
begin
  asm
    mov eax, Val
    shr eax, 1
    add ax, $1FC0
    mov Result, eax
  end;
end;

Vastaavaa assembly-pätkää voi ja kannattaa käyttää siis juuri siinä, missä neliöjuurta tarvitsee. Funktiossa hommaa hidastaa se, että luku pitää palauttaa liukulukurekisterissä, minne se ladataan ylimääräisillä komennoilla. Inlinenä tällaista hidastetta ei ole, kun luku yksinkertaisesti luetaan muistista ja kirjoitetaan muistiin.

Merri [13.07.2006 14:36:08]

#

Vielä löytyi yksi vaihde VB:llekin ja sillä lopputulos on melkein kolme kertaa nopeampi kuin natiivi Sqr. Menetelmään vaadittava koodi vain on jo sen verran pitkä, että se pitäisi lähetyttää omana koodivinkkinään. Lisäksi menetelmä ei vielä toimi IDE:n alla. Itse tehtävän tekevä koodi on kuitenkin varsin minimaalinen:

Public Sub SngSqrREAL(ByVal lngValue As Long, ByRef lngResult As Long)
    If lngValue > 0 Then lngResult = (lngValue \ 2) + &H1FC00000
End Sub

Menetelmä on siis sellainen, että korvataan proseduuri muistissa toisella :)

Metabolix [14.07.2006 11:37:01]

#

Tuo kyllä nopeutuu vielä, kun ottaa negatiivisuustarkistuksen pois. (Ja hassua, että tarkistat, onko vastaava Long-muuttuja negatiivinen, sen paikalla muistissahan on kuitenkin liukuluku. Sattumalta tuo tarkistus nyt toimii kuitenkin, mutta...) Eihän sillä kukaan mitään tee. Mahtaisiko mennä täydestä, jos tuosta tekisi DLL-funktion ja määrittelisikin sen sitten näin:

Public Declare Sub SngSqrREAL Lib "OmaLib" (ByVal Luku As Single, ByRef Tulos As Long)

Merri [15.07.2006 11:04:14]

#

DLL-funktion kutsuminen VB:stä on hitaampaa kuin luulisi. Esimerkiksi CopyMemoryn käyttö ei tässä tilanteessa kannata, koska yksittäisen sellaisen kutsuminen on jo puolta hitaampaa.

Miksi nollatarkistus ei toimisi? Mitä sattumaa asiassa näet? Jos Single on nolla, kaikki sen bitit ovat epäaktiivisia, joka on myös Longissa nolla. Vastaavasti jos ylin bitti on aktiivinen, silloin luku on negatiivinen niin Singlessä kuin Longissa. En siis näe mitään vikaa tai sattumaa nollatarkistuksessa vaan puhdasta logiikkaa.

Nopeustestissä nollatarkistuksen lisäys tai poisto ei vaikuttanut nopeuteen oikeastaan yhtään ja kaiken lisäksi sen lisääminen korjasi toimivuutta. Kun funktiolle syötti nollan kun siinä ei ollut tarkistusta, se antoi jonkin aivan oudon vastauksen; sama tietenkin negatiivisen luvun kanssa.

Metabolix [15.07.2006 14:34:23]

#

No juuri sitä sattumaa, että formaatit sattuvat osumaan kohdakkain tuolla tavalla. Tuo on vähän sama kuin sotkisi C++:ssa eri luokkien osoittimia keskenään. Kyllähän se toimii, jos kaikissa sattuu olemaan suunnilleen samat jäsenet, mutta niin ei silti ole kaunista tehdä.

Tuo tarkistus lisää kaksi assembly-käskyä, mikä on jo aika paljon, kun muuten toimitus on inlinenä vain puolisen tusinaa käskyä kooltaan (neljä, jos sen kirjoittaa itse). Varsinkin hyppykäsky on hidas, ja jopakäsin kirjoitetussa versiossa nopeus puolittuu (se ei nyt päässyt tähän asti). Nämä tulokset tulivat, kun juurrettava oli aina positiivinen:

538 Inline-sqrt
1380 Inline-sqrt tarkistuksella
3560 Funktio-sqrt
3597 Funktio-sqrt tarkistuksella

Tosiaan, funktiokutsusta aiheutuva hitaus ja säätö kumoaa tuonkin merkityksen, mutta sen sijaan inlinenä hidastuminen on huomattava. Ohjelmahan pitäisi tehdä niin, että neliöjuurelle ei voi tulla negatiivisia lukuja. Kyllä, negatiivinen etumerkki rikkoo kaiken, ja nollan neliöjuureksi tulee jotakin hieman epämääräistä. Kuitenkin jälkimmäisessä tapauksessa luku on niin pieni, ettei sillä ole merkitystä, ja ensimmäiseen tapaukseen ei pitäisi koskaan päätyä. Tyypillinen käyttötarkoitus neliöjuurelle on vektorien normalisointi, jolloin juuren alle tulee komponenttien neliöiden summa, joka ei voi olla negatiivinen. Ja onko jollakin tavalla parempi jättää tulokseksi se, mitä tulos-muuttujassa ennestää sattuu olemaan, kuin vaihtaa siihen se epämääräinen tulos, joka laskusta tulee? Lopputulos on kuitenkin väärä. Jos negatiivisiin lukuihin haluaa varautua, voi kehittää aivan erikseen virheentarkastuksen, laittaa vaikka funktiolle virhepalautusarvoksi -1, mutta tässä vaiheessa hyöty on jo menetetty, yhtä hyvin voisi käyttää oikeaa neliöjuurta.

Tässä vaiheessa lienee syytä mainita, että tässä tapauksessa kerrankin itse assemblyllä kirjoitettu versio on todistettavasti nopeampi kuin kääntäjän kääntämä (tai vähintään yhtä nopea, vaikka kuinka hyvän kääntäjän joku tuohon etsisi), eli kenenkään on turha tulla tänne kertomaan, kuinka asm-optimointi on nykyään tyhmää, "kun kääntäjä osaa optimoida niin hyvin".

Merri [15.07.2006 17:04:26]

#

En ymmärrä tätä kauneuskäsitettäsi. Koko homman ideahan on tehdä juuri tätä mielestäsi epäkaunista asiaa, käsitellä "väärää" muuttujatyyppiä sellaisilla tavoilla, jotka eivät normaalisti ole mahdollisia: bitshiftiä ja kokonaislukutason yhteenlaskua liukuluvulle. Sitten kun teen täysin toimivan asian, joka toimii bittitasolla täysin oikein ja virheettömästi, niin se on sitten yllättäen väärin? En tässä tosiaan ymmärrä kauneusihanteitasi tai käsitystäsi "sattumasta". Ymmärrän toki tarkistuksen vähäisen hyödyn, mutta en ymmärrä miten siitä voi vääntää noin ison numeron tai väittää että sen tekeminen tai menetelmä jolla se tehdään on jollakin tavalla väärin.


Asensin sitten ThunderVB:n. Sillä ja MASMilla sai aikaan tämmöistä:

Public Function SngSqr(ByVal Value As Single) As Single
    '#asm'  mov eax,[esp+4]
    '#asm'  shr eax,1
    '#asm'  add eax,532676608
    '#asm'  mov [esp+4],eax
    SngSqr = Value
End Function
Public Function DblSqr(ByVal Value As Double) As Double
    '#asm'  mov eax,[esp+8]
    '#asm'  shr eax,1
    '#asm'  add eax,536346624
    '#asm'  mov [esp+8],eax
    DblSqr = Value
End Function

Nopeutta en ole vielä vertaillut. On vielä puutteita ASMin ja VB:n yhteiselon ymmärtämisessä, en millään keksinyt tapaa saada funktion palautusarvo suoraan ASMilla, joten kiersin ongelman sitten toisella tavalla. Ensimmäisiä kertoja ylipäätään kirjoittelen ASMia.

Merri [15.07.2006 19:08:57]

#

Tämä pyyhkii sitten lattiat:

Public Sub AsmSqr(ByVal Value As Double, ByRef Answer As Double)
    '#asm'  mov ebx,[esp+12]
    '#asm'  mov eax,[esp+4]
    '#asm'  mov [ebx],eax
    '#asm'  mov eax,[esp+8]
    '#asm'  shr eax,1
    '#asm'  add eax,536346624
    '#asm'  mov [ebx+4],eax
End Sub

Seitsemän kertaa nopeampi kuin Sqr.

Metabolix [16.07.2006 11:18:05]

#

Long-tyyppinen palautusarvo palautetaan yleensä eax-rekisterissä, liukuluvut liukulukurekisterissä. Liukulukurekisteriä on kuitenkin hyvä välttää, koska se taas kerran hidastaa asioita. Tarvittaessa arvon sijoitus kuitenkin onnistuu fld-komennolla:
fld qword ptr [doublen_sijainti]
fld dword ptr [singlen_sijainti]
Eli tähän mennessä nähdyistä voi kaikista oikeastaan muuttaa palautustavaksi tuon ByRef-systeemin tai vastaavan, niin nopeus kasvaa taas kerran. Huonona puolena on, että juurta ei voi enää käyttää lausekkeen keskellä, mutta ehkäpä se ei ole liian suuri menetys.

Juu, enpä ajatellut tuolta kannalta. :) Mutta siitä huolimatta tarkistus on hidastus.

Tässä on kerrankin kunnon vinkin ainesta, kun tuli näin pitkä keskustelu ja paljon testailua ja säätöä. ^^ Jos joku jaksaa tämän kaiken lukea läpi, niin tuleepa ainakin asia ymmärretyksi. Merri voisikin heittää nuo omansa VB-puolelle vinkiksi.

Merri [16.07.2006 13:45:04]

#

Eipä tässä kukaan muu näyttänyt olevan tästä kiinnostunut. Ehkä "oheistuotoksista" sun muista voisi vääntää jotakin kasaan. Sqr:ää tarvitsee niin kovin harvoin, että sillä ei ole niin väliä kuin noilla muilla kokeilluilla asioilla.

Metabolix [16.07.2006 22:06:19]

#

Noh noh, äläs nyt. Peleissä neliöjuurta todella tarvitsee, ja eikös noita "3-D sota peli projekteja" tehdä VB-alueella harva se päivä? ;) Jos jokin niistä vaikka pääsisi joskus neliöjuureen asti. :P

Nyt pääsin vihdoin oman koneen ääreen ja testasin tätä kunnolla käytännössä. Muutin yhden neliöjuurikohdan peliprojektini valaistusfunktiossa käyttämään tätä versiota, lopputulos on varsin huima:

* Keskiarvo: 205.20 fps  * (Ennen)
* Keskiarvo: 237.88 fps  * (Jälkeen)

Ei mennyt hukkaan. Tähän tarkoitukseen tarkkuus ainakin riitti. Jos joskus oikein jaksaisin yrittää, voisin kirjoittaa assemblyllä koko valaistusfunktion ja katsoa, mitä siitä tulisi, mutta toistaiseksi se taitaa olla minulle hieman liikaa...


Sivun alkuun

Vastaus

Aihe on jo aika vanha, joten et voi enää vastata siihen.

Tietoa sivustosta