[C++] Meccaniche segrete del compilatore

[Zero-Giulio]

Ok, premetto che non sono un grande esperto di programmazione.
Magari quello che mi succede è normalissimo e nessuno di voi si stupirebbe. Ma a me sembra inspiegabile e quindi vi chiedo lumi.

Guardate un attimo queste 5 righe di codice:

double value = 0.0;
value += (ptr->_pdd * ptr->_ndd->_temp3);
value += (ptr->_pdu * ptr->_ndu->_temp3);
value += (ptr->_pud * ptr->_nud->_temp3);
value += (ptr->_puu * ptr->_nuu->_temp3);

Niente di che. La variabile double value viene inizializzata a zero, e successivamente incrementata 4 volte.

Ora guardate queste 3 righe di codice:

double temp = 0.0;
double ttt = (ptr->_pdd * ptr->_ndd->_temp3)
              + (ptr->_pdu * ptr->_ndu->_temp3)
              + (ptr->_pud * ptr->_nud->_temp3)
              + (ptr->_puu * ptr->_nuu->_temp3);
temp = ttt;

Anche qui, nulla di strano. Nella variabile double ttt sommo 4 addendi, e poi a temp assegno ttt.

Notate che i 4 addendi nei due estratti di codice sono gli stessi.

Io sarei portato a pensare che value e temp assumano lo stesso valore. E invece NO!!!

Qualcuno sa spiegarmi perchè?

C'è qualcosa che mi sfugge, ma non so cosa. I miei colleghi non riescono ad aiutarmi.

Ovviamente posso darvi tutte le informazioni che volete sul codice in questione, anche se non credo siano molto rilevanti.
Per la cronaca, ptr è un puntatore a un oggetto, il quale oggetto ha, tra gli altri attributi, il double temp3, 4 double (pdd, pdu, pud, puu) e 4 puntatori a oggetti come lui.
Quindi le variabili temp, ttt, value, ptr->pxy e ptr->nxy->temp3 (per x, y = u, d) sono tutte double.

La situazione in sostanza è questa.
Immaginate di avere un albero. I nodi dell'albero sono le strutture di sopra (hanno 4 probabilità e 4 puntatori ad altrettanti nodi, più vari double temporanei).
A questo punto io faccio un ciclo sui nodi dell'albero, e in corrispondenza di ogni nodo mi calcolo il value (vedi codice sopra) e il temp (vedi codice sopra), e li confronto. Ovviamente tra il value e il temp non faccio niente. Sono uno dopo l'altro. Il puntatore ptr viene utilizzato solo in lettura.
Quindi confronto value e temp. E sono diversi.

Qualcuno sa spiegarmi cosa succede?

P.S. ovviamente posso allegare il mio codice se qualcuno non crede a quello che mi succede.

[marco.r]

