[Alg] vermijden/optimaliseren van 64bit instructies

[Zoijar]

Ik zit met de voglende regels code:

C++:

return static_cast<uint32>(((static_cast<uint64>(a)*key + b) % p_) % r_);

Waar 'p' typisch zo'n aantal honderduizend groot is (meer dan 16bit) en a en b tussen 0 en p liggen, dus neem even aan ook groter dan 2^16. Key kan van alles zijn, waarschijnlijk ook groter dan 2^16. r_ is relatief klein.

Het punt is nu dus dat er altijd een 32bit antwoord uitkomt. Maar voor de vermenigvuldiging heb ik 64bits nodig, anders komt er een overflow en is het antwoord fout. Deze code werkt, met de cast naar uint64, maar mijn profiler geeft aan dat er += 40% van de tijd in _aullrem() en consorten wordt doorgebracht op mijn 32bit amd athlon xp.

Is er een manier om dit te omzeilen? Door de modulo anders te schrijven?

[madwizard]

Het probleem komt waarschijnlijk omdat de div operator zoals de intel x86 die heeft om (unsigned) integers te delen, als deeltal een 64-bit getal neemt (edx:eax), als deler een 32-bit getal en als resultaat een 32-bit quotient (eax) en een 32-bit rest (edx). Dit geeft een probleem als het resultaat niet in 32-bit past, je krijgt dan een overflow.

Omdat een standaard deling met 64-bits getallen zoals in jouw code ook gedaan wordt ook moet werken met quotienten groter dan 32-bit is er een library functie nodig om dit te realiseren, vandaar de aullrem call. aullrem is een functie die van twee 64-bits getallen de 64-bit deler bepaald.

Voor sommige gevallen, zoals een deler van < 32-bits (jouw geval bijvoorbeeld) zijn snelle oplossingen mogelijk, alleen moet aullrem dit eerst bepalen omdat die functie zo'n garantie natuurlijk niet heeft maar gewoon altijd z'n werk moet doen. Je kunt de snelle oplossing ook zelf maken, het scheelt je de check die aullrem doet om te kijken of dit mogelijk is (omdat je het zelf garandeert), en de functie call natuurlijk.

Het hele probleem komt dus neer op dat de div-opcode alleen een 64-bits getal kan delen waarvan het quotient in 32-bits past. Als de deler ook in 32-bits past is dit met een truc mogelijk met 2 div's.

Een 64-bits getal kun je opsplitsen in een hoog gedeelte (H) van 32-bits en een laag gedeelte (L) van evenveel bits. Je wilt de rest van dit getal hebben, dus modulo P:

rest = (H * 2^32 + L) mod P

Delen door P en de rest uitlezen kan soms niet omdat de waarde te groot is. Je kunt dit echter oplossen door H te vervangen door (H mod P), omdat dat modulo P gewoon hetzelfde getal is. Dus je neemt H2 = H % P en krijgt:

rest = (H2 * 2^32 + L) mod P

H2 is gegarandeerd kleiner dan P (het is immers de rest van een deling met P), wat betekent dat als je H2 * 2^32 deelt door P je ook nooit over de 2^32 grens heen kunt komen. Zelfs als L z'n maximale waarde heeft en erbij opgeteld wordt kun je na deling met P nooit een quotient groter dan 32-bit krijgen. De rest zal echter nog steeds hetzelfde zijn als bij de eerste formule, omdat H2 en H gelijk zijn binnen modulo P.

Tot zo ver het wiskundige gedoe , misschien niet nodig om de implementatie te maken maar wel aardig om te weten lijkt me. De implementatie kun je in assembly op de volgende manier maken: neem het hoge gedeelte van het deeltal (H) als 32-bit getal, en bepaal de rest met P. Gebruik deze rest als nieuwe H. Nu kun je het 64-bits getal H:L veilig delen door P zonder overflow te krijgen.

GAS:

; edx:eax is het 64-bits getal waarvan
; de rest bepaald moet worden. [p] is de
; 32-bits deler.
mov ebx, [p]    ; ebx = deler
mov ecx, eax    ; ecx = low 32-bit van deeltal (L)
mov eax, edx    ; eax = high 32-bit van deeltal (H)
xor edx, edx    ; edx = 0 (edx:eax = eax)
div ebx         ; H / P, eax = quotient, edx = deler (H2)
mov eax, ecx    ; edx:eax = H2:L
div ebx         ; H2:L / P, eax = quotient, edx = deler

