Benchmark SSE2 logaritmu

4.6.2009 20:52 | Přečteno: 1179× | Jiné | poslední úprava: 4.6.2009 23:08

Tak jak jsem psal v komentářích předminulého zápisu, spatlal jsem logaritmus v SSE2 instrukcích. Trošku jsem otestoval jeho rychlost a tu jsou výsledky + zdroják.

Měření jsem prováděl 5x a vzal průměrnou hodnotu (použil jsem time). Počet iterací byl 500000000. SSE2 logaritmus to zvládnul za 21.52s a logaritmus z C math knihovny 35.06s.

; log_x86_64.asm

BITS 64

SECTION .text
global sse2_log_base_e

; SSE2 implementation of natural logarithm, x86_64 version
; based on that found at http://gruntthepeon.free.fr/ssemath/
; Modified to operate on doubles instead of floats
;
; input xmm9
; output xmm9
;
; eats xmm10-xmm15
; 

sse2_log_base_e:
    movapd xmm9, xmm0 ; Kvůli C testu a gcc, který první double šoupá do xmm0
                      ; tohle tu nebude v kresliči grafů
                      ; kterej to bude šoupat do xmm9
    movapd xmm15, xmm9
    maxpd  xmm9, [rel double_min_norm_pos]
    movdqa xmm12, xmm9
    
    andpd    xmm9, [rel double_inv_mantisa_mask]
    movapd   xmm14, [rel double_0p5]
    xorpd    xmm10, xmm10
    psrlq    xmm12, 52
    psubq    xmm12, [rel double_exponent_bias]
    movapd   xmm11, [rel double_1]
    
    pshufd   xmm12, xmm12, 0xd8
    cvtdq2pd xmm12, xmm12
    
    orpd     xmm9, xmm14
    cmplepd  xmm15, xmm10
    movapd   xmm10, xmm9
    movapd   xmm13, xmm9
    cmpltpd  xmm10, [rel double_cephes_SQRTHF]
    
    subpd  xmm9, xmm11
    andpd  xmm13, xmm10
    addpd  xmm12, xmm11
    andpd  xmm11, xmm10
    movapd xmm10, [rel double_cephes_log_p0]
    addpd  xmm9, xmm13
    subpd  xmm12, xmm11
    mulpd  xmm10, xmm9
    movapd xmm13, xmm9
    movapd xmm11, xmm12
    
    mulpd  xmm12, [rel double_cephes_log_q2]
    mulpd  xmm13, xmm9
    
    addpd  xmm10, [rel double_cephes_log_p1]
    mulpd  xmm11, [rel double_cephes_log_q1]
    mulpd  xmm1, xmm9
    addpd  xmm10, [rel double_cephes_log_p2]
    mulpd  xmm10, xmm9
    addpd  xmm10, [rel double_cephes_log_p3]
    mulpd  xmm10, xmm9
    addpd  xmm10, [rel double_cephes_log_p4]
    mulpd  xmm10, xmm9
    addpd  xmm10, [rel double_cephes_log_p5]
    mulpd  xmm10, xmm9
    addpd  xmm10, [rel double_cephes_log_p6]
    mulpd  xmm10, xmm9
    addpd  xmm10, [rel double_cephes_log_p7]
    mulpd  xmm10, xmm9
    addpd  xmm10, [rel double_cephes_log_p8]
    mulpd  xmm10, xmm9
    
    mulpd xmm10, xmm13
    mulpd xmm13, xmm14
    addpd xmm10, xmm11
    subpd xmm10, xmm13
    addpd xmm9, xmm10
    addpd xmm9, xmm12
    orpd  xmm9, xmm15
    movapd xmm0, xmm9 ; kvůli gcc a C testu
    ret
    
