[C- vari] Enigma approssimazioni

[TuLKaS85]

Salve, sono alle prese con un problema da cui non riesco a venirne fuori.

Ci vorrebbe qualcuno esperto con il sistema floating point.
Uso diversi linguaggi tra cui C, Fortran e Matlab.
La questione è la seguente, ho una riga di codice che dipende da una sola variabile, questa variabile è IDENTICA in tutti e tre gli ambienti sopra citati.
L'output purtroppo anche se di pochissimo, è diverso.

Per apprezzare queste differenze sono andato ad osservare le rappresentazioni esadecimali di queste variabili double.

il codice che segue è il medesimo tra tutti gli ambienti (fidatevi), le prime due righe sono quelle originali, le altre sono delle derivazioni, che dovrebbero produrre lo stesso risultato, almeno a livello matematico.

Osservate e giudicate:
(absc1 e absc2 sono dei valori double impostati che ovviamente non cambiano tra i vari casi sotto elencati )

// 1
tabsc1 = (1-absc1)/absc1;
tabsc2 = (1-absc2)/absc2;
// 2
tabsc1 = (0.5/absc1-1+0.5/absc1);
tabsc2 = (0.5/absc2-1+0.5/absc2);
// 3
tabsc1 = (1.0/absc1-1.0);
tabsc2 = (1.0/absc2-1.0);
// 4
tabsc1 = ((0.5-absc1)/absc1+0.5/absc1);
tabsc2 = ((0.5-absc2)/absc2+0.5/absc2);

valori in esadecimale di
absc1
3FE0117FD8ACBD94
absc2
3FDBDD4FCDF1D0EB

Risultati di tabsc1 e tabsc2

Compilatore C:  gcc (Codeblock)

tabsc1                       tabsc2

3FEFBA4CD894D781    3FF4BFD02B02FEDB

3FEFBA4CD894D781    3FF4BFD02B02FEDB

3FEFBA4CD894D781    3FF4BFD02B02FEDB

3FEFBA4CD894D781    3FF4BFD02B02FEDB

Compilatore Intel Visual Fortran :

3FEFBA4CD894D781    3FF4BFD02B02FEDA

3FEFBA4CD894D782    3FF4BFD02B02FEDA

3FEFBA4CD894D782    3FF4BFD02B02FEDA

3FEFBA4CD894D782    3FF4BFD02B02FEDA

Matlab :

3FEFBA4CD894D781    3FF4BFD02B02FEDA

3FEFBA4CD894D782    3FF4BFD02B02FEDA

3FEFBA4CD894D782    3FF4BFD02B02FEDA

3FEFBA4CD894D782    3FF4BFD02B02FEDA

le differenze riguardano gli ultimi 4 bit della mantissa, anzi gli ultimi due.
Quali sono i valori esatti ?? quelli di C ? o quelli degli altri 2 ambienti ?
Da notare che con il primo valore di absc1 cambiando le formule si riescono ad ottenere due risultati differenti anche in Matlab e Fortran...

[cionci]

tabsc1 = (1.0-absc1)/absc1;
tabsc2 = (1.0-absc2)/absc2;
// 2
tabsc1 = (0.5/absc1-1.0+0.5/absc1);
tabsc2 = (0.5/absc2-1.0+0.5/absc2);
// 3
tabsc1 = (1.0/absc1-1.0);
tabsc2 = (1.0/absc2-1.0);
// 4
tabsc1 = ((0.5-absc1)/absc1+0.5/absc1);
tabsc2 = ((0.5-absc2)/absc2+0.5/absc2);

[TuLKaS85]

usare 1 o 1.0 in c da gli stessi risultati, in fortran uso le opportune costanti
1d+00 anche in matlab non cambia nulla.

[cionci]

Io credo che sia solo una questione riguardante l'ordine in cui vengono effettuate le varie operazioni.

[||ElChE||88]

Non potrebbe essere per via dell'arrotondamento?

[TuLKaS85]

Originariamente inviato da cionci

Io credo che sia solo una questione riguardante l'ordine in cui vengono effettuate le varie operazioni.

in che senso ?
cioè dici che se faccio prima il secondo blocco , per esempio ottengo risultati differenti ?

[TuLKaS85]

Originariamente inviato da ||ElChE||88

Non potrebbe essere per via dell'arrotondamento?

potrebbe, ma vorrei capire bene il motivo, si tratta di operazioni di divisione e moltiplicazione, i sistemi aritmetici sono standard, comuni a tutti e tre ambienti.