; rest zit nu in edx

Overigens gebruikt aullrem soortgelijke code, maar dus wel met een extra check of de deler 32-bits is.

[Verwijderd]

C++:

 return  (  (  (a%p_) * (key%p_)  + b) % p_ )  % r_;

zou volgens mij moeten werken

[madwizard]

Verwijderd schreef op 11 juli 2004 @ 22:54:

C++:

 return  (  (  (a%p_) * (key%p_)  + b) % p_ )  % r_;

zou volgens mij moeten werken

Als p_ een getal van enkele 100.000en is en a en key liggen daarbij in de buurt krijg je dan nog steeds een overflow bij de vermenigvuldiging (100.000^2 is al groter dan 2^32).

[Zoijar]

Ja die extra mod p had ik ook al aan gedacht, en dat werkt idd niet omdat p te groot is.

Verder moet ik ook nog deze berekenen trouwens:

C++:

return static_cast<uint32>(((static_cast<uint64>(a)*key + b) % p_) / r_);

Dus met een deling aan het einde. Die twee functies samenvoegen zal al een hoop tijd schelen...

Dit is trouwens mijn profile info:



Waarom lijken die aullxxx() calls zo ontzettend traag? Ok ze worden vaak uitgevoerd, maar samen zo'n 40% van de totale executie tijd?

Bedankt voor die code madwizard, maar ik gebruik liever geen asm om het portable te houden. Misschien dat het sneller is in C met wat AND 0xFFFF0000s etc.? Het zelfde principe dus maar dan in C, zal het eens proberen.

Ik zit ook nog steeds te denken of ik het niet kan herschrijven dat ik al eerder die mod/div door 'r' kan doen. Ik weet wel dat 'p' een priemgetal is en er is dus een modulo inverse die ik kan precomputen met een gcd algoritme. Zou ik daar wat aan hebben hmm...

[EfBe]

En als je de windows only (alleen om te testen) __int64 VC++ type gebruikt, of win32's ULONGLONG ? Loopt het dan rapper?

[Zoijar]

EfBe schreef op 12 juli 2004 @ 13:07:
En als je de windows only (alleen om te testen) __int64 VC++ type gebruikt, of win32's ULONGLONG ? Loopt het dan rapper?

Ohja, uint64 is getypedefed (mooi woord ) als __int64 onder windows/VC (en anders long long). uint32 is gewoon unsigned int.