SECTION .data align=16
    
    ; double
    ;   1b - sign
    ;  11b - exponent (mantisa)   -1023
    ;  52b - significand
    ;  64b - total
    
    double_mantisa_mask:     dq 0x7FF0000000000000,0x7FF0000000000000 ; mask for exponent
    double_inv_mantisa_mask: dq 0x800FFFFFFFFFFFFF,0x800FFFFFFFFFFFFF ; inverted mask for exponent
    double_sign_mask:        dq 0x8000000000000000,0x8000000000000000 ; mask for sign
    double_inv_sign_mask:    dq 0x7FFFFFFFFFFFFFFF,0x7FFFFFFFFFFFFFFF ; inverted nasj for sign
    double_exponent_bias:    dq 0x00000000000003FF,0x00000000000003FF ; 1023 (subtract from exponent to get real one)
    double_min_norm_pos:     dq 0x0010000000000000,0x0010000000000000 ; the smallest double
    
    double_0p5:              dq 0x3FE0000000000000,0x3FE0000000000000 ; 0.5
    double_1:                dq 0x3FF0000000000000,0x3FF0000000000000 ; 1.0
    
    double_cephes_SQRTHF:    dq 0x3FE6A09E667F3BCD,0x3FE6A09E667F3BCD ; 0.707106781186547524
    double_cephes_log_p0:    dq 0x3FB204376245245A,0x3FB204376245245A ; 7.0376836292E-2
    double_cephes_log_p1:    dq 0xBFBD7A370B138B4B,0xBFBD7A370B138B4B ; -1.151461031e-1
    double_cephes_log_p2:    dq 0x3FBDE4A34D098E98,0x3FBDE4A34D098E98 ; 1.1676998740000000e-1
    double_cephes_log_p3:    dq 0xBFBFCBA9DB73ED2C,0xBFBFCBA9DB73ED2C ; -1.2420140846000000e-1
    double_cephes_log_p4:    dq 0x3FC23D37D4CD3339,0x3FC23D37D4CD3339 ; 1.4249322787000000e-1
    double_cephes_log_p5:    dq 0xBFC555CA04CB8ABB,0xBFC555CA04CB8ABB ; -1.6668057665000000e-1
    double_cephes_log_p6:    dq 0x3FC999D58F0FBE3E,0x3FC999D58F0FBE3E ; 2.0000714765000000e-1
    double_cephes_log_p7:    dq 0xBFCFFFFF7F002B13,0xBFCFFFFF7F002B13 ; -2.4999993993000000e-1
    double_cephes_log_p8:    dq 0x3FD555553E25CD96,0x3FD555553E25CD96 ; 3.3333331174000000e-1
    double_cephes_log_q1:    dq 0xBF2BD0105C4242EA,0xBF2BD0105C4242EA ; -2.1219444000000000e-4
    double_cephes_log_q2:    dq 0x3FE6300000000000,0x3FE6300000000000 ; 6.9335937500000000e-1

/* test.c */

#include <stdio.h>
#include <math.h>

#define SSE2_LOG

#ifdef SSE2_LOG
double sse2_log_base_e(double b);
#endif

int main()
{
	size_t i;
    double a;
    double b;
    
    
    for (i = 0; i < 500000000; ++i)
    {
    	a = (double)i;
        #ifdef SSE2_LOG
        b = sse2_log_base_e(a);
        #else
        b = log(a);
        #endif
    }
    
    return 0;
}

yasm -f elf64 -o log.o log_x86_64.asm
gcc -o test log.o test.c -lm

Jediné, v čem sse2 logaritmus moc neválí je přesnost výpočtu, C log to zvládne trochu lépe. Ještě to zkusím trochu vytunit

Leoši,

Prosím, prosím, mohl bys opravit BUG, kdy je všechno moje < a > po náhledu v textovém okně nahrazeno špičatejma závorkama a musím to zas přepisovat na < a >?

Update

A ještě jsem zapomněl zmínit, že ten SSE2 logaritmus je schopen počítat 2 najednou (double v horní i dolní 64bit části xmm registru), takže ten čas C logaritmu vynásobte tak ... 1,8x, ať nežeru

BTW pokud někdo přehlídl kreslič funkce na screenu desktopu, tak tu. Je to sice starší verze, ale moc se toho zatim nezměnilo

       

Komentáře

[RS]

Ten bug je nahlasenej, pokud se nepletu. Tak to priste muzes pripsat na bugzillu. BTW: to tak fakt zatim nikomu tolik nevadilo, ze to neni opraveny?

BTW: Dobry reseni je zkopirovat si text pred odeslanim do nejakyho editoru, provest nahled a zkontrolovat a kdyz je to OK, tak text z textovyho pole vymazat a nakopirovat zpatky original z textovyho editoru:)

