RISOLTO Eseguire ricerca in un algoritmo

Pubblicità
Stato
Discussione chiusa ad ulteriori risposte.
- seguendo invece l'idea che ho proposto nei precedenti post:

C++:
#include <iostream>
#include <cstdlib>
#include <ctime>

using namespace std;

int main()
{
    srand(time(0));

    const int MIN = 38;
    const int MAX = 55;
    const int DIM_1 = 12;
    const int DIM_2 = MAX - MIN + 1;

    int v[DIM_1];
    int u[DIM_2] = {0};

    for(int i = 0; i < DIM_1; ++i)
    {
        v[i] = rand() % DIM_2 + MIN;
        u[v[i] - MIN] = 1;
        cout << v[i] << " ";
    }
    cout << endl;
    for(int i = MIN; i <= MAX; ++i)
    {
        if(u[i - MIN] == 0)
        {
            cout << i << " ";
        }
    }
}

A sto punto anche direttamente così, senza il vettore v, che non mi sembra sia tra i requisiti:

C++:
#include <iostream>
#include <cstdlib>
#include <ctime>

using namespace std;

int main()
{
    srand(time(0));
    const int MIN = 38;
    const int MAX = 55;
    const int DIM_1 = 12;
    const int DIM_2 = MAX - MIN + 1;

    int u[DIM_2] = {0};

    for(int i = 0; i < DIM_1; ++i)
    {
        int v  = rand() % DIM_2 + MIN;
        u[v - MIN] = 1;
        cout << v << " ";
    }

    cout << endl;

    for(int i = MIN; i <= MAX; ++i)
    {
        if(!u[i - MIN])
        {
            cout << i << " " ;
        }
    }
}

Tralaltro interessante vedere come viene trasformato quel rand() % DIM_2 + MIN pur di non usare la divisione:

Codice:
        call    rand
        mov     edi, OFFSET FLAT:std::cout
        movsx   rsi, eax
        cdq
        imul    rsi, rsi, 954437177
        sar     rsi, 34
        sub     esi, edx
        lea     edx, [rsi+rsi*8]
        add     edx, edx
        sub     eax, edx
        mov     esi, eax
        cdqe
        add     esi, 38
        mov     DWORD PTR [rsp+rax*4], 1

ottimo, o anche in C

Codice:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

static int v[56];

int main()
{
        int i, c = 0;

        srand(time(0));

        do {
                i = rand() % 100;
                if (i < 38 || i > 55)
                        continue;
                v[i]++;
                c++;
        } while (c < 10);

        for (i = 38; i <= 55; i++) {
                if (!v[i])
                        printf("il numero %d non e' nella serie\n", i);
        }

        return 0;
}

Brutto però quell'if li dentro... 😁 farà molti cicli inutili. Non ho capito perchè stai usando 100 anche, tanto i valori sopra al 55 vengono scartati comunque. Mi sono perso qualcosa?
 
A sto punto anche direttamente così, senza il vettore v, che non mi sembra sia tra i requisiti:
Nel post iniziale viene menzionato, ma per quello che deve fare il programma è effettivamente superfluo.

Tralaltro interessante vedere come viene trasformato quel rand() % DIM_2 + MIN pur di non usare la divisione: ...
Che tradotto per noi umani sarebbe? 😅
 
A sto punto anche direttamente così, senza il vettore v, che non mi sembra sia tra i requisiti:

C++:
#include <iostream>
#include <cstdlib>
#include <ctime>

using namespace std;

int main()
{
    srand(time(0));
    const int MIN = 38;
    const int MAX = 55;
    const int DIM_1 = 12;
    const int DIM_2 = MAX - MIN + 1;

    int u[DIM_2] = {0};

    for(int i = 0; i < DIM_1; ++i)
    {
        int v  = rand() % DIM_2 + MIN;
        u[v - MIN] = 1;
        cout << v << " ";
    }

    cout << endl;

    for(int i = MIN; i <= MAX; ++i)
    {
        if(!u[i - MIN])
        {
            cout << i << " " ;
        }
    }
}

Tralaltro interessante vedere come viene trasformato quel rand() % DIM_2 + MIN pur di non usare la divisione:

Codice:
        call    rand
        mov     edi, OFFSET FLAT:std::cout
        movsx   rsi, eax
        cdq
        imul    rsi, rsi, 954437177
        sar     rsi, 34
        sub     esi, edx
        lea     edx, [rsi+rsi*8]
        add     edx, edx
        sub     eax, edx
        mov     esi, eax
        cdqe
        add     esi, 38
        mov     DWORD PTR [rsp+rax*4], 1



Brutto però quell'if li dentro... 😁 farà molti cicli inutili. Non ho capito perchè stai usando 100 anche, tanto i valori sopra al 55 vengono scartati comunque. Mi sono perso qualcosa?
importante che funzioni, era un misto, %100 era per velocizzare, cicli sprecati in un moderno pc non sono un dramma, ma era cosi, per scriverlo in poche righe., Tanto avra' gia consegnato il compitino al prof, fatto benino anche da lui.
 
Che tradotto per noi umani sarebbe? 😅

Giusta osservazione 😛

Codice:
        call    rand
        mov     edi, OFFSET FLAT:std::cout
        movsx   rsi, eax
        cdq
        imul    rsi, rsi, 954437177
        sar     rsi, 34
        sub     esi, edx
        lea     edx, [rsi+rsi*8]
        add     edx, edx
        sub     eax, edx
        mov     esi, eax
        cdqe
        add     esi, 38
        mov     DWORD PTR [rsp+rax*4], 1

Provo a scomporlo in parti (salto ciò che possiamo "ignorare").

Il risultato della rand (il nr generato insomma) viene piazzato in EAX che a sua volta viene spostato in RSI:

Codice:
 movsx   rsi, eax
 cdq
 imul    rsi, rsi, 954437177

rsi viene poi moltiplicato, in pratica è: RSI = RSI * 954437177. La costante è sicuramente qualcosa di noto, è un valore intermedio usato per velocizzare il calcolo.

Codice:
        sar     rsi, 34
        sub     esi, edx

