Nei titoli e nei testi troverete qualche rimando cinematografico (ebbene si, sono un cinefilo). Se non vi interessano fate finta di non vederli, già che non sono fondamentali per la comprensione dei post...

Di questo blog ho mandato avanti, fino a Settembre 2018, anche una versione in Spagnolo. Potete trovarla su El arte de la programación en C. Buona lettura.

lunedì 20 dicembre 2021

Double Express
come si comparano i double e i float in C

Ponchia: Ma non ti rovini i denti con tutta quella cioccolata?
Teresa: Ma ho un amico che è dentista.
Ponchia: Cioè, uno che ha un amico dentista allora si deve rovinare i denti. E se io avessi un amico patologo cosa dovrei fare?

La programmazione è un viaggio, un lungo viaggio alla ricerca del codice perfetto... esattamente come il viaggio nostalgico e catartico descritto nel bel Marrakech Express, un film che, metaforicamente, ci mostra la fuga dalla realtà e dal trantran quotidiano di un gruppo di amici. Un film brillantissimo che è ambientato in un momento storico e spensierato che mette un po' di nostalgia visto l'andazzo attuale. Ma bando alla tristezza: siamo qui per trattare di codice di alta qualità! Guardiamo avanti...

...ma non ti rovini il codice con tutti quei double?...

Nell'ultimo articolo avevamo proseguito il nostro viaggio tra le nuove keyword del C99/C11 (abbiamo già visto la prima e la seconda parte, quindi toccherebbe la terza). Ma oggi voglio divagare un po', mi prendo una piccola pausa per cambiare di tema (sento già i sospiri di sollievo: ma siete così sicuri che non sarà una pizza questo nuovo argomento?).

E allora vediamo in po': mi sono accorto, mio malgrado, che nel nostro piccolo-grande mondo dei programmatori ci sono alcuni peccati originali veramente duri a morire. In alcuni articoli passati (ad esempio qui, qui o qui) ho cercato di smitizzarne un po', e oggi cercherò di scrivere la parola fine (si fa per dire...) su un argomento che, navigando in rete, ho notato che continua a generare molti dubbi, e non è che sia un tema poi così nuovo! E quindi, tanto per chiarire:

I DOUBLE (E I FLOAT) NON SI COMPARANO MAI USANDO L'OPERATORE "=="

