momente şi schiţe de informatică şi matematică
To attain knowledge, write. To attain wisdom, rewrite.

Subrutină performantă pentru calculul factorialului

factorial
2007 oct

În Hello world! ne-am ocupat şi de o subrutină amifac.s în limbaj de asamblare, pentru calculul factorialului; s-au conturat acolo câteva idei de optimizare şi le probăm aici.

Optimizări asupra algoritmului

N! s-a calculat în tabloul tw[] modelând mot-a-mot definiţia factorialului:

     iniţial tw[] = 1, F = N; 
     repetă cât timp F > 1: 
            înmulţeşte factorialul curent tw[] cu factorul F;
            F = F - 1  // următorul factor

"înmulţeşte factorialul curent tw[] cu factorul F" - înseamnă desigur, că se înmulţeşte F cu fiecare "cifră" înscrisă în tabloul tw[], propagând corespunzător "transportul".

Numărul de înmulţiri se poate reduce foarte mult, înlocuind factorii individuali cu care se înmulţeşte repetat tw[] curent, prin produs parţial de aceşti factori; de exemplu, în loc de a înmulţi tw[] întâi cu F, apoi a înmulţi rezultatul tw[] cu (F-1) şi apoi tw[] obţinut cu (F-2) - putem înmulţi tw[] direct cu produsul F*(F-1)*(F-2).

Această îmbunătăţire asupra algoritmului trebuie corelată cu posibilităţile microprocesorului în ceea ce priveşte operaţia de înmulţire. Instrucţiunea mul operand înmulţeşte valoarea existentă în registrul acumulator (AL, AX sau EAX după caz) cu operandul respectiv, producând rezultatul în "acumulatorul extins" (AH:AL, DX:AX, respectiv EDX:EAX după cum operandul este o valoare de 8 / 16 / 32 biţi); dacă operandul se află în memorie (nu într-un registru), atunci timpul necesar înmulţirii este mai mare. Pentru operanzi de 16 biţi, instrucţiunea mul se execută în 13 - 26 tacţi ("clock cycles") pe I80486 şi respectiv 11 tacţi în cazul unui Pentium; pentru operanzi de 32 biţi - în 13 - 42 tacţi pe I80486, respectiv 10 tacţi la Pentium.

Vom angaja înmulţiri de forma mull %esi: valoarea de 32 biţi existentă în EAX se înmulţeşte cu valoarea de 32 biţi din ESI şi rezultatul de 64 biţi se obţine în perechea de regiştri EDX, EAX (în 13 tacţi pe I80486, repectiv 10 pe Pentium). Produsul parţial de factori vizat mai sus trebuie să "încapă" pe 32 de biţi (în ESI).

Prin urmare, pentru obţinerea valorii binare a factorialului în tabloul tw[], avem următorul plan de lucru:

— ESI = produsul parţial curent de factori individuali.
EAX = F, factorul individual curent; mul (F-1) dă rezultatul F*(F-1) în EDX:EAX, setând flagul Carry (în registrul de flaguri al microprocesorului) dacă rezultatul depăşeşte 32 biţi; cât timp Carry == 0, coboară F (= F - 1) şi repetă mul. Când Carry devine 1, atunci produsul parţial necesar este rezultatul penultimei înmulţiri (salvează în ESI).

— EBX = adresa de bază a tabloului tw[].
Vom adresa tabloul tw[] pe cuvinte de câte 32 biţi ("cifre" de câte 4 octeţi); EBX referă prima cifră, EBX + 4 referă a doua cifră, etc.
Obs. Îmbunătăţim astfel maniera de lucru din prima versiune, unde tw[] era adresat pe cuvinte de numai 16 biţi (adică se opera în baza 2^16; acum vom modela operaţiile necesare în baza 2^32).