Ok Dit zorgt al voor 40% performance increase (ten opzichte van dezelfde functie die beide in een keer berekent, anders is het zo'n 60%)

C++:

inline void LossyDictionary::hashquotient(uint32 key, uint32 a, uint32 b, uint32& h, uint32& q) const {
    uint64 mul = static_cast<uint64>(a)*key + b;
    uint32 m = ((static_cast<uint32>(mul >> 32) % p_) + static_cast<uint32>(mul & 0xFFFFFFFF)) % p_;
    h = m % r_;
    q = m / r_;
}

Nooit geweten dat dit zoveel uitmaakt...


Nog steeds wel de moeite om lookup() in assembler te schrijven. 10 regels code die 60% van de tijd inneemt....hmm

[ Voor 59% gewijzigd door Zoijar op 12-07-2004 13:38 ]

[madwizard]

Zoijar schreef op 12 juli 2004 @ 13:21:

C++:

inline void LossyDictionary::hashquotient(uint32 key, uint32 a, uint32 b, uint32& h, uint32& q) const {
    uint64 mul = static_cast<uint64>(a)*key + b;
    uint32 m = ((static_cast<uint32>(mul >> 32) % p_) + static_cast<uint32>(mul & 0xFFFFFFFF)) % p_;
    h = m % r_;
    q = m / r_;
}

Deze code doet toch een hele andere berekening? Hier neem je de modulus van het high en low gedeelte apart en telt ze bij elkaar op, maar dat is niet hetzelfde als mul % p_. Of mis ik iets?

Mijn code was in assembler omdat het daar makkelijk in kan, in C weet ik niet of je hetzelfde voor elkaar kan krijgen, waarschijnlijk niet zo efficient. Maar ik zal er nog even over denken, portability is wel goed natuurlijk.

[Zoijar]

madwizard schreef op 12 juli 2004 @ 19:53:
Deze code doet toch een hele andere berekening? Hier neem je de modulus van het high en low gedeelte apart en telt ze bij elkaar op, maar dat is niet hetzelfde als mul % p_. Of mis ik iets?

Ja klopt dit was een beetje iets doms van me, ik kwam er ook later achter dat dit helemaal niet klopt. Het werkte wel omdat het een hash functie is, en er dus alsnog steeds dezelfde getallen werden berekent. Maar de clustering was erg slecht. Ik kwam erachter toen de kwaliteit ineens behoorlijk achteruit was gegaan.

Moet ik trouwens wel twee keer delen in asm? Dat het quotient overflowed is niet zo erg toch? De rest is altijd kleiner dan 32bits. Of werkt de hele div dan niet goed, en staat de rest ook niet goed?

[ Voor 15% gewijzigd door Zoijar op 12-07-2004 21:21 ]

[madwizard]

Zoijar schreef op 12 juli 2004 @ 20:35:
Moet ik trouwens wel twee keer delen in asm? Dat het quotient overflowed is niet zo erg toch? De rest is altijd kleiner dan 32bits. Of werkt de hele div dan niet goed, en staat de rest ook niet goed?

Je moet wel 2 keer delen ja want het zijn echt 2 aparte stappen. Ik heb er nog eens over nagedacht maar volgens mij is dit echt de enige oplossing, VC wil iig met geen mogelijkheid een div gebruiken. Ik heb nog gedacht ook 2^32 eerst door mod p te halen (= (-p)%p) en die dan te vermenigvuldigen met H om zo een 64-bits deling te voorkomen maar het probleem is dat het resultaat van die vermenigvuldiging ook weer groter dan 32-bits kan zijn en je dus nog geen stap verder komt. Assembler is volgens mij wel de makkelijkste en snelste oplossing.

edit: oja en het quotient overflowen is dus wel erg omdat je dan een exception krijgt .

[ Voor 6% gewijzigd door madwizard op 12-07-2004 21:31 ]

[Zoijar]

madwizard schreef op 12 juli 2004 @ 21:31:
Ik heb nog gedacht ook 2^32 eerst door mod p te halen (= (-p)%p) en die dan te vermenigvuldigen met H om zo een 64-bits deling te voorkomen maar het probleem is dat het resultaat van die vermenigvuldiging ook weer groter dan 32-bits kan zijn en je dus nog geen stap verder komt. Assembler is volgens mij wel de makkelijkste en snelste oplossing.

Ja, dacht ik ook net aan, had hier al 2^32 % p op papier staan zelfde conclusie idd.

edit: oja en het quotient overflowen is dus wel erg omdat je dan een exception krijgt .

Ja...zie het hier nu ook, m'n 386 ref manual maar weer eens opgezocht. Tijd geleden zeg...maar idd een int 0 bij overflow of deling door 0.
Ik inline asm het wel voor visual C++ en GCC, dan zou het op de meeste platforms wel moeten compilen.

[madwizard]

Ik ben alleen bang dat alleen die functie in assembler schrijven nog steeds niet veel sneller is. VC zal in ieder geval die functie vervolgens niet meer inlinen en alle parameters enzo via de stack doorgeven (weet niet of GCC nog iets slims hiervoor heeft). Zo heb je alsnog een functie call. Ik weet niet hoe je deze functie gebruikt maar als de loop eromheen niet te ingewikkeld is zou ik dat stuk ook in assembler doen, anders is waarschijnlijk de overhead van het inline asm stuk te groot om het efficient te maken. Bovendien kun je dan p_ en r_ mooi in registers blijven houden (mits ze hetzelfde blijven natuurlijk).

[Zoijar]

Het valt mee, VC optimized redelijk goed.

Deze code:

C++:

inline uint32 modmul(uint32 key, uint32 a, uint32 p) {
    _asm {
        mov ebx, key    ;
        mov eax, a      ;
        mul ebx         ;
        mov ebx, p      ;
        mov ecx, eax    ;
        mov eax, edx    ;
        xor edx, edx    ;
        div ebx         ;
        mov eax, ecx    ;
        div ebx         ;
        mov eax, edx    ;
    }
}

inline void LossyDictionary::hashquotient(uint32 key, uint32 a, uint32 b, uint32& h, uint32& q) const {
    uint32 m = (modmul(key, a, p_) + b) % p_;
    h = m % r_;
    q = m / r_;
}

Worden beide inlined in de lookup() tot het volgende (ook die 2 delingen voor rem/quo van r_ worden met 1 div gedaan)

code:

; 227  :    uint32 hash;
; 228  :    uint32 quotient;
; 229  : 
; 230  :    uint32 q;
; 231  :    sint32 d;
; 232  : 
; 233  :    hashquotient(key, a1_, b1_, hash, quotient);

    mov eax, DWORD PTR [esi+8]
    mov ecx, DWORD PTR [esi]
    push    ebx
    push    edi
    mov edi, DWORD PTR [esi+12]
    mov DWORD PTR $T78891[esp+16], eax
    mov DWORD PTR $T78883[esp+16], ecx

; hier begint de modmul code
    mov ebx, DWORD PTR _key$[esp+12]
    mov eax, DWORD PTR $T78891[esp+16]
    mul ebx
    mov ebx, DWORD PTR $T78883[esp+16]
    mov ecx, eax
    mov eax, edx
    xor edx, edx
    div ebx
    mov eax, ecx
    div ebx
    mov eax, edx
; einde modmul

; de add en laatste % p_
    mov ecx, DWORD PTR [esi]
    add eax, edi
    xor edx, edx
    div ecx

;het bepalen van / r_ en % r_
    mov edi, DWORD PTR [esi+4]
    mov eax, edx
    xor edx, edx
    div edi
    mov edi, eax

; 234  :    unpack(data_[hash], q, d);

    mov eax, DWORD PTR [esi+32]

; 241  :    if (q == (quotient&quotient_mask)) {
; 242  :        return d;
; 243  :    }
; 244  : 
; 245  :    return 0;

    movzx   eax, WORD PTR [eax+edx*2]
    mov edx, eax
    xor edx, edi
    test    dl, 7
    je  SHORT $L78931

// etc.

Anders wordt het dit. Daar mixed hij wel mooi die andere tests en unpack erdoorheen...

code:

; 227  :    uint32 hash;
; 228  :    uint32 quotient;
; 229  : 
; 230  :    uint32 q;
; 231  :    sint32 d;
; 232  : 
; 233  :    hashquotient(key, a1_, b1_, hash, quotient);

    mov ecx, DWORD PTR _key$[esp+12]
    mov eax, DWORD PTR [esi+8]
    push    ebx
    push    ebp

; 241  :    if (q == (quotient&quotient_mask)) {
; 242  :        return d;
; 243  :    }
; 244  : 
; 245  :    return 0;

    mov ebp, DWORD PTR [esi]
    push    edi
    xor edi, edi
    push    edi
    push    ecx
    push    edi
    push    eax
    call    __allmul
    mov ecx, eax
    mov eax, DWORD PTR [esi+12]
    xor ebx, ebx
    add ecx, eax
    adc edx, ebx
    xor eax, eax
    push    eax
    push    ebp
    push    edx
    push    ecx
    call    __aullrem
    mov ecx, DWORD PTR [esi+4]
    push    ebx
    push    ecx
    push    edx
    push    eax
    call    __aulldvrm
    mov DWORD PTR tv212[esp+32], ebx
    mov DWORD PTR tv72[esp+32], edx
    mov edx, DWORD PTR [esi+32]
    movzx   ecx, WORD PTR [edx+ecx*2]
    mov DWORD PTR tv280[esp+28], edx
    mov edx, ecx
    xor edx, eax
    test    dl, 7
    je  SHORT $L77843
    mov ecx, DWORD PTR _key$[esp+24]
    mov eax, DWORD PTR [esi+16]
    push    edi
    push    ecx
    push    edi
    push    eax
    call    __allmul
    mov ecx, eax
    mov eax, DWORD PTR [esi+20]
    add ecx, eax
    adc edx, edi
    xor eax, eax
    push    eax
    push    ebp
    push    edx
    push    ecx
    call    __aullrem
    mov ecx, DWORD PTR [esi+4]
    push    edi
    push    ecx
    push    edx
    push    eax
    call    __aulldvrm
    mov DWORD PTR tv203[esp+32], ebx
    mov DWORD PTR tv159[esp+32], edx
    mov edx, DWORD PTR [esi+4]
    add edx, ecx
    mov ecx, DWORD PTR tv280[esp+28]
    movzx   ecx, WORD PTR [ecx+edx*2]
    mov edx, ecx
    xor edx, eax
    test    dl, 7
    jne SHORT $L9433
$L77843:

; 234  :    unpack(data_[hash], q, d);
; 235  :    if (q == (quotient&quotient_mask)) {
; 236  :        return d;
; 237  :    }
; 238  : 
; 239  :    hashquotient(key, a2_, b2_, hash, quotient);
; 240  :    unpack(data_[hash + r_], q, d);

    mov eax, ecx
    shr eax, 3
    pop edi

; 241  :    if (q == (quotient&quotient_mask)) {
; 242  :        return d;
; 243  :    }
; 244  : 
; 245  :    return 0;

    shl eax, 19                 ; 00000013H
    pop ebp
    sar eax, 19                 ; 00000013H
    pop ebx

; 246  : }

    add esp, 16                 ; 00000010H
    ret 4