[Jardík]

No ono by to možná šlo nechat i ve schránce

[Jardík]

já nevim, proč mi dáváte takový hodnocení, sám jsem si dal špatný. To tak milujete assembler? Já vim, vy milujete mě

[Miloslav Ponkrác]

Problém SSE věcí je, že dělají většinou jenom část věcí proti tomu co dělá FPU, či math knihovna.

Klasická math knihovna, či FPU dělá:

1) test definičního oboru a infinite x

2) výpočet hodnoty funkce (s přesností o pár bitů vyšší, než je výstup)

Problém je, že čím větší přesnost, tím je výpočet pomalejší. FPU počítá na 80 bitů plus několik dalších bitů interně, takže přesnost je velmi vysoká. Těmi několika bity navíc se velmi potlačuje vliv zaokrouhlovacích chyb při výpočtu.

Dále C standardní knihovna bohužel kromě výpočtu ještě nastavuje různé math error kódy, takže je to dost brzdí.

Zatím jsem ještě nepotkal žádný SSE algoritmus, který byl použit pro rychlostní srovnání, který by dělal plný rozsah akcí toho, s čím se srovnává.

[Jardík]

Sem ještě nezkoušel, co to udělá na inf/nan a 0, ale při ln(x) x < 0 to vrací NaN

[Jardík]

Při 0 mi to vrátí -inf, C knihovna NaN, správně je -inf. Tak neotravujte

[Jardík]

BTW když s tím kreslím funkci, do okna, který má pár stovek pixelů, tak opravdu nepotřebuju přesnost na 100 desetinných míst (leda při nějakém extra zoomu, kterej bude max. 1px ~~ 1, default mám 40px ~~ 1 a pro každý pixel počítám 4 hodnoty, abych "vylepšil" nespojité funkce)

[Jardík]

Sakra, zase vedu monolog. Ještě, aby tu ten první komentář nebyl a měl jsem nějaký halucinace. ǔz konečně přestaň hulit :*)

[Jardík]

Ono jde vlastně o druhý zoom, takže kecám sračky, ten je zatím neomezený, ale jelikož to přibližuje po pixelu, tak snad nikdo do nekonečna neodscrolluje

[Jardík]

Né fatk, asi bych už měl jít spát. Jsem se trochu překoukl, -inf vrací C knihovna, takže to má dobře ta :-/ každopádně stejně jako do grafu nenakreslim -inf, tak tam nenakreslim NaN, tak mě to vyjde na stejno ... Dobrou noc, Járo

[Miloslav Ponkrác]

Při 0 mi to vrátí -inf, C knihovna NaN, správně je -inf. Tak neotravujte 

Základní problém je, že ani netušíte, co Váš algoritmus dělá, ale musíte to zkusit, namísto toho, abyste věděl, co jste navrhnul.

Takže já bych měl důvěru ve Vaše algoritmy takovou, že bych raději použil cokoli jiného.

Možná na běžných algoritmech se dá hádat, ale v numerické matematice musíte vědět. Jinak vycházejí blbosti a testováním podprogramu jako black box na skutečnou příčinu samozřejmě nepřijdete.

Jinak matematické funkce nad float čísly se neřídí čistou matematickou, protože počítačový float má velmi daleko k ideálním reálným číslům. Pro běžného programátora je rozdíl mezi „pseudoreálnými čísly typu floating point“ a dokonalými reálnými čísly, které používá matematika celkem zanedbatelný. Ale při takových výpočtech je nutné s ním počítat. Například čísla typu floating point nejsou spojitou množinou, není u nich zaručena asociativnost jako u dokonalých reálných čísel, atd..

Skutečné výsledky nad floating point čísly se proto v limitních případech neřídí vždy podle matematiky, ale podle celosvětově a jednotně přijaté normy IEEE-754, která přesně a naprosto do posledního detailu definuje, jak se to má chovat. Tudíž správně je log(0) = NaN.

 

[Grunt]

Můžu se zeptat kde se toto učí a nebo se s tím prostě musí člověk narodit a mít to v genech?

[Jardík]