il risultato ottenuto in precedenza viene shiftato di 34 bit verso destra. Ogni spostamento di 1 bit equivale a una divisione per 2.
La particolarità di SAR è che mantiene anche il bit del segno; quindi se era negativo prima dello shift, viene reinserito il bit del segno a sx.
Viene poi presa solo la parte bassa di RSI, cioè ESI, e avviene un ESI = ESI - EDX (in teoria la parte alta dovrebbe comunque non essere valorizzata visto lo shift di prima verso dx).
Il valore in EDX dipende dal bit del segno di EAX; la modifica è stata fatta dall'istruzione CDQ che si vede li sopra.

Codice:
        lea     edx, [rsi+rsi*8]
        add     edx, edx
        sub     eax, edx
        mov     esi, eax
        cdqe
        add     esi, 38

il valore appena calcolato in ESI (che sono i 32bit bassi di RSI) viene poi usato per calcolare RSI + RSI*8; questo risultato viene memorizzato in EDX che viene poi nell'istruzione sotto sommato di nuovo a EDX (eg. ha appunto evitato anche qui di moltiplicare).

Questo risultato finale viene sottratto al valore di EAX che viene memorizzato in ESI; a ESI viene sommato 38, mentre EAX contiene di fatto l'indice dell'array, la posizione alla quale accedere. In pratica il contenuto di v[i] è in EAX.
Mi sfugge a prima vista perchè avvenga quell' ADD ESI + 38, mi sembra non venga usato e faccia già probabilmente parte del calcolo precedente... ma dovrò guardare meglio.

Di sicuro non viene usato il vettore "v" però.
Per altro, il vettore "u" è in realtà molto più grande di come scritto nel codice (il valore DIM_2). Nel programma compilato si vede infatti:

Codice:
        pxor    xmm0, xmm0
        mov     QWORD PTR [rsp+64], 0
        movaps  XMMWORD PTR [rsp], xmm0
        movaps  XMMWORD PTR [rsp+16], xmm0
        movaps  XMMWORD PTR [rsp+32], xmm0
        movaps  XMMWORD PTR [rsp+48], xmm0

che è l'inizializzazione dell'array.

Il codice completo disassemblato è questo, per chi fosse interessato:

Codice:
.LC0:
        .string " "
main:
        push    rbx
        xor     edi, edi
        mov     ebx, 12
        sub     rsp, 80
        call    time
        mov     edi, eax
        call    srand
        pxor    xmm0, xmm0
        mov     QWORD PTR [rsp+64], 0
        movaps  XMMWORD PTR [rsp], xmm0
        movaps  XMMWORD PTR [rsp+16], xmm0
        movaps  XMMWORD PTR [rsp+32], xmm0
        movaps  XMMWORD PTR [rsp+48], xmm0
.L4:
        call    rand
        mov     edi, OFFSET FLAT:std::cout
        movsx   rsi, eax
        cdq
        imul    rsi, rsi, 954437177
        sar     rsi, 34
        sub     esi, edx
        lea     edx, [rsi+rsi*8]
        add     edx, edx
        sub     eax, edx
        mov     esi, eax
        cdqe
        add     esi, 38
        mov     DWORD PTR [rsp+rax*4], 1
        call    std::basic_ostream<char, std::char_traits<char> >::operator<<(int)
        mov     edx, 1
        mov     esi, OFFSET FLAT:.LC0
        mov     rdi, rax
        call    std::basic_ostream<char, std::char_traits<char> >& std::__ostream_insert<char, std::char_traits<char> >(std::basic_ostream<char, std::char_traits<char> >&, char const*, long)
        sub     ebx, 1
        jne     .L4
        mov     rax, QWORD PTR std::cout[rip]
        mov     rax, QWORD PTR [rax-24]
        mov     rbx, QWORD PTR std::cout[rax+240]
        test    rbx, rbx
        je      .L15
        cmp     BYTE PTR [rbx+56], 0
        je      .L6
        movsx   esi, BYTE PTR [rbx+67]
.L7:
        mov     edi, OFFSET FLAT:std::cout
        mov     ebx, 38
        call    std::basic_ostream<char, std::char_traits<char> >::put(char)
        mov     rdi, rax
        call    std::basic_ostream<char, std::char_traits<char> >::flush()
        jmp     .L9
.L8:
        add     rbx, 1
        cmp     rbx, 56
        je      .L16
.L9:
        mov     eax, DWORD PTR [rsp-152+rbx*4]
        test    eax, eax
        jne     .L8
        mov     esi, ebx
        mov     edi, OFFSET FLAT:std::cout
        add     rbx, 1
        call    std::basic_ostream<char, std::char_traits<char> >::operator<<(int)
        mov     edx, 1
        mov     esi, OFFSET FLAT:.LC0
        mov     rdi, rax
        call    std::basic_ostream<char, std::char_traits<char> >& std::__ostream_insert<char, std::char_traits<char> >(std::basic_ostream<char, std::char_traits<char> >&, char const*, long)
        cmp     rbx, 56
        jne     .L9
.L16:
        add     rsp, 80
        xor     eax, eax
        pop     rbx
        ret
.L6:
        mov     rdi, rbx
        call    std::ctype<char>::_M_widen_init() const
        mov     rax, QWORD PTR [rbx]
        mov     esi, 10
        mov     rax, QWORD PTR [rax+48]
        cmp     rax, OFFSET FLAT:_ZNKSt5ctypeIcE8do_widenEc
        je      .L7
        mov     esi, 10
        mov     rdi, rbx
        call    rax
        movsx   esi, al
        jmp     .L7
.L15:
        call    std::__throw_bad_cast()

Da debugger risulterebbe meno arcano comprendere 😂 e almeno si vedrebbe il contenuto della memoria e dei registri.

Io ho fatto la prova online usando -O2 per compilare. Senza ottimizzazioni è anche più divertente ancora... 🤣

Codice:
.LC0:
        .string " "
main:
        push    rbp
        mov     rbp, rsp
        sub     rsp, 160
        mov     edi, 0
        call    time
        mov     edi, eax
        call    srand
        mov     DWORD PTR [rbp-12], 38
        mov     DWORD PTR [rbp-16], 55
        mov     DWORD PTR [rbp-20], 12
        mov     DWORD PTR [rbp-24], 18
        pxor    xmm0, xmm0
        movaps  XMMWORD PTR [rbp-160], xmm0
        movaps  XMMWORD PTR [rbp-144], xmm0
        movaps  XMMWORD PTR [rbp-128], xmm0
        movaps  XMMWORD PTR [rbp-112], xmm0
        movq    QWORD PTR [rbp-96], xmm0
        mov     DWORD PTR [rbp-4], 0
        jmp     .L2