$L9433:
    pop edi
    pop ebp

; 241  :    if (q == (quotient&quotient_mask)) {
; 242  :        return d;
; 243  :    }
; 244  : 
; 245  :    return 0;

    xor eax, eax
    pop ebx

; 246  : }

    add esp, 16                 ; 00000010H
    ret 4
?lookup@LossyDictionary@PTRS@@QBEHI@Z ENDP      ; PTRS::LossyDictionary::lookup

Sorry voor de lange post...

[ Voor 47% gewijzigd door Zoijar op 12-07-2004 22:54 ]

[madwizard]

Hmm dat doet ie beter dan ik had gedacht ja, hoe snel is het nu?

[Zoijar]

madwizard schreef op 12 juli 2004 @ 23:03:
Hmm dat doet ie beter dan ik had gedacht ja, hoe snel is het nu?

Ja het is een stuk sneller. Dit is de timing nu (hoewel die niet zoveel meer zegt, omdat alles inlined is in lookup)



Ik ging ongeveer van 30 naar 38 fps zoiets met die asm code. Dat is 26% sneller. En voor dit alles was het nog trager. (begon op 9fps voor het optimizen, maar dat zat ook nog ergens anders )

[Zoijar]

Hmm eerste keer dat ik gcc asm schrijf... dit lijkt niet te werken? Het compiled wel zonder fouten, maar de data is corrupt.