— ECX = rangul cifrei curente din tw[] (sărind însă peste zerourile iniţiale din tw[]).
Cifra curentă va fi copiată în acumulator prin movl (%ebx, %ecx, 4), %eax (EAX = valoarea de 32 biţi de la adresa EBX + 4*ECX).
Obs. Îmbunătăţim prima versiune: de data aceasta vom ignora cifrele zero nesemnificative (aflate în tw[] pe rangurile inferioare) - a vedea secţiunea marcată FAC11 din amifac.s redat mai jos.

— înmulţeşte cifra curentă EAX cu produsul parţial curent ESI (şi adaugă transportul precedent).

— înscrie în tw[] cifra EAX rezultată prin înmulţire şi salvează transportul EDX rezultat;
Adresează următoarea cifră (ECX ++) din tw[] şi continuă înmulţirea asupra acesteia (copiază noua cifră în EAX, înmulţeşte cu ESI, etc.).

În amifac.s redat mai jos, acest plan este realizat în secvenţele marcate PR_par_fct (obţinerea produsului parţial de factori) şi FAC1 - FAC 3.

Optimizări asupra conversiei la forma zecimală

Obţinând factorialul tw[] în forma binară (în baza 2^32), urmează să-l convertim la forma zecimală obişnuită. Experimentele au arătat că această porţiune finală a programului solicită mai mult de 90% din timpul total de execuţie. De exemplu, am inserat o secvenţă de afişare a unui mesaj "Gata" imediat după secvenţa FAC0 - FAC3 care determină forma binară a factorialului şi am lansat ./amifac 100000; s-a afişat "Gata" (adică s-a terminat calculul lui 100000!, în formă binară) după circa 5 secunde, apoi s-a afişat şi în baza 10 după încă 90 secunde (în condiţiile în care am păstrat conversia la baza 10^4, din cadrul versiunii "amifac.s" redate în Hello world!).

Putem reduce mult numărul de împărţiri (prin noua bază) necesare algoritmului de conversie, mărind baza de la 10^4 (cum folosisem anterior, când forma binară tw[] era văzută ca reprezentare în baza 2^16, adică cu cifre de 16 biţi) la 10^9; într-adevăr, se vede uşor că 10^9 este cea mai mare bază (care să fie o putere a lui 10) ale cărei cifre "încap" pe 32 biţi (avem 231 < 109 < 232 < 1010). Făcând modificările necesare (vezi secvenţele P2 - P4 şi QQ) şi efectuând în noile condiţii (baza 10^9 în loc de 10^4) experimentul pomenit mai sus, am constatat reducerea timpului necesar conversiei la 30 secunde (de la 90 secunde).

Adaptările necesare angajării bazei 10^9 în locul celei vechi 10^4 sunt relativ simple (şi n-a trebuit să rescriem secvenţa respectivă, ci doar să facem unele înlocuiri). Dificultatea a apărut din cu totul altă parte; obţinând timpi de trei ori mai buni (cum am arătat mai sus), ne-am pus evident şi problema corectitudinii rezultatului şi am constatat că acestea (cam pe oricare N!) sunt greşite!

Investigaţiile făcute (folosind GNU gdb) au arătat până la urmă că forma binară tw[] (în baza 2^32) se obţine corect prin secvenţa FAC0 - FAC3; dar ea este apoi "stricată" de către secvenţa SWAP - secvenţă pe care nu o modificasem faţă de versiunea iniţială 'amifac.s'.

SWAP interschimba o cifră de 16 biţi de la un capăt al tabloului tw[] cu cea corespunzătoare de la celălalt capăt, încât în final cifrele (de câte 16 biţi) din tw[] ajungeau în ordinea lor uzuală (de la cifra cea mai semnificativă pe rangul 0, spre cifra "unităţilor" pe rangul ultim).

Însă a interschimba cifre de 16 biţi (cum fusese necesar când angajasem operaţii în baza 2^16) nu este totuna cu a interschimba cifre de 32 biţi (cum este necesar acum, când se operează în baza 2^32). Prin urmare, am modificat (cum se vede în programul de mai jos) şi secvenţa SWAP, încât ea să asigure interschimbarea de cuvinte de câte 32 biţi.