Ten algoritmus, jak řekl níže __dark__ jsem si nevymyslel, ve drojáku máte odkaz, kde najdete původní implementaci v C a mmx/sse a sse2 intrinsic. Já to pouze přeplácal do assembleru a upravil, aby to počítalo s doublem místo singlu/floatu. Taky někde píšu, že jsem se překouk a vrací to NaN a že -inf vrací ten céčkový log. Stěžujte si tedy na vývojáře, že to nemají podle IEEE-754, když to tak vyžadujete (kdo ví, jestli to tam vůbec je, já ji nezkoumal, akorát jsem na wiki sebral reprezentaci doublu v binární soustavě a spoléhal, že je tam správně).

[Deleted [8409]]

Ooo tady se zase někdo ozval. Doufám, že nám tu nebudete opět tvrdit, že FPU vypočítá ten logaritmus za stejný čas jako push instrukci :)

A pokud dělá standardní C knihovna něco navíc, co zrovna nepotřebuji, tak nevidím na vlastním logaritmu nic špatného. Dále kdyby jste se v oblasti trochu pohyboval, asi by vám byly povědomé konstanty s prefixem cephes, které si Jarda opravdu nevymyslel a tu aproximaci také ne. Takže ty relativní chyby co ten logaritmus produkuje jsou podle mě v normě a podporu pro nějaký infinitiv ve výpočtu, který se vyplatí tímto způsobem paralelizovat asi opravdu potřebovat nikdo nebude (přesněji nikdo kdo sáhne právě po takovém logaritmu).

[Jardík]

FPU vypočítá ten logaritmus za stejný čas jako push instrukci

[Daniel Kvasnička ml.]

Co mas presne za HW?

[Jardík]

Něakýho starýho pomalýho Athlona 64 X2 5200+ EE s 512K L2 cache do socketu AM2.

[Daniel Kvasnička ml.]

A proc nepouzivas SSE3, tenhle procesor ho prece umi, ne?

[Jardík]

To je sice hezký, ale řekněte mi, co ze SSE3 tam využuju .... tam je pár instrukcí a jsou mi k ničemu.

[Daniel Kvasnička ml.]

Jo takhle. Ja zil v predstave, ze SSE3 je dalsi generace, ktera SSE2 nahrazuje.

[trekker.dk]

Leoši,

Prosím, prosím, mohl bys opravit BUG, kdy je všechno moje < a > po náhledu v textovém okně nahrazeno špičatejma závorkama a musím to zas přepisovat na < a >?

Skoro bych si tipnul, že Leoš tě blokuje...

[Jardík]

[Grunt]

Jednoduchá odpověď. Tam už tě ignorovat nemůže.

[Jardík]

BTW tobě se to taky nahradilo, nebo jsi místo &amp;lt; napsal jen [ampersand lt;]?

[trekker.dk]

Prostě jsem to nabral do schránky a vrazil to sem prostředním tlačítkem, takže se to nahradilo samo. A to nepoužívám ten hnusiwig editor.

[Daniel Kvasnička ml.]

No hele, nechci tobe ani klukum z AMD kazit radost... ale muj stary Mac Mini s jednojadrovym PowerPC 7447A (G4) na 1.25 GHz to pri pouziti optimalizace pro AltiVec spocita za 26 az 29 vterin... jo a jeste malickost: v Jave ...s pouzitim JNIAltivec, pricemz cas je meren pres System.currentTimeMillis() na zacatku a na konci main(), tzn. pocita se do toho i inicializace JNI knihovny, alokace poli pro vysledky apod. (pri spousteni z konzole jako JAR a mereni s pomoci time je to asi 31 vterin, protoze se tam pocita cely start VM apod.)

Ani radsi nechtej vedet jak by to dopadlo, kdybych vzal PowerPC procesor, ktery je chronologicky a taktem ekvivalentni tvemu (napr. 970MP na 2.5 GHz, uvedeny zhruba rok pred tvym procesorem)...

[Jardík]

To ví každej, že jsou x86 a x86_64 procesory šitky Ale taky ten můj procík je úsporná verze s 512kB L2 cache, takže je to takový pomalý kšunt. Věřím, že kdybych to testoval na Intel Core 2 Duo na NB, tak to dopadne líp.