e vi dirò di più: si dovrebbe anche diffidare degli operatori "<", "<=", ">", ">="  e "!="  ma questo proprio se vogliamo fare gli schizzinosi (e anche questo punto lo vedremo più avanti nell'articolo).

Come detto questa è una vecchia storia, vecchia come i Computer e il Software, ma è ancora sconosciuta a molti. Comunque, per i dubbiosi, facciamo cantare il codice, dove proveremo a comparare due numeri che sono certamente uguali (o no?):

#include <stdio.h>

// main() - funzione main
int main()
{
// a e b dovrebbero valere entrambi 1.3...
double a = (0.3 * 3.0) + 0.4;
double b = 1.3;

// ...e invece...
printf("a = %.32f\n", a);
printf("b = %.32f\n", b);

// ...ecco la prova del disastro!
if (a == b)
printf("i numeri a e b sono uguali!\n");
else
printf("i numeri a e b NON sono uguali...\n");

return 0;
}

ecco, questo semplicissimo codice stampa questo:

a = 1.29999999999999982236431605997495
b = 1.30000000000000004440892098500626
i numeri a e b NON sono uguali...

Sorpresa! a  e b  NON  sono uguali! Beh, in realtà non è una sorpresa, basta indagare un po su come sono rappresentati i numeri double e float in un computer, e come vengono trattati da compilatore e CPU, e quindi ragionare un po' sulle approssimazioni di calcolo... ma non è questa la sede per approfondire questi dettagli, non voglio annoiare nessuno, quindi vi rimando a una delle tante descrizioni che si trovano in rete, ad esempio questa.

Allora: abbiate on non abbiate letto la documentazione consigliata sopra, fidatevi: il problema esiste e, come anticipato più sopra, non è affatto una novità: per esempio il grande Donald Knuth (che cito spesso su queste pagine) ne parlava già nel suo mitico libro The art of computer programming  (Volume II, 1981), dove proponeva alcune interessanti soluzioni del problema. E allora, facendo tesoro delle soluzioni proposte da Knuth, ho scritto questo semplice codice che vi invito a compilare ed eseguire:

#include <stdio.h>
#include <math.h>
#include <stdbool.h>

// prototipi locali
bool isEqual(double a, double b, double myepsilon);

// main() - funzione main
int main()
{
// a e b dovrebbero valere entrambi 1.3...
double a = (0.3 * 3.0) + 0.4;
double b = 1.3;

// ...e invece...
printf("a = %.32f\n", a);
printf("b = %.32f\n", b);

// ...ecco la prova del disastro!
if (a == b)
printf("i numeri a e b sono uguali!\n");
else
printf("i numeri a e b NON sono uguali...\n");

// ma c'è una soluzione!
if (isEqual(a, b, 0.001))
printf("ma con la soluzione i numeri a e b sono uguali!\n");
else
printf("anche con la soluzione i numeri a e b NON sono uguali...\n");

return 0;
}

// isEqual() - funzione di comparazione per variabili double (ispirata dalle
// soluzioni di D.Knuth in "The art of computer programming" (vol.2, 1981))
bool isEqual(
double a, // primo double da comparare
double b, // secondo double da comparare
double myepsilon) // un numero piccolo per riferimento di comparazione
{
if (islessequal(
fabs(a - b),
myepsilon * (isgreater(fabs(a), fabs(b)) ? fabs(a) : fabs(b)))) {

// se sono (approssimativamente) uguali ritorna true
return true;
}

// se sono (approssimativamente) diversi ritorna false
return false;
}

Dopo l'esecuzione vedrete questo:

a = 1.29999999999999982236431605997495
b = 1.30000000000000004440892098500626
i numeri a e b NON sono uguali...
ma con la soluzione i numeri a e b sono uguali!

Ecco, la funzione isEqual() è un buon compromesso tra affidabilità e semplicità per eseguire comparazione tra double (e si può `facilmente scriverne una versione per i float). In realtà, il problema è più complesso di quello che sembra, e le soluzioni possibili sono molteplici, ma quella che vi ho appena proposto è interessante e flessibile (anche se alcuni puristi dicono che non è sufficiente).

Analizziamo il codice: per eseguire il confronto si usa un numero molto piccolo di riferimento che ho chiamato myepsilon per distinguerlo da eventuali epsilon definiti in alcuni linguaggi e compilatori come "numero più piccolo trattabile" (ad esempio nel C++ c'è: std::numeric_limits<T>::epsilon): in realtà non abbiamo bisogno di un epsilon assoluto  ma possiamo usarne uno definito da noi (da qui il nome myepsilon) in base alla precisione che ci serve (e per questo ho definito questa soluzione come flessibile). 

Notare che myepsilon non viene usato come comparatore assoluto ma come moltiplicatore, così la isEqual() adatta il suo comportamento ai numeri in uso (che potrebbero anche essere molto piccoli, più piccoli dello stesso myepsilon).

Avrete visto, poi, che non ho usato direttamente gli operatori ">" e "<=" per effettuare due comparazioni, ma ho invece usato due macro, fornite dalla glibc, che sono state scritte ad-hoc per eseguire queste operazioni con variabili floating-point: e visto che le abbiamo a disposizione usiamole! Sono sicuramente più affidabili degli operatori quando si eseguono operazioni in virgola mobile (anche perché tengono conto di eventuali operandi di tipo NaN). La lista di equivalenza che ci interessa è questa:

operatore macro corrispondente
-----------------------------------------
> int isgreater(x, y);
>= int isgreaterequal(x, y);
< int isless(x, y);
<= int islessequal(x, y);
!= int islessgreater(x, y);

E vi ricordo che, nel C++, le stesse macro hanno comportamento ed uso identici ma ritornano bool invece di int.

Comunque, dopo tutta questa sbrodolata su come fare e non fare le operazioni floating point, vi aggiungo un "consiglio da amico" che resetta il tutto (anzi: se avessi cominciato l'articolo così vi avrei risparmiato tutto il tempo della lettura). Tenetevi forte: la maniera migliore di fare le operazioni in virgola mobile è non farle proprio!

E mi spiego meglio: tornando con i piedi per terra e parlando di casi pratici, ad esempio la scrittura di Software industriale (un classico del C e C++) è sempre una buona idea trattare i dati (chessoio: misure di sensori: tensione, corrente, pressione e chi più ne ha più ne metta) come interi associati a una precisione: quindi il nostro Sottware conterrà un parametro di configurazione che è la precisione, che verrà usata come divisore solo quando si dovranno presentare i dati su un monitor o si dovranno stampare dei report. Ad esempio: supponendo di avere una precisione 1000 (che corrisponde a 3 decimali), una tensione di 1,525V per noi sarà, internamente, un intero che vale 1525, che poi diventerà 1,525 dividendolo per la precisione.

Semplice, no? Eseguendo internamente tutti i calcoli usando solo interi si eliminano automaticamente i problemi descritti in questo articolo. Evidentemente ci sono casi in cui usare direttamente double e float  è indispensabile, ma nel Software industriale il trucco qui sopra è abbastanza usuale, ve l'assicuro.

Beh, per oggi può bastare. Abbiamo asfaltato un altro argomento dubbioso e prometto che ritorneremo (nel prossimo articolo) con l'argomento che avevamo lasciato pendente, quello delle nuove keyword del C (salvo ripensamenti dell'ultima ora, ah, ah, ah).

Ciao, e al prossimo post!