Optimizări privind operaţia de împărţire

Implicând baza 10^9, am redus foarte mult numărul total de împărţiri necesare (faţă de cazul când foloseam baza 10^4). Iarăşi, trebuie să corelăm această îmbunătăţire (de nivel algoritmic) cu posibilităţile microprocesorului în ceea ce priveşte operaţia de împărţire.

Instrucţiunea div operand împarte conţinutul acumulatorului (sau acumulatorului extins) la operand şi furnizează câtul şi restul în acumulatorul extins; necesită cel puţin 16 / 24 / 40 tacţi pe I80486 şi 17 / 25 / 41 tacţi pe Pentium, după cum împărţitorul are 8 / 16 / 32 biţi (numai una sau două instrucţiuni din setul standard de instrucţiuni Pentium, sunt mai costisitoare decât div sau idiv).

## Exemplu de folosire a instrucţiunii DIV ##
movl High, %edx       # EDX = partea High a deîmpărţitului 
                      #       (= 0, pentru deîmpărţit de 32 biţi)
movl Low, %eax        # EAX = partea Low a deîmpărţitului; 
                      # EDX:EAX = deîmpărţitul (de 64 biţi)
movl divisor, %ebx    # EBX = împărţitor 
                      # (poate fi utilizat orice registru disponibil)
divl %ebx             # rezultă EAX = câtul (EDX:EAX / EBX) şi 
                      #         EDX = restul (EDX:EAX % EBX)

Se poate obţine un câştig important de viteză, evitând folosirea instrucţiunii div. Ideea constă în a înlocui împărţirea la operand prin înmulţirea cu inversul operandului (mul fiind mai rapidă ca div).

De exemplu, pentru a obţine câtul şi restul împărţirii unei valori de 32 biţi prin 10 avem următoarea secvenţă, mult mai rapidă decât div:

   movl $0xCCCCCCCD, %edi    # EDI = 10^(-1) deplasat stânga cu 3 biţi
   movl (%esi), %eax         # preia deîmpărţitul (fie y) din memorie, în EAX 
   movl %eax, %ebx           # salvează y în EBX (necesar pentru calculul restului)
   mull %edi                 # EDX:EAX = y * 10^(-1) 
                             #           adică EDX = câtul deplasat stânga 3 biţi
   shr $3, %edx              # deplasează invers 3 biţi# EDX = câtul împărţirii (fie q)
   leal (%edx,%edx,4), %eax  # EAX = EDX + 4*EDX = 5*q
   addl %eax, %eax           # EAX = 10 * q
   subl %eax, %ebx           # EBX = y - 10*q = restul împărţirii

Socotind mot-a-mot (pentru Pentium) mul necesită 10 tacţi, iar celelalte 7 instrucţiuni implicate necesită câte 1 tact. Evident, în cazul când împărţirile trebuie repetate asupra câtului obţinut (ca în contextul conversiei la baza 10), numărul total de tacţi se va diminua, datorită scoaterii în afara ciclului a primelor două instrucţiuni.

De unde provine valoarea 0xCCCCCCCD? Inversul lui 10 este 0.1 în baza 10, dar în baza 2 este fracţia 0.0001100110011... (de perioadă 0011). Neglijând cele trei zerouri de după virgulă (adică deplasând stânga cu 3 biţi), considerând numai primii 32 biţi ai fracţiei rezultate şi ignorând "virgula", obţinem numărul întreg de 32 de biţi 11001100...11001100 (şi avem 1100 = 8 + 4 = 0xC), care aproximează prin lipsă fracţia respectivă.
Rotunjind în sus (ca să ţinem cont în calcul şi de biţii 1100... care urmează după cei 32 biţi reţinuţi din fracţie), obţinem 0xCCCCCCCD. Când se înmulţeşte y cu această valoare (y*0.1) se obţine valoarea "rotunjită" (şi deplasată spre stânga cu 3 biţi) a câtului y/10; eliminând cei 3 biţi finali (prin deplasarea rezultatului spre dreapta) obţinem valoarea exactă a câtului întreg.