.L3:
        call    rand
        mov     ecx, eax
        movsx   rax, ecx
        imul    rax, rax, 954437177
        shr     rax, 32
        mov     edx, eax
        sar     edx, 2
        mov     eax, ecx
        sar     eax, 31
        sub     edx, eax
        mov     eax, edx
        sal     eax, 3
        add     eax, edx
        add     eax, eax
        sub     ecx, eax
        mov     edx, ecx
        add     edx, 38
        mov     eax, DWORD PTR [rbp-4]
        cdqe
        mov     DWORD PTR [rbp-80+rax*4], edx
        mov     eax, DWORD PTR [rbp-4]
        cdqe
        mov     eax, DWORD PTR [rbp-80+rax*4]
        sub     eax, 38
        cdqe
        mov     DWORD PTR [rbp-160+rax*4], 1
        mov     eax, DWORD PTR [rbp-4]
        cdqe
        mov     eax, DWORD PTR [rbp-80+rax*4]
        mov     esi, eax
        mov     edi, OFFSET FLAT:_ZSt4cout
        call    std::basic_ostream<char, std::char_traits<char> >::operator<<(int)
        mov     esi, OFFSET FLAT:.LC0
        mov     rdi, rax
        call    std::basic_ostream<char, std::char_traits<char> >& std::operator<< <std::char_traits<char> >(std::basic_ostream<char, std::char_traits<char> >&, char const*)
        add     DWORD PTR [rbp-4], 1
.L2:
        cmp     DWORD PTR [rbp-4], 11
        jle     .L3
        mov     esi, OFFSET FLAT:_ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_
        mov     edi, OFFSET FLAT:_ZSt4cout
        call    std::basic_ostream<char, std::char_traits<char> >::operator<<(std::basic_ostream<char, std::char_traits<char> >& (*)(std::basic_ostream<char, std::char_traits<char> >&))
        mov     DWORD PTR [rbp-8], 38
        jmp     .L4
.L6:
        mov     eax, DWORD PTR [rbp-8]
        sub     eax, 38
        cdqe
        mov     eax, DWORD PTR [rbp-160+rax*4]
        test    eax, eax
        jne     .L5
        mov     eax, DWORD PTR [rbp-8]
        mov     esi, eax
        mov     edi, OFFSET FLAT:_ZSt4cout
        call    std::basic_ostream<char, std::char_traits<char> >::operator<<(int)
        mov     esi, OFFSET FLAT:.LC0
        mov     rdi, rax
        call    std::basic_ostream<char, std::char_traits<char> >& std::operator<< <std::char_traits<char> >(std::basic_ostream<char, std::char_traits<char> >&, char const*)
.L5:
        add     DWORD PTR [rbp-8], 1
.L4:
        cmp     DWORD PTR [rbp-8], 55
        jle     .L6
        mov     eax, 0
        leave
        ret


un misto, %100 era per velocizzare, cicli sprecati in un moderno pc non sono un dramma, ma era cosi, per scriverlo in poche righe., Tanto avra' gia consegnato il compitino al prof, fatto benino anche da lui.

No era più un'osservazione riguardante il tempo perso ciclando; vero che sono poche esecuzioni. Ma se la probabilità di pescare un valore da 1 a 100 è la medesima, ci sono l'80% (arrotondo) di probabilità di non averlo nel range. Giusto?
 
l'ho buttato giu in 2 minuti, onestamente non mi son posto troppi quesiti accademici :)

Ho solo valutato al vol che rand da fino a un minimo di 32767, quindi ho usato il modulo. Poi ovviamente, si, due cicli su 3, o 3 su 4 circa andranno sprecati, che per 10 numeri in un moderno pc e' un tempo trascurabile. Mi piaceva ridurre il codice al massimo.
 
rsi viene poi moltiplicato, in pratica è: RSI = RSI * 954437177. La costante è sicuramente qualcosa di noto, è un valore intermedio usato per velocizzare il calcolo...
Visto che mi incuriosiva come cosa, ho deciso di rifletterci un po' su. Alla fine non ho chiarito la questione del tutto, ma qualcosina penso di averlo capito.

Innanzitutto ho pensato che matematicamente per "trasformare" una divisione in una moltiplicazione si può fare qualcosa del genere:

A/B = A*(1/B)

Ho provato poi ad effettuare la divisione 1/9 in binario:

1/9 = 0.000111 = 0.000111000111000111...

divisione.png

Posto

A = 6589 = 1100110111101
B = 9 = 1001

avremo quindi che

A/B = A*(1/B) = 6589*(1/9) = 1100110111101*(1/1001) = 1100110111101 * 0.000111

Dal momento che 1/9 non può essere rappresentato in binario in modo esatto, e tenendo presente il discorso da te fatto sui 32 bit alti e bassi, ho pensato di approssimare la parte dopo la virgola in 32 bit:

0.000111 0.000(11100011100011100011100011100011) = .11100011100011100011100011100011 / 2^3

Dalla classica moltiplicazione in colonna che si fa a scuola sappiamo che si moltiplicano i due fattori ignorando le virgole e poi alla fine si divide il risultato per 10 elevato al numero di cifre dopo le virgole; similmente la parte intera della suddetta moltiplicazione (che poi è anche la parte intera Q della nostra divisione) si può calcolare come:

Q = 1100110111101 * 11100011100011100011100011100011 >> (3+32) = 1011011100 = 732

Essendo

1/18 = (1/9)/2

tornando al caso del topic si ricava:

Q = A * 11100011100011100011100011100011 >> (3+32+1) =
= A * 11100011100011100011100011100011 >> 36 =
= A * (11100011100011100011100011100011 >> 2) >> 34 =
= A * 111000111000111000111000111000 >> 34

Ed ecco che spuntano il 34 e il numero

111000111000111000111000111000 = 954437176

che differisce di un'unità dal 954437177 sopra quotato.
Quasi sicuramente deve esserci qualche particolare che mi sfugge e di cui non ho tenuto conto nella fase di approssimazione, ma come disse un noto virologo: "funzionicchia"! 😁
Se avete qualche idea al proposito, fatemi sapere.