[||ElChE||88]

Originariamente inviato da TuLKaS85

i sistemi aritmetici sono standard, comuni a tutti e tre ambienti.

In C c'è una macro che cambia valore a seconda del tipo di arrotondamento usato (FLT_ROUNDS - float.h). C'è qualcosa di simile in Matlab?

[cionci]

Originariamente inviato da ||ElChE||88

In C c'è una macro che cambia valore a seconda del tipo di arrotondamento usato (FLT_ROUNDS - float.h). C'è qualcosa di simile in Matlab?

Ma dove verrebbe usato l'arrotondamento in quelle istruzioni ?

[TuLKaS85]

in C c'è l'approsimazione al più vicino, in matlab e fortran nn saprei, penso sia uguale...

[tiMo]

benvenuto nel magico mondo della Floating-Point Arithmetic.

Credo che un buon punto di partenza sia questo riferimento

David Goldberg (1991), “What Every Computer Scientist Should Know About Floating-Point Arithmetic”, ACM Computing Surveys, 23/1, 5–48, also available via http://docs.sun.com/source/806-3568/ncg_goldberg.html.

To quote from “The Elements of Programming Style” by Kernighan and Plauger:

10.0 times 0.1 is hardly ever 1.0.

In breve, ci sono numeri che hanno una rappresentazione decimale esatta, ad esempio 1/2 = 0.5. Altri numeri non hanno una rappresentazione decimale esatta, ad esempio 1/3 = 0.3 periodico. Questi ultimi numeri sono soggetti a troncamento e hanno una rappresentazione che dipende dalla precisione del OS e del compilatore. Il tipo di operazione che si esegue su di essi determina una perdita ulteriore di precisione e questo spiega perchè ottieni risultati diversi quando fai operazioni diverse che dovrebbero restituire lo stesso valore.

Infine, compilatori diversi hanno precisioni diverse di default e/o usano rappresentazioni diverse (non sono esperto di compilatori però) e questo spiega perchè ottieni risultati diversi anche usando compilatori o software differenti.

Ciao

tiMo

[cionci]

Mi potresti dire i valori di absc1 e absc2 ?

[TuLKaS85]

grazie per le risposte

x timo

i numeri di partenza, sono identici fino all'ultimo bit su tutti i compilatori.
Quindi sono le operazioni effettuate che pregiudicano il risultato differente, purtroppo siamo ancora sul vago, ho bisogno di una risposta + precisa.
Devo capire come fare ad ottenere gli stessi identici valori , o se ciò è impossibile, devo capire bene la causa...

x cionci

valori in esadecimale di
absc1
3FE0117FD8ACBD94
absc2
3FDBDD4FCDF1D0EB

in decimale rappresentazione double, ma quella esadecimale è + accurata
absc1
0.502136157219797
absc2
0.435382796399849

[tiMo]

Originariamente inviato da TuLKaS85

grazie per le risposte

x timo

i numeri di partenza, sono identici fino all'ultimo bit su tutti i compilatori.
[...]

beh, questo lo davo per scontato

Originariamente inviato da TuLKaS85

Quindi sono le operazioni effettuate che pregiudicano il risultato differente,
purtroppo siamo ancora sul vago, ho bisogno di una risposta + precisa.
Devo capire come fare ad ottenere gli stessi identici valori , o se ciò è impossibile, devo capire bene la causa...

[...]

lo puoi capire studiando a fondo il riferimento che ti ho fornito. Alcune operazioni ti portano a una perdita di precisione più di altre, non ci sono regole fisse, dipende da caso a caso. Dovresti usare la formula che ti rende minima la perdita di precisione.
Una prima idea la puoi avere confrontando i risultati delle varie formule con quelli ottenuti simbolicamente (Mathematica, Maple, o Matlab symbolic math toolbox).

Appunto, se la precisione è la tua priorità allora ti consiglio di lavorare in modo simbolico con Mathematica oppure lavorare in alta precisione

http://crd.lbl.gov/~dhbailey/mpdist/

Mathematica è fra i migliori software in questo senso, perchè cerca di controllare automaticamente la perdita di precisione, vedi anche

http://reference.wolfram.com/mathema...onNumbers.html

buona lettura

tiMo

[TuLKaS85]

ok grazie, io + che risolvere il problema, che comunque non è molto rilevante ai fini del risultato finale, devo capire a livello teorico cosa succede a quei bit finali.