Pare improbabil să putem folosi o secvenţă analogă celei de mai sus şi pentru împărţirile la 10^9 (de pe parcursul conversiei factorialului tw[] din baza 2^32, la baza 10^9). Avem nevoie nu numai de cât, dar şi de rest - ori restul se va determina (vezi liniile 7-9 în secvenţa de mai sus) scăzând din deîmpărţitul iniţial produsul dintre câtul obţinut şi împărţitorul 10^9; această înmulţire a câtului s-a putut înlocui prin adunări în cazul când împărţitorul era 10, dar dacă împărţitorul este 10^9 atunci pare obligatoriu de folosit mul (timpul necesar ar creşte la 27 tacţi, ceea ce încă ar fi bine faţă de div - numai că apar necesităţi suplimentare, de exemplu salvarea temporară şi recuperarea valorilor din regiştrii implicaţi EDX, EAX; în total n-ar fi mai avantajos, faţă de folosirea directă a instrucţiunii div).

Varianta finală amifac.s

De data aceasta punem toate subrutinele în acelaşi fişier (în versiunea din "Hello world!" folosisem fişiere separate şi un Makefile); asamblarea şi obţinerea fişierului executabil se face deci prin:

   vb@debian:~/lar/ASM as  -o amifac.o  amifac.s
   vb@debian:~/lar/ASM ld  amifac  -o amifac.o

Sau mai simplu, scriem aceste comenzi într-un fişier "Makefile" şi folosim make pentru a le executa (după fiecare modificare).

.data
        zq9: .byte 0x30,0x30,0x30,0x30,0x30,0x30,0x30,0x30,0x30
.align 4
        nw: .long 1     # numărul de cifre (de 32 biţi) ale factorialului
        tw: .long 0     # adresa reprezentării binare a factorialului
        Nxt: .long 0    # factorul următor (pentru produsul parţial viitor)

.text   
.align 8        

    ## subrutina principală:  vb@debian:~/lar/ASM/ ./amifac 20000   (N = "20000")
    ## rezultă N! (factorialul)  # dar se presupune N > 2     
.globl _start   
_start: 
   movl 8(%esp), %edi   # EDI = adresa şirului de cifre N 
   call edi_toi         # valoarea numerică (32 biţi) a lui N
   movl %eax, Nxt       # Nxt = factorul următor, iniţial N 
   pushl %eax           
   call lg_fact         # (ESP) = numărul de cifre 2^32 pentru tw[]
   call me_mmap         # alocă memorie pentru tw[] (4*(ESP) + 4 Bytes)
   movl %eax, tw        
   movl Nxt, %edi
   movl %edi, (%eax)    # iniţializează factorialul tw[] cu N
   popl %eax            
   decl Nxt             # se pp. că N > 2
   
   movl tw, %ebx        # fixează EBX pe adresa de bază a tw[], cu acces 
                        # la cifra de rang ECX prin (EBX + 4*ECX)
PR_par_fct:
        # următorul produs parţial de factori Nxt*(Nxt-1)*... (< 2^32)
   movl Nxt, %edi       # EDI = factorul până la care s-a coborât
   movl %edi, %eax
D1:
   movl %eax, %esi      # salvează ESI = produsul parţial curent (< 2^32)
   subl $1, %edi
   cmpl $1, %edi
   je D2
   mul %edi             # EAX = produsul parţial de factori EDI*(EDI-1)*...
   jnc D1               # extinde încă p.p., dacă el încă este < 2^32
D2:
   movl %edi, Nxt       # salvează factorul până la care s-a coborât
                        # ESI = produsul parţial curent (de factori)
FAC1:
   xorl %ebp, %ebp      # EBP va servi pentru propagarea transporturilor
   movl $0xFFFFFFFF, %ecx