[Daniel Kvasnička ml.]

S Core 2 Duo jsem to taky porovnaval, pomoci toho benchmarku dodavaneho k JNIAltivec. Dokonce i tam byl ten Mac Mini u nekterych funkci schopny to Duo (pouzivajici cistou Javu) porazit. Ovsem cista Java asi neni nijak optimalizovana pro SSE, nevim... takze mozna to neni uplne fer.

[Deleted [8409]]

No java je určitě optimalizovaná pro SSE2, i ty funkce v mac-mini budou optimalizované, existuje na to dokonce i nějaký framework.

x86 převálcuje v budoucnu všechny

[Daniel Kvasnička ml.]

Pokud uz java optimalizovana je, pak to s tim valcovanim teda porad nejak nevychazi .... Ted jeste nove ARM Cortexy, ktere taky maji svoji SIMD jednotku a ktere nejspis Atomum pekne zatopi u zadnice... Ja pevne doufam, ze to bude x86, kdo v budoucnu zemre. Ma na to vek a zdravotni stav a nazivu ji drzi jen Microsoft se svym monopolem. Ne nadarmo je soucasny nejrychlejsi pocitac na svete tvoren ze 2/3 PowerPC procesory...

Existuje Accelerate.framework, ale na ten bych si asi musel napsat JNI most, aby fungoval v Jave a to se mi nechce. Navic kdyz uz Java vyuziva SSE...

[Deleted [8409]]

Ty léta výzkumu a vývoji překladačů pro x86, to by bylo škoda vyhodit, ne? Já osobně mám tu architekturu rád, a radši si koupím x86 procesor, než nějaký powerpc nebo arm Mluvím samozřejmě o desktopech (na mobilech fandím armu, je to úsporný procesor).

[Daniel Kvasnička ml.]

Ty léta výzkumu a vývoji překladačů pro x86, to by bylo škoda vyhodit, ne?

Ne To je jak kdybys rikal, ze by bylo skoda zapomenout Michaela Jacksona, kdyz uz ma za sebou tolik plastickych operaci a tak se snazil... x86 je proste neperspektivni a je zbytecne ji priohybat (napr. nebyla pri svem poceti na rozdil od PPC tvorena s planem na 64bitu), kdyz jsou tu jine, cistsi, architektury...

Já osobně mám tu architekturu rád, a radši si koupím x86 procesor, než nějaký powerpc nebo arm

Tak jiste, z pragmatickeho pohledu neni co resit. Apple se v pondeli s PPC rozloucil definitivne, Windows (uz) na tom taky nejedou a pocet linuxovych dister, ktera PPC berou alespon trosku vazne, se taky spise snizuje a kdyz, tak je stejne optimalizace pro tu platformu temer nulova... V embedded zarizenich pro divokou matematiku (medicina, armada, avionika, prumysl, supercomputing, herni konzole) ale PPC prospiva, protoze tam nepotrebuje za sebou tahat onu kouli ve forme "musis podporovat Windows, jinak jsi nezajimavy"...

[Deleted [8409]]

No tak až bude PPC architektura perspektivní, tak o ní můžu uvažovat:) Ale v současnosti neznám jediný důvod, proč bych se takovou architekturou měl zabývat. Nemyslím to tak, že bych tu architekturu nesnášel, ale z pohledu současnosti (a podle mě i budoucnosti) má smysl věnovat se architektuře, kterou má na desktopu každý. Pokud se bavíme o serverovém trhu, tak tam to bylo vždycky jiné, ale třebá i já mám server x86, a popravdě neznám důvod pro výběr jiné architektury. Ono ten graf nejvíce používaných architektur je podle mě celkem jasný (x86 válcuje i servery, a zachvilky bude jiná architektura v globálním měřítku ojedinělá).

[Daniel Kvasnička ml.]

Nemyslím to tak, že bych tu architekturu nesnášel, ale z pohledu současnosti (a podle mě i budoucnosti) má smysl věnovat se architektuře, kterou má na desktopu každý.

To je prece jak kdybys tvrdil, ze je lepsi se vykaslat na Linux a venovat se Windows, protoze to stejne ma na desktopu kazdej... :/