A colpo d'occhio non vedo problemi.
Che valori hanno le variabili incriminate ? E che valori ottieni ? (cosi' intanto capiamo in quale dei due pezzi di codice le cose vanno male).

[Zero-Giulio]

Non si tratta di valori giusti o sbagliati.
Non saprei neanche quale dei due sia quello giusto, giacchè le differenze sono sulla quindicesima/sedicesima cifra decimale.
Il punto è che non sono uguali, e io non capisco perchè.
Non ci sono conversioni, e il tipo e l'ordine delle operazioni dovrebbe essere lo stesso (per dire, lo so che in aritmetica a precisione finita le operazioni non sono associative, o altri dettagli di questo tipo... Ma nel mio caso dovrebbe essere tutto perfettamente identico).

[Unrue]

double temp = 0.0;
double ttt = (ptr->_pdd * ptr->_ndd->_temp3)
              + (ptr->_pdu * ptr->_ndu->_temp3)
              + (ptr->_pud * ptr->_nud->_temp3)
              + (ptr->_puu * ptr->_nuu->_temp3);
temp = ttt;

In questo caso non conosci l'ordine delle somme, che non necessariamente sono eseguite nell'ordine che hai scritto se attivi le ottimizzazioni ( -O3 ad esempio). Quindi, essendo la somma tra floating point non commutativa, puoi avere valori diversi tra i due casi.

Nel tuo caso potrebbe bastare disabilitare TUTTI i flags di ottimizzazione ( -O0 ad esempio), per verificare in fase di debug se i risultati tra i due casi coincidono.

A+B != B+A con A e B floating point.

[Zero-Giulio]

Ma quindi secondo te è un problema di ordine delle operazioni?
Io ho sempre pensato che l'ordine delle operazioni lo decidessi al moemento di scrivere materialmente il codice.
Cmq flag di ottimizzazione non ne ho.
Sviluppo con un ide, CodeBlocks, ma nella finestra dei flag ho solo -Wall

[Unrue]

Originariamente inviato da Zero-Giulio

Ma quindi secondo te è un problema di ordine delle operazioni?

Supponendo che prima di quelle somme non ci siano altri bug nel codice, si.

Originariamente inviato da Zero-Giulio

Ma quindi secondo te è un problema di ordine delle operazioni?
Io ho sempre pensato che l'ordine delle operazioni lo decidessi al moemento di scrivere materialmente il codice.

No, quello che è garantito è la correttezza delle operazioni anche se vengono scambiate l'ordine di alcune di esse. Ovviamente non ti può modificare la logica del codice.

Originariamente inviato da Zero-Giulio

Cmq flag di ottimizzazione non ne ho.
Sviluppo con un ide, CodeBlocks, ma nella finestra dei flag ho solo -Wall

Controlla che il compilatore non inserisca flags di ottimizzazione e poi metti esplicitamente -O0

[marco.r]

Originariamente inviato da Unrue

No, quello che è garantito è la correttezza delle operazioni anche se vengono scambiate l'ordine di alcune di esse. Ovviamente non ti può modificare la logica del codice.

Uhm, ma l'operatore + in c++ associa a sinistra, quindi l'ordine dovrebbe essere ben definito

[Zero-Giulio]

Domani proverò con l'O0...
Vedremo...

[marco.r]

Originariamente inviato da Unrue

No, quello che è garantito è la correttezza delle operazioni anche se vengono scambiate l'ordine di alcune di esse. Ovviamente non ti può modificare la logica del codice.

La mia ipotesi e' un po' diversa: nel primo caso il valore temporaneo viene sempre riportato ogni volta nella variabile temporanea, per cui in memoria
Nel secondo caso il compilatore puo' aver benissimo deciso di lasciare il risultato intermedio nel registro della FPU che nel caso di CPU x86 ha precisione interna maggiore.
Poi provo a verificare guardando il codice generato nei due casi.

[LMCH]

Originariamente inviato da Zero-Giulio

Qualcuno sa spiegarmi cosa succede?

Sei entrato nel fantabosco del floating point, dove vivono i NaN(senza i), gli zeri mannari ed altre bestie strane ...
e dove gli errori di troncamento ti possono far uscire pazzo.

Più semplicemente con i valori in floating point ad ogni operazione che fai, hai un risultato che in generale è un approssimazione del "vero" risultato.

Ad esempio i double hanno una mantissa normalizzata a 56bit (in realtà 55 con il 56-esimo implicitamente sempre ad 1), se fai una moltiplicazione o una divisione ti servirebbe una mantissa di 112bit (56+56 bit) per avere il risultato esatto-esatto e con somma sottrazione è anche peggio (se uno dei due numeri ha un esponente troppo piccolo rispetto all'altro, è come sommare o sottrarre zero a causa del "troncamento implicito della mantissa").

Per questo motivo sugli x86 se si usa la FPU (invece delle istruzioni SSE) le operazioni vengono eseguiti su registri a precisione estesa (64bit di mantissa e 16bit di esponente, se ricordo bene).
Quando poi scrivi il risultato in memoria su una variabile double per troncamento o arrotondamento perdi di brutto 8 bit di mantissa.

Nel primo spezzone di codice, dopo ogni moltiplicazione (fatta su registri interni) sommavi il risultato e lo memorizzavi su una variabile double, quindi producevi un troncamento/arrotondamento dopo ogni moltiplicazione e somma quando ri-memorizzi il risultato (per un totale di 4 troncamenti che ogni volta tranciano via 8 bit "residui").

Nel secondo spezzone di codice hai fatto tutte le moltiplicazioni e somme in un unica espressione, tenendo i valori intermedi sui registri a precisione estesa
ed alla fine fai un solo troncamento finale per memorizzare il risultato.

Se il compilatore non fa altri giochini di ottimizzazione ecc. ecc. è quindi quello a produrre valori differenti, ma di solito ci si mette pure il compilatore (specialmente se si usano opzioni che accelerano le operazioni a discapito della precisione).

Esistono comunque degli standard su "come gestire correttamente i floating point secondo gli standard IEEE" in modo da ridurre le "differenze di risultato" (non si tratta di avere più precisione, ma di avere risultati abbastanza coerenti.

Con i compilatori Microsoft puoi usare le opzioni /fp ( http://msdn.microsoft.com/en-us/library/e7s85ffb.aspx ).

Con GCC ci sono -ffloat-store -fexcess-precision=... oppure puoi forzare l'uso di solo sse2 in modo da non usare la precisione estesa (-mfpmath=sse -msse2 ecc. ecc.).
ecc. ecc.
Per maggiori dettagli dai un occhiata qui: http://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html
(per trovare quello che ti interessa, vai avanti fino a quando con -ffloat-store inizia la lista delle opzioni che controllano precisione ed ottimizzazioni del floating point).

Normalmente le "differenze di comportamento" sono trascurabili, ma ovviamente in certi casi bisogna fare molta attenzione per evitare che succedano cose strane.
E' per questo che in certi casi per simulazioni scientifiche o per elaborazioni in cui l'accumulazione di errori potrebbe essere pericolosa, si usano librerie che garantiscono maggiore "stabilita dei risultati" tipo MPFR (http://www.mpfr.org/ ) o GMP ( http://gmplib.org/ )

[Unrue]

Originariamente inviato da marco.r

Uhm, ma l'operatore + in c++ associa a sinistra, quindi l'ordine dovrebbe essere ben definito

A meno di trucchi di ottimizzazione .

La cosa migliore in questi casi è analizzare l'assembly generato.

[Zero-Giulio]

Uh, mi sa che sono proprio gli errori che troncamento.
Che noia 'sti double, se i numeri non sono uguali non posso usare WinMerge per fare i test di non regressione (trovo mille mila ddifferenze, anceh se tutte nelle ultime cifre).

Vabbè, oggi non ho avuto molto tempo per testare se veramente è questo il problema (il mio boss mi ammazza se rimango fermo sullo stesso punto per la sedicesima cifra decimale). Magari prox ci ritorno su.

Cmq grazie mille a tutti per le dritte.

[Unrue]

Originariamente inviato da Zero-Giulio

Uh, mi sa che sono proprio gli errori che troncamento.
Che noia 'sti double, se i numeri non sono uguali non posso usare WinMerge per fare i test di non regressione (trovo mille mila ddifferenze, anceh se tutte nelle ultime cifre).

Vabbè, oggi non ho avuto molto tempo per testare se veramente è questo il problema (il mio boss mi ammazza se rimango fermo sullo stesso punto per la sedicesima cifra decimale). Magari prox ci ritorno su.

Cmq grazie mille a tutti per le dritte.

In realtà se il tuo algoritmo è immune ad errori dalla sedicesima cifra decimale puoi anche ignorarli. Ovviamente se invece dà instabilità devi contenerli.

[Zero-Giulio]

Qualche instabilità la da.
L'albero mi serve per la valutazione di alcuni prodotti finanziari, e questa procedura è inserita in un algoritmo di ottimizzazione per la calibrazione dei parametri del modello sui dati di mercato.
Poichè l'algoritmo di ottimizzazione lavora con le derivate numeriche, errori alla sedicesima cifra sui prezzi diventano errori alla ottava cifra nella matrice dello jacobiano e da li in poi...
Nulla di che, alla fine l'output della calibrazione è lo stesso alla settima cifra (e a me tipicamente servono solo le prime 4-5 cifre).

Diciamo che era più che altro curiosità intellettuale.
Io non riuscivo a spiegarmi in cosa potessero differire le due versioni del codice.
Io non sapevo, per dire, dei registri a precisione estesa...

[LMCH]

Originariamente inviato da Zero-Giulio

Nulla di che, alla fine l'output della calibrazione è lo stesso alla settima cifra (e a me tipicamente servono solo le prime 4-5 cifre).

Se vuoi maggior precisione puoi usare le librerie che avevo indicato in precedenza, in particolare la GMP
ed il suo port per MSVC++ MPIR.

[Zero-Giulio]

Ai tempi dell'uni avevo lavorato con la CLN (class library for numbers), e ricordo non ne ero rimasto troppo affascinato.
Credo che la CLN si appoggiasse alla GMP, o sbaglio?
Boh, son passati parecchia anni, non ricordo...