P.S.1
Se la parte dopo la virgola di 1/B occupa 32 bit, anche A dovrà occupare al massimo 32 bit, affinché appunto il risultato della moltiplicazione possa essere espresso attraverso una variabile intera a 64 bit.

P.S.2
Relativamente al funzionamento di queste ottimizzazioni, ipotizzo che il compilatore abbia già a disposizione una sorta di tabella che ad ogni numero dispari B associ un intero che approssimi 1/B. Che ne pensate?
 
Ultima modifica:
Visto che mi incuriosiva come cosa, ho deciso di rifletterci un po' su. Alla fine non ho chiarito la questione del tutto, ma qualcosina penso di averlo capito.

Innanzitutto ho pensato che matematicamente per "trasformare" una divisione in una moltiplicazione si può fare qualcosa del genere:

A/B = A*(1/B)

Ho provato poi ad effettuare la divisione 1/9 in binario:

1/9 = 0.000111 = 0.000111000111000111...

divisione.png

Posto

A = 6589 = 1100110111101
B = 9 = 1001

avremo quindi che

A/B = A*(1/B) = 6589*(1/9) = 1100110111101*(1/1001) = 1100110111101 * 0.000111

Dal momento che 1/9 non può essere rappresentato in binario in modo esatto, e tenendo presente il discorso da te fatto sui 32 bit alti e bassi, ho pensato di approssimare la parte dopo la virgola in 32 bit:

0.000111 0.000(11100011100011100011100011100011) = .11100011100011100011100011100011 / 2^3

Dalla classica moltiplicazione in colonna che si fa a scuola sappiamo che si moltiplicano i due fattori ignorando le virgole e poi alla fine si divide il risultato per 10 elevato al numero di cifre dopo le virgole; similmente la parte intera della suddetta moltiplicazione (che poi è anche la parte intera Q della nostra divisione) si può calcolare come:

Q = 1100110111101 * 11100011100011100011100011100011 >> (3+32) = 1011011100 = 732

Essendo

1/18 = (1/9)/2

tornando al caso del topic si ricava:

Q = A * 11100011100011100011100011100011 >> (3+32+1) =
= A * 11100011100011100011100011100011 >> 36 =
= A * (11100011100011100011100011100011 >> 2) >> 34 =
= A * 111000111000111000111000111000 >> 34

Ed ecco che spuntano il 34 e il numero

111000111000111000111000111000 = 954437176

che differisce di un'unità dal 954437177 sopra quotato.
Quasi sicuramente deve esserci qualche particolare che mi sfugge e di cui non ho tenuto conto nella fase di approssimazione, ma come disse un noto virologo: "funzionicchia"! 😁
Se avete qualche idea al proposito, fatemi sapere.


P.S.1
Se la parte dopo la virgola di 1/B occupa 32 bit, anche A dovrà occupare al massimo 32 bit, affinché appunto il risultato della moltiplicazione possa essere espresso attraverso una variabile intera a 64 bit.

P.S.2
Relativamente al funzionamento di queste ottimizzazioni, ipotizzo che il compilatore abbia già a disposizione una sorta di tabella che ad ogni numero dispari B associ un intero che approssimi 1/B. Che ne pensate?

Il ragionamento mi sembra filare alla grande, domani lo rileggo bene (ho finito di lavorare alle 19.40, cercando informazioni nello stack con gdb...).

In merito al PS2, si molto probabilmente si. So anche di librerie che si occupano di questi bei magheggi comunque: https://libdivide.com/ Se guardi il file header che viene linkato, vedrai tante belle cose.
Sicuramente il compilatore nella fase di ottimizzazione si occupa di generare del codice il meglio ottimizzato possibile anche in base all'architettura; quindi è probabile che qualcosa avvenga solo in quel momento, e non tutto sia come dire "harcodato" nel compilatore.
Posso provare a raccogliere qualche info in più su GCC comunque.

Ho loggato anche per linkarti/linkare questi due "topic":

- suboptimal 'division by constant' optimization: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=28417
- Consider improving division and modulo by constant if highpart multiply is cheap https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89845

Non ho provato a leggere più di tanto perchè non ce la faccio più oggi, però mi sembra ci siano info interessanti. Come noti sono riferite a versioni molto vecchie di GCC, ma almeno danno una buona idea.
 
Non ho ancora chiarito del tutto la questione dal punto di vista matematico (quindi non sono sicuro che funzioni al 100% sempre), ma nel frattempo ho implementato il procedimento descritto nel mio precedente post:

C:
#include <stdio.h>
#include <stdint.h>

uint32_t optimized_division(const uint32_t D, uint32_t d)
{
    uint32_t D_1 = 1;
    uint8_t n = 0;
    for(; !(D_1 & d); ++n, d >>= 1);
    if(d == 1)
    {
        return D >> n;
    }
    for(; D_1 < d; ++n, D_1 <<= 1);
    uint32_t q = 0;
    for(uint8_t i = 0; i < 32; ++i)
    {
        q <<= 1;
        if(D_1 > d)
        {
            q |= 1;
            D_1 -= d;
        }
        D_1 <<= 1;
    }
    return (uint64_t)D * q >> 31 + n;
}

int main()
{
    uint32_t D = 941787252;
    uint32_t d = 36248;

    printf("%u / %u = %u  DIVISIONE FUNZIONE\n", D, d, optimized_division(D, d));
    printf("%u / %u = %u  DIVISIONE NATIVA\n", D, d, D / d);
}

Ovviamente l'ottimizzazione si riferisce al fatto che, una volta noti i valori di q e n per ogni possibile valore assunto dal divisore d, la divisione si riduce di fatto all'ultima riga della funzione.
 
Bisognerebbe provare a testare i tempi usando un loop.
Ho provato al volo a creare due sorgenti diversi, uno per ciascuna operazione, per richiamarli usando time (sono sotto Linux) e i tempi sono praticamente identici.

Per avere un risultato un pò più affidabile meglio leggere in input i due valori che sono ora nel codice, mi riferisco ovviamente a D e d.

Il compilatore ora come ora sta barando se si compila con -O2:

Codice:
main:
        sub     rsp, 8
        mov     esi, 36248
        mov     edi, 941787252
        call    optimized_division
        mov     edx, 36248
        mov     esi, 941787252
        mov     edi, OFFSET FLAT:.LC0
        mov     ecx, eax
        xor     eax, eax
        call    printf
        mov     ecx, 25981
        mov     edx, 36248
        xor     eax, eax
        mov     esi, 941787252
        mov     edi, OFFSET FLAT:.LC1
        call    printf
        xor     eax, eax
        add     rsp, 8
        ret