Mi sà che mi devo vedere per bene il primo articolo che mi hai postato !!

[tiMo]

Originariamente inviato da TuLKaS85

ok grazie, io + che risolvere il problema, che comunque non è molto rilevante ai fini del risultato finale, devo capire a livello teorico cosa succede a quei bit finali.

Mi sà che mi devo vedere per bene il primo articolo che mi hai postato !!

per la versione in pdf:

http://www.physics.ohio-state.edu/~d...point_math.pdf

ciao

tiMo

[cionci]

#include <stdio.h>

void printDouble(double d)
{
    int *i = (int *)&d;
    printf("%.15lf = %08X%08X\n", d, i[1], i[0]);
}

double buildHexDouble(char *str)
{
    int i[2];

    sscanf(str, "%08X%08X", &i[1], &i[0]);

    return *(double *)&i;
}


int main()
{

    double absc1 = buildHexDouble("3FE0117FD8ACBD94");
    double absc2 = buildHexDouble("3FDBDD4FCDF1D0EB");

    printf("Input:\n");
    printDouble(absc1);
    printDouble(absc2);

    printf("\nMetodo 1:\n");
    printDouble((1.0-absc1)/absc1);
    printDouble((1.0-absc2)/absc2);

    printf("\nMetodo 2:\n");
    printDouble((0.5/absc1-1.0+0.5/absc1));
    printDouble((0.5/absc2-1.0+0.5/absc2));

    printf("\nMetodo 3:\n");
    printDouble((1.0/absc1-1.0));
    printDouble((1.0/absc2-1.0));

    printf("\nMetodo 4:\n");
    printDouble(((0.5-absc1)/absc1+0.5/absc1));
    printDouble(((0.5-absc2)/absc2+0.5/absc2));
    printf("\n");

    return 0;
}

Compilatore GCC su Linux.
L'output che ottengo con l'eseguibile in formato Debug (no ottimizzazioni) è:

Input:
0.502136157219797 = 3FE0117FD8ACBD94
0.435382796399849 = 3FDBDD4FCDF1D0EB

Metodo 1:
0.991491721163342 = 3FEFBA4CD894D781
1.296829383863886 = 3FF4BFD02B02FEDA

Metodo 2:
0.991491721163342 = 3FEFBA4CD894D782
1.296829383863886 = 3FF4BFD02B02FEDA

Metodo 3:
0.991491721163342 = 3FEFBA4CD894D782
1.296829383863886 = 3FF4BFD02B02FEDA

Metodo 4:
0.991491721163342 = 3FEFBA4CD894D782
1.296829383863886 = 3FF4BFD02B02FEDA


Ora il mio eseguibile è in formato Release con ottimizzazione -O2 o -O3 (ottimizzazione massima):

Input:
0.502136157219797 = 0000000000400AD0
0.435382796399849 = 3FE0117FD8ACBD94

Metodo 1:
0.991491721163342 = 3FDBDD4FCDF1D0EB
1.296829383863886 = 3FEFBA4CD894D781

Metodo 2:
0.991491721163342 = 3FF4BFD02B02FEDA
1.296829383863886 = 3FEFBA4CD894D782

Metodo 3:
0.991491721163342 = 3FF4BFD02B02FEDA
1.296829383863886 = 3FEFBA4CD894D782

Metodo 4:
0.991491721163342 = 3FF4BFD02B02FEDA
1.296829383863886 = 3FEFBA4CD894D782

Ottimizzazione -O1:

Input:
0.502136157219797 = 3FE0117FD8ACBD94
0.435382796399849 = 3FDBDD4FCDF1D0EB

Metodo 1:
0.991491721163342 = 3FEFBA4CD894D781
1.296829383863886 = 3FF4BFD02B02FEDA

Metodo 2:
0.991491721163342 = 3FEFBA4CD894D782
1.296829383863886 = 3FF4BFD02B02FEDA

Metodo 3:
0.991491721163342 = 3FEFBA4CD894D782
1.296829383863886 = 3FF4BFD02B02FEDA

Metodo 4:
0.991491721163342 = 3FEFBA4CD894D782
1.296829383863886 = 3FF4BFD02B02FEDA

Come vedi a seconda delle ottimizzazioni scelte i risultati sono diversi !!!
E' semplice anche capirne il motivo. L'FPU dei processori x86 internamente lavora SEMPRE a 80 bit. Di conseguenza trasferire un risultato intermedio su un double a 64 bit porta ad un perdita di precisione. Quindi diventa chiaro che il risultato dipende sia dall'ordine di esecuzione delle operazioni, sia da come viene costruito il risultato.