Pokud se bavíme o serverovém trhu, tak tam to bylo vždycky jiné, ale třebá i já mám server x86, a popravdě neznám důvod pro výběr jiné architektury.

Tak ono popravde to PPC ani pro ten serverovy trh neni prilis zajimave, protoze na serveru (ted myslim web, posta, enterprise) tezko nejak casto vyuzijes vypocty v plovouci radove carce... Kdyz mas ale aplikaci, ktera ty vypocty nekde potrebuje, tak je dobre mit smiseny cluster a vypocty prenechavat nodum s PPC procesorem, na kterych bezi specialne optimalizovany SW (tady se dostavame zas k jine veci: je dost mozne, ze cena PPC reseni + optimalizace SW nakonec bude vyssi, nez kdyz se to same vyresi hrubou tupou silou x86 serveru a neoptimalizovaneho SW.... coz je debilni, ale nejspis to tak bude :/). A to se dotyka i desktopu. Hudba, video, desktopove efekty. Apple nelhal, kdyz rekl, ze Power Mac je prvni pocitac, ktery umoznuje lidem mit doma na stole 1 GFLOPS. Muj staricky Power Mac na 733 MHz ma 1.3 GFLOPS, pokud se pouzije AltiVec. Intel tohohle byl schopen pouze s nekolikanasobkem taktu procesoru...

Ovsem je fakt, ze v dobe OpenCL se CPU celkove posouva do trochu jine role... zustava mu ta "nudna" prace a k zajimavym vecem se dostane cim dal mene. Otazka je, jestli to prave nebylo iniciovano tim, ze x86 -- jakozto vsudyprtomna a de fakto diky MS nenahraditelna architektura -- konstantne nezvladala v FP vypoctech nejak vic zazarit...

[Deleted [8409]]

Nepřirovnával bych architekturu PC k Windows. Architektura x86 se asi prosadila díky Windows (nebo na to měl Microsoft vliv), ale musíte uznat, že se jedná o levné procesory, které jsou poměrně kvalitní (alespoň podle mě jo). A tady je zakopaný pes, proč by měl někdo podporovat dražší architekturu, když existuje mnohem levnější a z časového hlediska naprosto vyzkoušená architektura, pro kterou jsou psané skoro všechny programy?

Intel zkusil vymyslet novou architekturu, a jak to dopadlo ?

Jinak s tím posledním odstavcem souhlasím, ale na druhou stranu, jaký by v tom případě mělo smysl vyvíjet procesory? Výrobci se vždycky budou snažit o vývoj a hlavně prodej, takže se karta může i obrátit, a budou výkonné procesory s možností provádět grafiku tam (myslím, že Intel přímo něco takového dělá, ale neznám podrobnosti).

[Daniel Kvasnička ml.]

A tady je zakopaný pes, proč by měl někdo podporovat dražší architekturu, když existuje mnohem levnější a z časového hlediska naprosto vyzkoušená architektura, pro kterou jsou psané skoro všechny programy?

Ale jak vis, ze PPC architektura je drazsi? Pocitace Apple v dobe PPC staly hodne, to jo. Ale to prece neznamena, ze ta architektura je draha. Pokud bys objednaval ve velkych mnozstvich, tak si myslim, ze bys klidne mohl postavit PPC desktop se stejnym pomerem vykon/cena jako x86, ne-li lepsim. Otazka je, co bys na tom provozoval OS X tezko, Win ne, Linux s nulovou optimalizaci...?

V soucasny dobe je to fakt dobry tak akorat na vypocty a nebo casem, az se Linux vic rozsiri a bude stat za to do toho investovat, tak pro desktopy a hlavne pro n(et|ote)booky, jelikoz vykon/watt byl u PPC (kdyz nepocitam serverove POWER a z nich odvozene G5) vzdycky dobry...

[trekker.dk]

protoze tam nepotrebuje za sebou tahat onu kouli ve forme "musis podporovat Windows, jinak jsi nezajimavy"...

"Koule" podpory Windows je spíš než co jiného jenom výmluva. Takový Efficeon Windows podporoval, výkonově byl na tom podobně jako odpovídající procesory x86, spotřeba podstatně menší, ale ani tak se neuchytil.