C++:

inline uint32 modmul(uint32 key, uint32 a, uint32 p) {
    register uint32 result;
    __asm__ ("mull %%ebx;\n\t"
             "movl %%ecx, %%ebx;\n\t"
             "movl %%eax, %%ecx;\n\t"
             "movl %%edx, %%eax;\n\t"
             "xorl %%edx, %%edx;\n\t"
             "divl %%ebx;\n\t"
             "movl %%ecx, %%eax;\n\t"
             "divl %%ebx;\n\t"
             :"=d"(result)
             :"a"(a), "b"(key), "c"(p)
             );
    return result;
}

[RickN]

Euh, als je

(( expr ) % p ) % r

wilt uitrekenen en r is kleiner dan p, dan reken je toch gewoon

( expr ) % r

uit, eventueel met nog wat strategisch geplaatste %r's om het resultaat klein te houden....

[Zoijar]

Dat is niet hetzelfde. bv (17 % 11) % 5 = 6%5 = 1, maar 17%5 = 2.

[madwizard]

Zoijar schreef op 13 juli 2004 @ 15:21:
Hmm eerste keer dat ik gcc asm schrijf... dit lijkt niet te werken? Het compiled wel zonder fouten, maar de data is corrupt.

Ik werk er ook nooit mee (geloof dat ik alleen ooit 'int 3' in gcc asm heb geschreven), maar gebruikt ie AT&T of intel syntax (voor de volgorde van de operands)? En bewaart ie wel automatisch ebx?