Come noti c'è un 25981 harcodato in ECX. E' il risultato che calcola lui in compile-time, così deve solo stamparlo senza fare la divisione vera.

Senza ottimizzazione invece esegue davvero una DIV.

Codice:
        mov     eax, DWORD PTR [rbp-4]
        mov     edx, 0
        div     DWORD PTR [rbp-8]
        mov     ecx, eax

Se leggi con una scanf, almeno non può "barare".
 
Ricapitolando... dal ragionamento fatto nel mio precedente post, si deduce che la divisione intera D/d (con d < D ) si riduce a:
C:
q ? (uint64_t)D * q >> 31 + n : D >> n
con q e n che dipendono esclusivamente dal valore di d.

Ovviamente, affinché abbia senso parlare di ottimizzazione, i coefficienti q e n devono essere noti per ogni valore assumibile dal divisore d.
Nel caso in esame ho posto d_max = 100000 e ho utilizzato gli array v_q e v_n per contenere rispettivamente i valori di q e n calcolati (dalla funzione fun() ) per valori di d appartenenti all'intervallo [1;100000]. Quindi per esempio per accedere ai valori di q e n per d = 41 basterà accedere rispettivamente agli elementi v_q[41] e v_n[41].

Per testare il tutto ho disabilitato le ottimizzazioni del compilatore e ho misurato i tempi necessari a calcolare N volte una stessa divisione in diversi scenari:

- NATIVA: si riferisce alla divisione nativa del compilatore;
- OTTIMIZZATA CON ARRAY: si riferisce alla versione da me implementata con i valori di q e n prelevati rispettivamente dagli array v_q e v_n. Inoltre non ho utilizzato la funzione optimized_division() per velocizzare il tutto evitando le varie chiamate a funzione (visto che la direttiva inline non sembra funzionare);
- OTTIMIZZATA SENZA ARRAY: come sopra, ma i valori di q e n sono stati estratti dai due array per risparmiare sui tempi di accesso agli array.

Un altro test sulla velocità l'ho fatto misurando i tempi necessari a calcolare N divisioni in cui il divisore è dato da rand() + 1 (che nel mio caso va bene essendo RAND_MAX < d_max ).

Infine ho eseguito un test simile al precedente, ma questo volta non per calcolare i tempi, ma per valutare se la divisione nativa e la divisione da me implementata restituissero sempre gli stessi risultati.

Di seguito riporto il codice e l'output che ottengo lanciandolo:

C:
#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>
#include <time.h>

#define MAX_d 100000

uint32_t v_q[MAX_d + 1];
uint32_t v_n[MAX_d + 1];

void fun(const uint32_t max_d)
{
    for(uint32_t d = 2; d <= max_d; ++d)
    {
        uint32_t D = 1;
        uint32_t d_2 = d;
        uint32_t q = 0;
        uint8_t n = 0;
        for(; !(D & d_2); ++n, d_2 >>= 1);
        if(d == d_2)
        {
            for(; D < d_2; ++n, D <<= 1);
            for(uint8_t i = 0; i < 32; ++i)
            {
                q <<= 1;
                if(D > d_2)
                {
                    q |= 1;
                    D -= d_2;
                }
                D <<= 1;
            }
            v_q[d] = q;
            v_n[d] = n;
        }
        else
        {
            v_q[d] = v_q[d_2];
            v_n[d] = v_n[d_2] + n;
        }
    }
}

uint32_t optimized_division(const uint32_t D, const uint32_t d)
{
    if(d < D)
    {
        return v_q[d] ? (uint64_t)D * v_q[d] >> 31 + v_n[d] : D >> v_n[d];
    }
    return d == D ? 1 : 0;
}

int main()
{
    uint32_t D = 941787252;
    uint32_t d = 36248; // <= MAX_d

    fun(MAX_d);
    printf("OK\n\n");

    uint32_t N = 50000000;
    srand(time(0));
    clock_t begin;
    clock_t end;

    begin = clock();
    for(uint32_t i = 0; i < N; ++i)
    {
        D / d;
    }
    end = clock();
    printf("%ux [ %u / %u = %u ] IN %fs (NATIVA)\n", N, D, d, D / d, (double)(end - begin) / CLOCKS_PER_SEC);

    begin = clock();
    for(uint32_t q = v_q[d], n = v_n[d], i = 0; i < N; ++i)
    {
        d < D ? (q ? (uint64_t)D * q >> 31 + n : D >> n) : (d == D ? 1 : 0);
    }
    end = clock();
    printf("%ux [ %u / %u = %u ] IN %fs (OTTIMIZZATA SENZA ARRAY)\n", N, D, d, d < D ? (v_q[d] ? (uint64_t)D * v_q[d] >> 31 + v_n[d] : D >> v_n[d]) : (d == D ? 1 : 0), (double)(end - begin) / CLOCKS_PER_SEC);

    begin = clock();
    for(uint32_t i = 0; i < N; ++i)
    {
        d < D ? (v_q[d] ? (uint64_t)D * v_q[d] >> 31 + v_n[d] : D >> v_n[d]) : (d == D ? 1 : 0);
    }
    end = clock();
    printf("%ux [ %u / %u = %u ] IN %fs (OTTIMIZZATA CON ARRAY)\n", N, D, d, d < D ? (v_q[d] ? (uint64_t)D * v_q[d] >> 31 + v_n[d] : D >> v_n[d]) : (d == D ? 1 : 0), (double)(end - begin) / CLOCKS_PER_SEC);

    begin = clock();
    for(uint32_t i = 0; i < N; ++i)
    {
        D / (rand() + 1);
    }
    end = clock();
    printf("\n%ux [ %u / (rand() + 1) ] IN %fs (NATIVA)\n", N, D, (double)(end - begin) / CLOCKS_PER_SEC);

    begin = clock();
    for(uint32_t i = 0; i < N; ++i)
    {
        d = rand() + 1;
        d < D ? (v_q[d] ? (uint64_t)D * v_q[d] >> 31 + v_n[d] : D >> v_n[d]) : (d == D ? 1 : 0);
    }
    end = clock();
    printf("%ux [ %u / (rand() + 1) ] IN %fs (OTTIMIZZATA CON ARRAY)\n", N, D, (double)(end - begin) / CLOCKS_PER_SEC);

    printf("\n------");
    for(uint32_t cont = 0, i = 0; i < N && cont < 5; ++i)
    {
        d = rand() + 1;
        if(D / d != optimized_division(D, d))
        {
            ++cont;
            printf("\n%u / %u = %u (NATIVA)", D, d, D / d);
            printf("\n%u / %u = %u (OTTIMIZZATA)\n", D, d, optimized_division(D, d));
        }
    }
    printf("------\n");
}