FAC11:                  # EBX + 4*ECX va referi cifra curentă (32 biţi) din tw[],
   incl %ecx            # evitând însă cifrele zero de la începutul lui tw[]
   movl (%ebx,%ecx,4), %eax
   testl %eax, %eax
   jz FAC11             # dacă cifra EAX == 0, reia ECX++
FAC2:
   movl (%ebx,%ecx,4), %eax     # EAX = cifra curentă (32 biţi) din tw[]
   mull %esi                    # EDX:EAX = pr_par_fct * cifra curentă   
   addl %ebp, %eax              # şi adaugă transportul
   adcl $0, %edx        # propagă eventualul transport al adunării precedente
   movl %edx, %ebp      # EBP = noul transport
   movl %eax, (%ebx,%ecx,4)     # înscrie în tw[] cifra curentă obţinută 
   incl %ecx            # referă următoarea cifră (de 32 biţi) din tw[]
   cmpl nw, %ecx
   jne FAC2             # reia, dacă în tw[] mai există cifre neoperate

   test %ebp, %ebp      # există transport final?
   je FAC3
   movl %ebp, (%ebx,%ecx,4)     # înscrie transportul final în tw[] şi
   incl nw                      # incrementează contorul de cifre (32 biţi) nw
   
FAC3:
   cmpl $1, Nxt
   jne PR_par_fct       # "reapelează" factorial(n-1) (dacă Nxt > 1)

        # inversează tw[] ('swap' între stânga şi dreapta)
   movl tw, %edi        # EDI referă primul BYTE din tw[]
   movl nw, %esi
   leal -4(%edi,%esi,4), %esi   # ESI referă ultima cifră din tw[] (EDI+4*nw-4)
SWAP:   # (EDI) <--> (ESI)# EDI += 4 şi ESI -= 4# până când EDI >= ESI
   movl (%edi), %eax    # EAX = (EDI) 
   xchgl (%esi), %eax   # (ESI) = EAX, EAX = (ESI)
   stosl                # (EDI) = EAX şi EDI += 4
   subl $4, %esi        # ESI -= 4
   cmpl %esi, %edi      # repetă interschimbarea dacă EDI < ESI
   jb SWAP

## converteşte tw[] la baza 10^9# z9[] (cu cifre de 32 biţi) ##
   pushl nw
   call lg_nw_z9        # estimarea numărului de cifre 10^9 (de 32 biţi)
   call me_mmap         # alocă memorie pentru z9[]
   popl %ebx

   movl %eax, %edi      # EDI = adresa de bază a tabloului z9[] alocat
   movl nw, %ecx        # ECX = nw (va contoriza 'nw' pt. deîmpărţitul tw[] curent)
   movl tw, %esi        # ESI referă deîmpărţitul tw[] curent, de la prima cifră (HIGH)
   movl $1000000000, %ebx       # se va împărţi prin noua bază, 10^9
   xorl %ebp, %ebp      # EBP va contoriza cifrele înregistrate în z9[]
P2:
   pushl %ecx           # determină câtul şi restul împărţirii (schema lui Horner)
   pushl %esi           # lui tw[] (deîmpărţitul "curent"), la noua bază 10^9
   xorl %edx, %edx      # EDX = restul curent (care înmulţit cu vechea bază 2^32 şi
                        # adunat cu cifra curentă EAX = deîmpărţitul parţial curent)
P3:
   movl (%esi), %eax    # deîmpărţitul parţial = EDX:EAX = restul_curent * 2^32 + EAX
   divl %ebx            # rezultă EDX = restul_curent şi EAX = câtul curent 
   movl %eax, (%esi)    # câtul curent intră în constituţia viitorului deîmpărţit tw[]
   addl $4, %esi
   loop P3
   movl %edx, (%edi)    # înscrie restul final al împărţirii curente a lui tw[] prin 10^9 
   addl $4, %edi
   incl %ebp            # contorizează cifrele înscrise
   popl %esi            # reconstituie adresa de bază a deîmpărţitului curent tw[]
   popl %ecx            # şi lungimea lui