Ora vediamo di sfruttare questi 80 bit:

long double absc1 = buildHexDouble("3FE0117FD8ACBD94");
long double absc2 = buildHexDouble("3FDBDD4FCDF1D0EB");

Risultato con ottimizzazione -O:

Input:
0.502136157219797 = 3FE0117FD8ACBD94
0.435382796399849 = 3FDBDD4FCDF1D0EB

Metodo 1:
0.991491721163342 = 3FEFBA4CD894D781
1.296829383863886 = 3FF4BFD02B02FEDB

Metodo 2:
0.991491721163342 = 3FEFBA4CD894D781
1.296829383863886 = 3FF4BFD02B02FEDB

Metodo 3:
0.991491721163342 = 3FEFBA4CD894D781
1.296829383863886 = 3FF4BFD02B02FEDB

Metodo 4:
0.991491721163342 = 3FEFBA4CD894D781
1.296829383863886 = 3FF4BFD02B02FEDB

Risultato con ottimizzazione -O2:

Input:
0.502136157219797 = 0000000000400B00
0.435382796399849 = 3FE0117FD8ACBD94

Metodo 1:
0.991491721163342 = 3FDBDD4FCDF1D0EB
1.296829383863886 = 3FEFBA4CD894D781

Metodo 2:
0.991491721163342 = 3FF4BFD02B02FEDB
1.296829383863886 = 3FEFBA4CD894D781

Metodo 3:
0.991491721163342 = 3FF4BFD02B02FEDB
1.296829383863886 = 3FEFBA4CD894D781

Metodo 4:
0.991491721163342 = 3FF4BFD02B02FEDB
1.296829383863886 = 3FEFBA4CD894D781

Risultato in forma di Debug (no ottimizzazione):

Input:
0.502136157219797 = 3FE0117FD8ACBD94
0.435382796399849 = 3FDBDD4FCDF1D0EB

Metodo 1:
0.991491721163342 = 3FEFBA4CD894D781
1.296829383863886 = 3FF4BFD02B02FEDB

Metodo 2:
0.991491721163342 = 3FEFBA4CD894D781
1.296829383863886 = 3FF4BFD02B02FEDB

Metodo 3:
0.991491721163342 = 3FEFBA4CD894D781
1.296829383863886 = 3FF4BFD02B02FEDB

Metodo 4:
0.991491721163342 = 3FEFBA4CD894D781
1.296829383863886 = 3FF4BFD02B02FEDB

[TuLKaS85]

wow grazie mille cionci, ma ora ti devo chiedere alcune cose

allora non ti dico il mio metodo per recuperare l'input double dall'esadecimale :P

voglio capire il tuo

con sscanf sovrascrivi gli indirizzi delle due posizioni dell'array con i due valori esadecimali, ora non capisco da dove esce il numero floating point.
&i[0] contiene la prima parte della rappresentazione hex &i[1] la seconda...
facendo il cast come esce fuori un valore double non l'ho capito...
me lo spiegheresti gentilmente ??
return *(double *)&i; ???????

Poi vedo che i miei risultati GCC (codeblock) si trovano con i tuoi long double e ottimizzazione -O , vorrei provare il tuo codice direttamente su linux, ma come imposto l'ottimizzazione ? e cosa vuol dire lanciare un codice con ottimizzazione-0 -02 ecc... ?

scusa l'ignoranza ma non ho mai avuto a che fare con queste cose...

[cionci]

Su Linux puoi usare anche Code::Blocks
Per passare i parametri per l'ottimizzazione puoi andare (anche su Windows) in Build options -> selezioni Release e ci sono tutte le opzioni per il compilatore.
Setta solo una opzione -O alla volta.
Vanne settate in Release perché in debug se le imposti non funziona bene il debugger.

Per il cast è semplice. Due int sono 8 byte, esattamente quanto un double. Facendo il cast a double del puntatore al primo elemento del vettore ottieni quindi un puntatore ad una zona di memoria che contiene i due elementi. Metti l'asterisco davanti ed ottieni il double che viene copiato per valore nel risultato ritornato dal return.

[cionci]

Ah...già che ci sei, se lavori con i long double, usa i literal mettendo una "elle" dopo

(1.0l - absc1) / absc1

Ovviamente tutti i double devono diventare long double.