Codice:
OK

50000000x [ 941787252 / 36248 = 25981 ] IN 0.134000s (NATIVA)
50000000x [ 941787252 / 36248 = 25981 ] IN 0.124000s (OTTIMIZZATA SENZA ARRAY)
50000000x [ 941787252 / 36248 = 25981 ] IN 0.140000s (OTTIMIZZATA CON ARRAY)

50000000x [ 941787252 / (rand() + 1) ] IN 0.711000s (NATIVA)
50000000x [ 941787252 / (rand() + 1) ] IN 0.782000s (OTTIMIZZATA CON ARRAY)

------
941787252 / 882 = 1067786 (NATIVA)
941787252 / 882 = 1067785 (OTTIMIZZATA)

941787252 / 147 = 6406716 (NATIVA)
941787252 / 147 = 6406715 (OTTIMIZZATA)

941787252 / 98 = 9610074 (NATIVA)
941787252 / 98 = 9610073 (OTTIMIZZATA)

941787252 / 882 = 1067786 (NATIVA)
941787252 / 882 = 1067785 (OTTIMIZZATA)

941787252 / 42 = 22423506 (NATIVA)
941787252 / 42 = 22423505 (OTTIMIZZATA)
------

Process returned 0 (0x0)   execution time : 1.931 s
Press any key to continue.


Volendo tirare le somme:

- ho disabilitato le varie ottimizzazioni del compilatore per fare in modo che quest'ultimo effettuasse realmente ogni divisione e ogni iterazione nei cicli. Ipotizzo però che anche con le ottimizzazione disabilitate il compilatore effettui la divisione intera nativa con una moltiplicazione e uno shift come mostrato da @DispatchCode in questo post, o sbaglio?

- Dall'output si osserva che relativamente ai tempi:
Codice:
OTTIMIZZATA CON ARRAY > NATIVA > OTTIMIZZATA SENZA ARRAY
Il che mi fa pensare che il compilatore prelevi comunque i coefficienti da qualche parte, ma che questo prelievo sia più veloce rispetto a quello fatto dagli array globali v_q e v_n. Che ne pensate al riguardo?

- Dal terzo test si osserva che in alcuni casi (ossia per alcuni valori di d ), la divisione da me implementata restituisce valori che sono più piccoli di un'unità rispetto al risultato reale. Come già ipotizzato in precedenza deve esserci qualche particolare che mi sfugge e di cui non ho tenuto conto nella fase di approssimazione del coefficiente q, il quale dovrebbe essere più grande di un'unità.
Magari ci rifletto un po' su, ma se avete qualche idea al proposito fatemi sapere.
 
Ultima modifica:
Gran bel lavoro, provo intanto a risponderti a qualcosa.

Ricapitolando... dal ragionamento fatto nel mio precedente post, si deduce che la divisione intera D/d (con d < D ) si riduce a:
C:
q ? (uint64_t)D * q >> 31 + n : D >> n
con q e n che dipendono esclusivamente dal valore di d.

Ovviamente, affinché abbia senso parlare di ottimizzazione, i coefficienti q e n devono essere noti per ogni valore assumibile dal divisore d.
Nel caso in esame ho posto d_max = 100000 e ho utilizzato gli array v_q e v_n per contenere rispettivamente i valori di q e n calcolati (dalla funzione fun() ) per valori di d appartenenti all'intervallo [1;100000]. Quindi per esempio per accedere ai valori di q e n per d = 41 basterà accedere rispettivamente agli elementi v_q[41] e v_n[41].

Per testare il tutto ho disabilitato le ottimizzazioni del compilatore e ho misurato i tempi necessari a calcolare N volte una stessa divisione in diversi scenari:

- NATIVA: si riferisce alla divisione nativa del compilatore;
- OTTIMIZZATA CON ARRAY: si riferisce alla versione da me implementata con i valori di q e n prelevati rispettivamente dagli array v_q e v_n. Inoltre non ho utilizzato la funzione optimized_division() per velocizzare il tutto evitando le varie chiamate a funzione (visto che la direttiva inline non sembra funzionare);
- OTTIMIZZATA SENZA ARRAY: come sopra, ma i valori di q e n sono stati estratti dai due array per risparmiare sui tempi di accesso agli array.

In che senso non sembra funzionare, usi sempre GCC o altro?
Se ti dice "undefined optimized_division" e qualcosa di analogo, devi solo dichiarare la funzione in cima:

Codice:
uint32_t optimized_division(const uint32_t D, const uint32_t d);

e poi quando definisci la funzione la indichi inline, eg:

Codice:
inline uint32_t optimized_division(const uint32_t D, const uint32_t d)
{
......

Quando dici "disabilitando le ottimizzazioni", che flag hai usato? Perchè non mi sembra che GCC abbia qualcosa, se non -O0. E anche così se non erro non disabilita completamente tutto.

Un altro test sulla velocità l'ho fatto misurando i tempi necessari a calcolare N divisioni in cui il divisore è dato da rand() + 1 (che nel mio caso va bene essendo RAND_MAX < d_max ).

Infine ho eseguito un test simile al precedente, ma questo volta non per calcolare i tempi, ma per valutare se la divisione nativa e la divisione da me implementata restituissero sempre gli stessi risultati.

Di seguito riporto il codice e l'output che ottengo lanciandolo:

C:
#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>
#include <time.h>

#define MAX_d 100000

uint32_t v_q[MAX_d + 1];
uint32_t v_n[MAX_d + 1];

void fun(const uint32_t max_d)
{
    for(uint32_t d = 2; d <= max_d; ++d)
    {
        uint32_t D = 1;
        uint32_t d_2 = d;
        uint32_t q = 0;
        uint8_t n = 0;
        for(; !(D & d_2); ++n, d_2 >>= 1);
        if(d == d_2)
        {
            for(; D < d_2; ++n, D <<= 1);
            for(uint8_t i = 0; i < 32; ++i)
            {
                q <<= 1;
                if(D > d_2)
                {
                    q |= 1;
                    D -= d_2;
                }
                D <<= 1;
            }
            v_q[d] = q;
            v_n[d] = n;
        }
        else
        {
            v_q[d] = v_q[d_2];
            v_n[d] = v_n[d_2] + n;
        }
    }
}

uint32_t optimized_division(const uint32_t D, const uint32_t d)
{
    if(d < D)
    {
        return v_q[d] ? (uint64_t)D * v_q[d] >> 31 + v_n[d] : D >> v_n[d];
    }
    return d == D ? 1 : 0;
}

int main()
{
    uint32_t D = 941787252;
    uint32_t d = 36248; // <= MAX_d

    fun(MAX_d);
    printf("OK\n\n");

    uint32_t N = 50000000;
    srand(time(0));
    clock_t begin;
    clock_t end;

    begin = clock();
    for(uint32_t i = 0; i < N; ++i)
    {
        D / d;
    }
    end = clock();
    printf("%ux [ %u / %u = %u ] IN %fs (NATIVA)\n", N, D, d, D / d, (double)(end - begin) / CLOCKS_PER_SEC);

    begin = clock();
    for(uint32_t q = v_q[d], n = v_n[d], i = 0; i < N; ++i)
    {
        d < D ? (q ? (uint64_t)D * q >> 31 + n : D >> n) : (d == D ? 1 : 0);
    }
    end = clock();
    printf("%ux [ %u / %u = %u ] IN %fs (OTTIMIZZATA SENZA ARRAY)\n", N, D, d, d < D ? (v_q[d] ? (uint64_t)D * v_q[d] >> 31 + v_n[d] : D >> v_n[d]) : (d == D ? 1 : 0), (double)(end - begin) / CLOCKS_PER_SEC);

    begin = clock();
    for(uint32_t i = 0; i < N; ++i)
    {
        d < D ? (v_q[d] ? (uint64_t)D * v_q[d] >> 31 + v_n[d] : D >> v_n[d]) : (d == D ? 1 : 0);
    }
    end = clock();
    printf("%ux [ %u / %u = %u ] IN %fs (OTTIMIZZATA CON ARRAY)\n", N, D, d, d < D ? (v_q[d] ? (uint64_t)D * v_q[d] >> 31 + v_n[d] : D >> v_n[d]) : (d == D ? 1 : 0), (double)(end - begin) / CLOCKS_PER_SEC);

    begin = clock();
    for(uint32_t i = 0; i < N; ++i)
    {
        D / (rand() + 1);
    }
    end = clock();
    printf("\n%ux [ %u / (rand() + 1) ] IN %fs (NATIVA)\n", N, D, (double)(end - begin) / CLOCKS_PER_SEC);

    begin = clock();
    for(uint32_t i = 0; i < N; ++i)
    {
        d = rand() + 1;
        d < D ? (v_q[d] ? (uint64_t)D * v_q[d] >> 31 + v_n[d] : D >> v_n[d]) : (d == D ? 1 : 0);
    }
    end = clock();
    printf("%ux [ %u / (rand() + 1) ] IN %fs (OTTIMIZZATA CON ARRAY)\n", N, D, (double)(end - begin) / CLOCKS_PER_SEC);

    printf("\n------");
    for(uint32_t cont = 0, i = 0; i < N && cont < 5; ++i)
    {
        d = rand() + 1;
        if(D / d != optimized_division(D, d))
        {
            ++cont;
            printf("\n%u / %u = %u (NATIVA)", D, d, D / d);
            printf("\n%u / %u = %u (OTTIMIZZATA)\n", D, d, optimized_division(D, d));
        }
    }
    printf("------\n");
}

Codice:
OK

50000000x [ 941787252 / 36248 = 25981 ] IN 0.134000s (NATIVA)
50000000x [ 941787252 / 36248 = 25981 ] IN 0.124000s (OTTIMIZZATA SENZA ARRAY)
50000000x [ 941787252 / 36248 = 25981 ] IN 0.140000s (OTTIMIZZATA CON ARRAY)

50000000x [ 941787252 / (rand() + 1) ] IN 0.711000s (NATIVA)
50000000x [ 941787252 / (rand() + 1) ] IN 0.782000s (OTTIMIZZATA CON ARRAY)

------
941787252 / 882 = 1067786 (NATIVA)
941787252 / 882 = 1067785 (OTTIMIZZATA)

941787252 / 147 = 6406716 (NATIVA)
941787252 / 147 = 6406715 (OTTIMIZZATA)

941787252 / 98 = 9610074 (NATIVA)
941787252 / 98 = 9610073 (OTTIMIZZATA)

941787252 / 882 = 1067786 (NATIVA)
941787252 / 882 = 1067785 (OTTIMIZZATA)

941787252 / 42 = 22423506 (NATIVA)
941787252 / 42 = 22423505 (OTTIMIZZATA)
------

Process returned 0 (0x0)   execution time : 1.931 s
Press any key to continue.
Su che OS lo stai runnando? Io sono su Linux al momento, però mi va in SIGSEGV.
Ho guardato al volo con gdb, il problema è qui:

Codice:
Core was generated by `./es'.
Program terminated with signal SIGSEGV, Segmentation fault.
#0 0x00000000004015e4 in main () at es.c:104
104             d < D ? (v_q[d] ? (uint64_t)D * v_q[d] >> 31 + v_n[d] : D >> v_n[d]) : (d == D ? 1 : 0);

Ho visto che il problema c'è anche compilando da godbolt, quindi presumo sia qualcosa che fa il compilatore.

Che versione stai usando (e che compilatore, lo so, chiedo ancora xD)?

Volendo tirare le somme:

- ho disabilitato le varie ottimizzazioni del compilatore per fare in modo che quest'ultimo effettuasse realmente ogni divisione e ogni iterazione nei cicli. Ipotizzo però che anche con le ottimizzazione disabilitate il compilatore effettui la divisione intera nativa con una moltiplicazione e uno shift come mostrato da @DispatchCode in questo post, o sbaglio?

In generale, se non usi -O2 / -O3, pare che GCC non produca gran ottimizzazioni con le divisioni, a quanto vedo. Se vuoi limitarle ulteriormente leggi in input i due valori che hai harcodato nel codice, così limiti ulteriormente quello che può fare non conoscendo i valori.

Nel caso del post che hai linkato l'operazione era il resto, quindi il compilatore ha ottimizzato comunque. Dipende dai casi, dipende da come si scrivono le cose e... il compilatore. Nel precedente caso era C++, possibile che si comporti diversamente rispetto a C, trattandosi di g++.

- Dall'output si osserva che relativamente ai tempi:
Codice:
OTTIMIZZATA CON ARRAY > NATIVA > OTTIMIZZATA SENZA ARRAY
Il che mi fa pensare che il compilatore prelevi comunque i coefficienti da qualche parte, ma che questo prelievo sia più veloce rispetto a quello fatto dagli array globali v_q e v_n. Che ne pensate al riguardo?

Sarà più veloce per una serie di fattori. Per esempio, l'array che hai dichiarato è di 100000 elementi, quindi occupa 400000 bytes. Questo implica che anche se continui ad accedere non sarà contenuto tutto nella cache; se accedi a posizioni distante e "randomiche" (non prevedibili insomma), nella cache ci saranno solo gli elementi adiacenti.

Ci sarebbero poi considerazione sulla memoria virtuale stessa che possono avere un impatto: è vero che tu accedi in O(1) agli array; ma viene trascurato un tempo che in realtà non è così trascurabile, ovvero gli accessi alla memoria che avvengono se il dato a cui accedi (o l'indirizzo dell'istruzione stessa... anche se penso non sia il caso, essendo molto piccolo come codice) non è in RAM. In tal caso il primo accesso richiede di percorrere la tabella delle pagine, che è su 4 o 5 livelli.

Quindi se sommi tutti questi ritardi può essere un motivo. A questo va aggiunto che viene calcolato sicuramente l'indice dell'array, quindi non sarà così immediato. Non ho avuto tempo di guardarmi il disasm ancora, l'ho buttato al volo su godbolt:

Codice:
v_q[d] = q;

in asm lo fa come:

Codice:
        mov     eax, DWORD PTR [rbp-4]
        mov     edx, DWORD PTR [rbp-16]
        mov     DWORD PTR v_q[0+rax*4], edx

quindi già solo qui ci sono 3 accessi alla memoria: 1 per leggere in rbp-4, 1 per leggere in rbp-16 e l'altro per calcolare 0+rax*4.

Se indichi "uint32_t d" come register, for(register uint32_t d = 2; d <= max_d; ++d), ottieni:

Codice:
        mov     edx, ebx
        mov     eax, DWORD PTR [rbp-20]
        mov     DWORD PTR v_q[0+rdx*4], eax

quindi risparmi (BEN lol) 1 accesso alla memoria.

Mi fermo qui al momento, verrò a rileggerti però appena ho del tempo e a riguardare il codice. 🙂
 
In che senso non sembra funzionare, usi sempre GCC o altro?
Se ti dice "undefined optimized_division" e qualcosa di analogo, devi solo dichiarare la funzione in cima:
Sì, uso GCC.
All'inizio non mi compilava a causa di un errore di undefined reference to 'optimized_division', ma avevo risolto utilizzando static inline (e a tal proposito ho constatato che anche il metodo che mi hai suggerito, di aggiungere il prototipo della funzione, fa scomparire quell'errore). In ogni caso con "non sembra funzionare" intendevo il fatto che i tempi con e senza inline sono gli stessi, quindi alla fine per eliminare i tempi delle chiamate a funzione ho tagliato la testa al toro andando a sostituire le chiamate nel main() con il corpo della funzione.


Quando dici "disabilitando le ottimizzazioni", che flag hai usato? Perchè non mi sembra che GCC abbia qualcosa, se non -O0. E anche così se non erro non disabilita completamente tutto.
Non saprei, mi sono limitato a disabilitare le spunte nella finestra delle ottimizzazioni del compilatore di Code::Blocks.


Io sono su Linux al momento, però mi va in SIGSEGV.
Non so, potrebbe essere dovuto a questo?
Un altro test sulla velocità l'ho fatto misurando i tempi necessari a calcolare N divisioni in cui il divisore è dato da rand() + 1 (che nel mio caso va bene essendo RAND_MAX + 1 < d_max ).


In generale, se non usi -O2 / -O3, pare che GCC non produca gran ottimizzazioni con le divisioni, a quanto vedo. Se vuoi limitarle ulteriormente leggi in input i due valori che hai harcodato nel codice, così limiti ulteriormente quello che può fare non conoscendo i valori.
Non conosco l'assembly, ma compilando con -O2/-O3 i cicli for nel main() venivano ignorati anche leggendo D e d da input. Per questo alla fine ho disabilitato tutti i flag relativi alle ottimizzazioni.
Alla fine lo scopo di quei test sulla velocità non è quello di battere il compilatore, anzi è quello di ottenere dei tempi sovrapponibili che mi lascino pensare che il procedimento che ho implementato è uguale, o comunque simile, a quello utilizzato dal compilatore.
Ovviamente affinché i test fatti abbiamo senso, anche la divisione nativa deve essere calcolata attraverso una moltiplicazione e uno shift come nel caso che ha aperto questa nostra discussione e che io ho cercato di capire e riprodurre... e dai tempi sono abbastanza convinto che l'algoritmo di fondo sia più o meno lo stesso.


Sarà più veloce per una serie di fattori...
Come già detto più volte, l'ottimizzazione del metodo che ho implementato ha senso solo se i coefficienti q e n per ogni possibile divisore d risultano già noti, al fine appunto di poter "trasformare" una divisione in un moltiplicazione e in uno shift. Nel mio caso calcolo questi coefficienti attraverso la funzione fun() e li salvo negli array v_q e v_n.
Detto ciò la domanda è la seguente: concordi sul fatto che anche il compilatore deve avere già a disposizione da qualche parte questi valori? E in tal caso dove e come ci accede?
Questo spiegherebbe per esempio che nonostante i tempi siano simili risulta:

t(OTTIMIZZATA CON ARRAY) > t(NATIVA) > t(OTTIMIZZATA SENZA ARRAY)
 
Stato
Discussione chiusa ad ulteriori risposte.
Pubblicità
Pubblicità
Indietro
Top