P4:
   movl (%esi), %eax    # testează dacă au apărut zerouri nesemnificative
   testl %eax, %eax     #          în deîmpărţitul tw[]
   jnz P2
   addl $4, %esi        # ignoră cifrele zero iniţiale din deîmpărţitul curent
   decl %ecx            # corectând lungimea deîmpărţitului curent tw[]
   jnz P4

## scrie în formă zecimală cifrele (32 biţi) din z9[] ##
   subl $4, %edi
   movl %edi, %esi              # ESI referă penultimul BYTE din z9[]
   movl %ebp, %ecx              # ECX = numărul cifrelor (de 32 biţi) din z9[]
   movl $0xCCCCCCCD, %edi       # EDI = 10^(-1) deplasat stânga cu 3 biţi
   
QQ:
   movl $zq9, %ebp      # EBP referă zona $zq9, unde se vor înscrie cifrele zecimale
   movl $0x30303030, (%ebp)     # iniţializează $zq9 cu 9 cifre '0' (cod ASCII 0x30)
   movl $0x30303030, 4(%ebp)
   movb $0x30, 8(%ebp)
   addl $9, %ebp        # EBP va adresa $zq9 începând de la cifra unităţilor
   movl (%esi), %ebx    # cifra curentă din z9[], de convertit la cele 9 cifre zecimale 
Q1:
   movl %ebx, %eax           # EAX = deîmpărţitul curent = y (câtul precedentei împărţiri)
   mull %edi                 # EDX:EAX = y * 10^(-1) adică EDX = câtul deplasat stânga
   shr $3, %edx              # deplasează invers# EDX = câtul împărţirii = q
   leal (%edx,%edx,4), %eax  # EAX = EDX + 4*EDX = 5*EDX
   addl %eax, %eax           # EAX = 10 * EDX
   subl %eax, %ebx           # EBX = y - 10*q = restul împărţirii
   addb $48, %bl        # transformă ASCII restul BL
   decl %ebp
   movb %bl, (%ebp)     # înscrie cifra zecimală
   movl %edx, %ebx      # salvează câtul curent în EBX
   testl %edx, %edx     # repetă împărţirea (lui EBX la 10) cât timp câtul curent > 0
   jnz Q1               

   pushl %ecx           # scrie cele 9 cifre zecimale înregistrate în zona $zq9
   movl $SYS_write, %eax        # EAX = numărul serviciului 'write'
   movl $STDOUT, %ebx           # EBX = descriptorul fişierului în care se scrie
   movl $zq9, %ecx              # ECX = adresa sursei (datele de scris)
   movl $9, %edx                # EDX = lungimea sursei
   int $0x80            # apelează la kernel pentru serviciul EAX
   popl %ecx

   subl $4, %esi        # repetă pentru cifra următoare din z9[]
   loop QQ

   movl $SYS_exit, %eax # apelează serviciul de încheiere a procesului
   xorl %ebx, %ebx
   int $0x80

################ subrutine auxiliare şi valori constante ################

.globl edi_toi  ## De la şir de cifre ASCII-Z (cu adresa în EDI) la "unsigned long" ##
edi_toi:        # EDI = adresa şirului cifrelor lui N. Returnează EAX = N.
   push %ebx    # salvează pe stivă EBX, ECX, EDX
   push %ecx
   push %edx
   movl $10, %ecx       # ECX = 10
   xorl %edx, %edx      # iniţializează EDX = EAX = EBX = 0
   xorl %eax, %eax
   xorl %ebx, %ebx