[RickN]

Zoijar schreef op 13 juli 2004 @ 15:58:
Dat is niet hetzelfde. bv (17 % 11) % 5 = 6%5 = 1, maar 17%5 = 2.

Oh jezus, ik wist wel dat ik niet had moeten reageren

[Zoijar]

RickN schreef op 13 juli 2004 @ 16:06:
Oh jezus, ik wist wel dat ik niet had moeten reageren

Je bracht me wel even aan het twijfelen Hmmm ja volgens mij wordt ebx wel bewaard...de output ziet er redelijk goed uit zo op het eerste gezicht. hmmm...

[madwizard]

Zoijar schreef op 13 juli 2004 @ 16:21:
Hmmm ja volgens mij wordt ebx wel bewaard...de output ziet er redelijk goed uit zo op het eerste gezicht. hmmm...

Ik denk dat er AT&T syntax gebruikt wordt waar de operands omgedraaid worden. mov eax, ebx schrijf je dan dus andersom als movl %%ebx, %%eax.

[Zoijar]

Ja ik heb alles ook omgedraaid toch?

[madwizard]

Ohja sorry dat had ik niet gezien. Maar werkt het dan nu wel? Want eerst zei je nog dat de data corrupt was.

[Zoijar]

madwizard schreef op 13 juli 2004 @ 19:14:
Ohja sorry dat had ik niet gezien. Maar werkt het dan nu wel? Want eerst zei je nog dat de data corrupt was.

Het werkt nog steeds niet nee. Met VC dus wel, maar op de cluster met GCC niet. Compiled zonder fouten, draait ook zonder fouten, maar dan start het display programma niet op met de gegenereerde dataset.

Als ik de asm vervang door die oude int64 dingen, dan werkt het wel op de cluster. Het zijn pentium3's die solaris draaien, dus het zou moeten werken volgens mij...Er zullen toch wel ergens registers niet goed bewaard worden ofzo, of de verkeerde input...

[madwizard]

Het zou gewoon moeten werken volgens mij, volgens de documentatie regelt gcc zelf de registers. Er is nog een extra 'clobbered registers' lijst die je op kunt geven aan het asm statement maar dat is alleen voor registers die geen input of output zijn en toch gebruikt worden en die hebt je niet.

Met Dev-C++ onder windows werkt je code goed hier volgens mij.

[Zoijar]

Aha, ik heb getest met gcc en als ik zo compile gaat het goed (ik test met een paar miljoen random keys of er hetzelfde uitkomt als met de 64bit instructies) gcc modmul.cc -o modmul -lstdc++ -O3 -fno-inline En zo gaat het fout: gcc modmul.cc -o modmul -lstdc++ -O3. Dus er is duidelijk iets mis met het bewaren/restoren/output/input van registers. Als de functie inlined wordt gaat het mis met gcc.

Ok, eindelijk opgelost. zal het nog even posten hier om de thread compleet te maken. Het bleek dat gcc aanneemt dat je input operands onveranderd blijven. Maar je kan ze ook niet als clobbered opgeven, want dan denkt gcc dat je ze niet als input mag gebruiken. Dus moet je ze maar ook als output opgeven. Zit op zich wel wat in, er wordt tenslotte een waarde output, ook al hebben we totaal geen interesse in die waarde.

C++:

inline uint32 modmul(uint32 key, uint32 a, uint32 p) { 
    register uint32 result; 
    register uint32 tmp;
    __asm__ ("mull %%ebx;\n\t" 
             "movl %%ecx, %%ebx;\n\t" 
             "movl %%eax, %%ecx;\n\t" 
             "movl %%edx, %%eax;\n\t" 
             "xorl %%edx, %%edx;\n\t" 
             "divl %%ebx;\n\t" 
             "movl %%ecx, %%eax;\n\t" 
             "divl %%ebx;\n\t" 
             :"=d"(result),"=a"(tmp),"=b"(tmp),"=c"(tmp) 
             :"a"(a), "b"(key), "c"(p) 
             ); 
    return result; 
}

[ Voor 56% gewijzigd door Zoijar op 16-07-2004 13:38 ]