ed1:
   movb (%edi), %bl     # BL = cifra ASCII pointată de EDI
   test %bl, %bl        # un şir ASCII-Z se încheie cu '\0' 
   jz ed2               # dacă BL=='\0' intră la adresa 'ed2'
   subl $48, %ebx       # BL = BL - '0'
   mul %ecx             # EDX:EAX = EAX * 10 + EBX
   addl %ebx, %eax
   inc %edi             # pointează următoarea cifră
   jmp ed1              # reia ("jump") de la adresa 'ed1'
ed2:
   popl %edx    # reconstituie EDX, ECX, EBX de pe stivă
   popl %ecx
   popl %ebx
ret             

.globl lg_fact  ## numărul cifre (32 biţi) ale reprezentării factorialului în baza 2^32 ##
lg_fact:           # (ESP+4) conţine iniţial N şi în final numărul cifre 2^32 pentru N! 
        fldz       # iniţial suma_lg(k!) este 0.0 (k=N..2)
qont:   
        fldlg2          # st(0)=lg(2), st(1)=suma_lg(i!) cu i=N..k-1
        fildl 4(%esp)   # st(0)=k (factorul curent, N..2), st(1)=lg(2), st(2)=suma_lg(i!)
        fyl2x           # st(0)=lg(k), st(1)=suma_lg(i!) cu i=N..k-1
        faddp %st, %st(1)       # st(0)=suma_lg(i!) cu i=N..k
        decl 4(%esp)    # următorul factor N..2 (k--)
        cmpl $1, 4(%esp)        
        jne qont        # reia, dacă k>1
                        # în final st(0) = suma_lg(k!) cu k=2..N
        fldlg2          # st(0) = lg(2), st(1) = suma_lg(k!)
        fadd %st, %st   # st(0) = 2*lg(2), st(1)=suma_lg(k!)
        fadd %st, %st   # st(0) = 4*lg(2), st(1)=suma_lg(k!)
        fadd %st, %st
        fadd %st, %st   # st(0) = 16*lg(2) = lg(2^16), st(1)=suma_lg(k!)
        fadd %st, %st   # st(0) = 32*lg(2) = lg(2^32), st(1)=suma_lg(k!)
        fxch            # interschimbă st(0) cu st(1)
        fdivp %st, %st(1)       # st(0) = suma_lg(k!) / lg(2^32) = nw =
        fistpl 4(%esp)          # = numărul de "cifre" în baza 2^32 pentru N!
ret

.globl lg_nw_z9 ## numărul cifre (32 biţi) ale reprezentării factorialului în baza 10^9 ##
                        # (ESP+4) conţine nw = numărul cifrelor 2^32 din tw[]
lg_nw_z9:               # se determina numărul de cifre 10^9 corespunzător lui nw
        fldz            # lungimea necesară rezultă din 10^(9*Z) = 2^(32*nw)
        fldlg2          # st(0) = lg(2)
        fadd %st, %st   # st(0) = 2*lg(2)
        fadd %st, %st   # st(0) = 4*lg(2)
        fildl 4(%esp)   # st(0) = nw, st(1) = 
        fmul %st(1), %st  # st(0) = nw * lg(16) = nr. dword-uri necesare
        fistpl 4(%esp)    # reprezentării 10^9, pentru tw[] (32*nw*lg(2) / 9 = 4*nw*lg(2))
        movl 4(%esp), %eax      # returnează rezultatul în EAX
ret

.globl me_mmap  ## apelează mmap(address, length, protect, flags, filedes, offset)# ##
me_mmap:
   movl %esp, %ebp
   subl $24, %esp       # mmap() va prelua cei 6 parametri din stivă

   movl 4(%ebp), %eax   # numărul de dWord-uri de alocat
   addl $2, %eax        # eventual, cu 2 mai mult decât s-a estimat
   sall $2, %eax        # un dWord = 4 BYTES
                        # se vor aloca EAX bytes, apelând mmap()
   xorl %edx,%edx
   movl %edx, (%esp)    # NULL (adresa de bază va fi alocată de sistem)
   movl %eax, 4(%esp)   # dimensiunea fişierului (în Bytes) 
   movl $3, 8(%esp)     # PROT_READ | PROT_WRITE
   movl $0x22, 12(%esp) # MAP_PRIVATE | MAP_ANONYMOUS
   movl $-1, 16(%esp)   # descriptorul fişierului (-1 sau 0, pentru anonim)
   movl %edx, 20(%esp)  # offest_2 (zero)
   
   movl $SYS_mmap, %eax 
   movl %esp, %ebx      # EBX referă cadrul stivă care conţine cei 6 parametri
   int  $0x80           # mmap() returnează adresa zonei alocate în EAX

   movl %ebp, %esp
ret

## apeluri de sistem utilizate ##
   .equ SYS_exit, 1     # numărul apelului de sistem pentru exit()
   .equ SYS_write, 4
   .equ STDOUT, 1       # descriptorul de fişier al ieşirii standard (la 'stdout')
   .equ SYS_mmap, 90    # pentru apelarea functiei de alocare mmap()

Fişierul sursă amifac.s are 250 rânduri (inclusiv liniile albe şi comentariile), dar fişierul executabil rezultat amifac măsoară doar 2.3 kB (n-am mai optat la asamblare, pentru includerea de informaţie necesară depanării).

Timpii de execuţie obţinuţi sunt (sau par a fi) onorabili: de la "instantaneu" până la 7-8 secunde pentru valori mai mici de 60000!; cam 30 secunde pentru valoarea lui 100000!

Timpii rezultaţi printr-un program C (compilat cu GNU gcc) care implică pentru obţinerea factorialului funcţia mpz_fac_ui() din biblioteca GMP (a vedea "Hello world!", secţiunea Programul de factoriale, folosind GMP) sunt totuşi, cam de două ori mai buni; vreo două motive sunt clare: în GMP factorii n, n-1, n-2, ..., 2 sunt grupaţi în produse parţiale de tipul celor care apar la descompunerea în factori primi (2^k * 3^h *...) iar produsul lor (pentru obţinerea finală a factorialului) se face "divide et impera" şi încă, folosind un algoritm de înmulţire (între numere mari) mai performant decât cel obişnuit folosit de noi (şi se folosesc numai deplasări în loc de înmulţiri, pentru factorul 2^k; se foloseşte un tablou precalculat cu valorile factorialelor până la 31!); în plus, se reuşeşte optimizarea împărţirilor (evitând complet div, ceea ce noi n-am reuşit aici).

Porţiunea din amifac.s care ar mai trebui optimizată (şi în principal este vorba de evitarea div) este aceea care face conversia de la forma binară a factorialului (rezultată în tw[] după ce s-a parcurs secvenţa FAC1 - FAC3 şi SWAP), la reprezentarea în baza 10^9 (adică secţiunea P2 - P4). Într-adevăr, inserând după SWAP o secvenţă suplimentară de scriere a unui mesaj şi executând programul pentru 100000!, am putut vedea că tw[] în forma binară se obţine în 4 - 5 secunde, în timp ce rezultatul în forma zecimală este afişat după încă 25 secunde.

Dar apare o întrebare firească: de ce să fim obligaţi (adesea) să optimizăm prin soft operaţii "elementare" precum div? Normal ar fi ca ele să fi fost realizate în cel mai bun mod posibil, direct la nivelul circuitelor hardware!

Realitatea este aceasta: div, chiar mul sunt operaţii rar folosite în programare, faţă de instrucţiunile de transfer de exemplu; un microprocesor care să facă în mod optim operaţii ca div (fără a mai necesita optimizări soft) ar costa probabil dublu (sau chiar triplu) - iar facilitatea de a face împărţiri în cel mai rapid mod posibil ar fi utilă unui număr destul de restrâns de programatori (pentru marea majoritate a utilizatorilor div-ul existent acum este chiar prea bun).

Pe de altă parte… mai trebuie lăsată apă la moară şi programatorilor (softişti)!

vezi Cărţile mele (de programare)

docerpro | Prev | Next