print(0.1 + 0.2 == 0.3)False11 L’aritmetica in virgola mobile
Cominciamo con una piccola provocazione, nella forma della singola riga di codice che segue:
Sorprese? Non vi sono dubbi che sul libro di analisi matematica questa eguaglianza \[ \frac{1}{10} + \frac{2}{10} = \frac{3}{10} \] sarebbe vera. Perché, invece, Python ci dice che è falsa? Si tratta di un bug di uno dei più linguaggi di programmazione più diffusi al mondo, oppure la cosa nasconde qualcosa di profondo?
L’aritmetica in virgola mobile è sempre una sorgente di sorprese per chi si avvicina per la prima volta alla programmazione. In questo capitolo cercheremo di mettere a fuoco le questioni principali. Tra le varie risorse sull’argomento presenti su web, segnaliamo questo tutorial, parte della documentazione ufficiale di Python.
11.1 Virgola mobile e rappresentazione in base 2
Nel capitolo Capitolo 7 abbiamo imparato come si rappresentano in base 2 i numeri interi. Se ricordate, avevamo lasciato deliberatamente da parte i numeri in virgola mobile, e adesso è arrivato il momento di riprendere e sviscerare la questione.
Partiamo da un esempio semplice. Che cosa intendiamo quando scriviamo il numero \(2.25\) (nell’ordinaria rappresentazione in base 10)? Evidentemente \[\begin{align*} 2.25 = 2 + 2 \times 10^{-1} + 5 \times 10^{-2}, \end{align*}\] il che vuol dire che, se non siamo interessate solamente ai numeri interi, dobbiamo generalizzare la nostra definizione nell’equazione Equazione 7.1 come \[\begin{align} (c_m\ldots c_2c_1c_0.c_{-1}c_{-2}\ldots)_b = c_m \times b^m + \cdots + c_2 \times b^2 + c_1 \times b + c_0 + c_{-1} b^{-1} + c_{-2} b^{-2} \ldots \end{align}\] In questo schema possiamo scrivere il nostro \(2.25\) nella rappresentazione in base 2 come \[\begin{align*} 0b10.01 = 1 \times 2 + 0 + 0 \times 2^{-1} + 1 \times 2^{-2} = 2 + \frac{1}{4} = 2.25. \end{align*}\]
Fino a qui niente di strano. Ma cosa succede se proviamo a fare la stessa cosa con il numero (apparentemente innocuo) \(0.1\)? Per rispondere possiamo provare ad approssimare (da sopra e da sotto) con un numero fissato di cifre binarie dopo il punto fino a che non otteniamo il risultato esatto:
| 0 | 0/2 | 0b0.0 | 0.1 | 0b0.1 | 1/2 | 0.5 |
| 0 | 0/4 | 0b0.00 | 0.1 | 0b0.01 | 1/4 | 0.25 |
| 0 | 0/8 | 0b0.000 | 0.1 | 0b0.001 | 1/8 | 0.125 |
| 0.0625 | 1/16 | 0b0.0001 | 0.1 | 0b0.0010 | 2/16 | 0.125 |
| 0.09375 | 3/32 | 0b0.00011 | 0.1 | 0b0.00100 | 4/32 | 0.125 |
| 0.09375 | 6/64 | 0b0.000110 | 0.1 | 0b0.000111 | 7/64 | 0.109375 |
| 0.09375 | 12/128 | 0b0.0001100 | 0.1 | 0b0.0001101 | 13/128 | 0.101562 |
| 0.0976562 | 25/256 | 0b0.00011001 | 0.1 | 0b0.00011010 | 26/256 | 0.101562 |
| 0.0996094 | 51/512 | 0b0.000110011 | 0.1 | 0b0.000110100 | 52/512 | 0.101562 |
Possiamo andare avanti quanto vogliamo, ma non è difficile convincersi che il numero decimale \(0.1\) non ha uno sviluppo finito nella sua rappresentazione binaria. Possiamo approssimare il valore con un numero arbitrario \(n\) di cifre (binarie) \[ \frac{\left\lfloor \frac{2^n}{10} \right\rfloor}{2^n} < 0.1 < \frac{\left\lfloor \frac{2^n}{10} \right\rfloor + 1}{2^n}, \] ma lo sviluppo binario di \(0.1\) è periodico: \[ 0.1 = 0b0.00011001100110011001100110011001100110011001100110011 \ldots = 0b0.0\overline{0}\overline{0}\overline{1}\overline{1} \]
Se ci pensate un attimo, la cosa non è sorprendente, perché anche nell’ordinaria rappresentazione decimale, molti numeri razionali non hanno un’espansione finita—ad esempio \[ \frac{1}{3} = 0.\overline{3} \] (E, in effetti, in un senso che non precisiamo, quai tutti i numeri razionali, in una qualsiasi base intera, non hanno una rappresentazione finita.)
Nota
La condizione necessaria e sufficiente perché un numero razionale, espresso nei suoi minimi termini dalla frazione \(p/q\), abbia uno sviluppo finito in base \(b\) è tutti i fattori primi del denominatore \(q\) siano anche fattori primi di \(b\). Così in base \(10\) un numero razionale ha uno sviluppo finito se e solo se il suo denominatore (dopo la riduzione ai minimi termini) è un prodotto di potenze arbitrarie di \(2\) e \(5\). In base \(2\) (e \(16\)), la condizione è più restrittiva, perché l’unico fattore primo della base è, appunto, il \(2\).
Notiamo, per completezza, che questo è uno dei motivi per cui, a partire dai Babilonesi e fino all’avvento dei calcolatori, la base \(60\) (che ha \(2\), \(5\) e \(3\) come fattori primi) è stata molto utilizzata per i calcoli—molti numeri razionali che non hanno rappresentazione finita in base \(10\) ce l’hanno in base \(60\).
11.2 Di nuovo sulla nostra provocazione iniziale
Abbiamo più o meno tutti i pezzi necessari per sbrogliare l’arcano della provocazione con cui abbiamo aperto questo capitolo. E la chiave per la soluzione è fatta da due fattori distinti, ugualmente importanti:
- la maggior parte dei numeri in virgola mobile, come abbiamo detto, non ha una rappresentazione finita in nessuna base intera per cui (indipendentemente dalla base utilizzata internamente dal calcolatore) è impossibile darne una rappresentazione esatta a meno di avere memoria infinita—nello stesso modo in cui non possiamo scrivere tutta l’espansione decimale della frazione \(1/3\);
- molti dei numeri in virgola mobile che possiamo rappresentare esattamente in base \(10\) non hanno uno sviluppo finito in base \(2\).
Prese singolarmente, queste due cose sono abbastanza ovvie. Ma quando le mettiamo insieme, e realizziamo che i numeri nei nostri programmi sono tipicamente rappresentati in base \(10\), mentre il computer li immagazzina internamente rappresentandoli in base \(2\), i risultati sono potenzialmente disastrosi.
In altre parole: quando scriviamo \(0.1\), \(0.2\) e \(0.3\), implicitamente stiamo assumendo che questi numeri si possano rappresentare esattamente, perché il loro sviluppo decimale è finito. Ma per il calcolatore, che lavora in base \(2\), nessuna tra le frazioni \[ \frac{1}{10},~\frac{1}{5}~\text{e}~\frac{3}{10} \] può essere rappresentata esattamente, perché entrambe hanno il fattore primo \(5\) al denominatore. Come abbiamo visto, essenzialmente il computer alloca in memoria solo numeri razionali in cui il denominatore è una potenza di \(2\), limitata in grandezza dal numero di bit che abbiamo a disposizione.
Possiamo rendere la cosa più quantitativa e capire cosa succede fino all’ultimo bit usando la funzione as_integer_ratio(), che ci dice esattamente qual è la frazione allocata dal computer
for x in (0.1, 0.2, 0.3):
numerator, denominator = x.as_integer_ratio()
print(f"{x} -> {numerator} / {denominator} = {x:.20f}")0.1 -> 3602879701896397 / 36028797018963968 = 0.10000000000000000555
0.2 -> 3602879701896397 / 18014398509481984 = 0.20000000000000001110
0.3 -> 5404319552844595 / 18014398509481984 = 0.29999999999999998890Il che vuol dire che la mappatura tra quello che noi scriviamo nel codice sorgente e quello che il computer alloca in memoria è \[ \begin{aligned} 0.1 \rightarrow \frac{3602879701896397}{36028797018963968} = \frac{3602879701896397}{2^{55}} \approx 0.10000000000000000555 \\[5pt] 0.2 \rightarrow \frac{3602879701896397}{18014398509481984} = \frac{3602879701896397}{2^{54}} \approx 0.20000000000000001110 \\[5pt] 0.3 \rightarrow \frac{5404319552844595}{18014398509481984} = \frac{5404319552844595}{2^{54}} \approx 0.29999999999999998890 \end{aligned} \]
Ah, allora abbiamo capito! Quando facciamo la somma dei primi due numeri otteniamo \[ 0.1 + 0.2 \rightarrow \frac{3602879701896397}{2^{55}} + \frac{3602879701896397}{2^{54}} = \frac{10808639105689191}{2^{55}} \] che, chiaramente, non è proprio uguale alla rappresentazione che conosciamo già di \(0.3\) \[ 0.3 \rightarrow \frac{5404319552844595}{2^{54}} = \frac{10808639105689190}{2^{55}} \] (Notate la differenza sull’ultima cifra.)
Giusto? Non proprio. Quando abbiamo fatto la somma ci siamo dimenticate che la precisione con cui possiamo scrivere la mantissa è limitata a \(53\) bit, per cui quando l’esponente è negativo, il numeratore della frazione ridotta che rappresenta il nostro numero in virgola mobile è limitata a \(2^{53} - 1 = 9007199254740991\) che, chiaramente, noi abbiamo superato. Ma ci siamo quasi. Cosa possiamo fare per riportare il numeratore entro l’intervallo permesso? Facile: dividiamo per \(2\) numeratore e denominatore—se non ché il numeratore è dispari per cui, nel processo di divisione, siamo costretti ad arrotondare all’intero (pari) più vicino \[ \frac{10808639105689191}{2^{55}} \rightarrow \frac{5404319552844596}{2^{54}} = \frac{1351079888211149}{2^{52}} \] (Nell’ultimo passaggio abbiamo ridotto la frazione ai minimi termini, perché il numeratore era multiplo di \(4\)). Riprova per le più scettiche
x = (0.1 + 0.2)
numerator, denominator = x.as_integer_ratio()
print(f"(0.1 + 0.2) -> {numerator} / {denominator} = {x:.20f}")(0.1 + 0.2) -> 1351079888211149 / 4503599627370496 = 0.30000000000000004441ovverosia \[ 0.1 + 0.2 \rightarrow \frac{1351079888211149}{4503599627370496} = \frac{1351079888211149}{2^{52}} \approx 0.30000000000000004441 \]
Adesso possiamo chiudere il cerchio. Il confronto con cui abbiamo aperto questo capitolo era falso semplicemente perché \[ (0.1 + 0.2) - 0.3 \rightarrow \frac{1351079888211149}{2^{52}} - \frac{5404319552844595}{2^{54}} = \frac{1}{2^{54}} \] che potete verificare da sole nell’interprete Python
x = (0.1 + 0.2) - 0.3
numerator, denominator = x.as_integer_ratio()
print(f"(0.1 + 0.2) - 0.3 -> {numerator} / {denominator} = {x:.20f}")(0.1 + 0.2) - 0.3 -> 1 / 18014398509481984 = 0.00000000000000005551Wow—questa è stata una sezione piuttosto densa, e vale la pena fermarsi un attimo per distinguere le cose importanti da quelle secondarie. Nella vita reale devo essere capace di replicare su carta e penna, bit per bit, quello che succede in ogni operazione tra numeri in virgola mobile su un calcolatore digitale? La risposta è: non fa male, ma non è strettamente necessario. (Anche se è un argomento interessante che, sotto opportune condizioni, può risolvere una conversazione a cena.)
Quello che dovete aver chiaro, invece, è che anche numeri apparentemente innocui come \(0.1\) sono rappresentati, all’interno del calcolatore, in modo intrinsecamente inesatto. E queste inesattezze fanno cose interessanti nelle normali operazioni aritmetiche, per cui è necessario fare un po’ di attenzione quando interpretiamo il risultato.
11.3 Virtù e limiti dell’aritmetica in virgola mobile
A questo punto il fatto che l’aritmetica in virgola mobile su un calcolatore digitale sia, per sua natura, inesatta, dovrebbe essere impresso in modo indelebile nella vostra testa. In questa sezione cercheremo di esplorare alcune delle implicazioni di questo fatto fondamentale—e nel farlo ci concentreremo sullo standard IEEE-754 a doppia precisione, che è quello che Python utilizza di default.
11.3.1 Range e risoluzione
Il range disponibile per i numeri in virgola è determinato principalmente dal numero di bit allocati per l’esponente. Nello standard a doppia precisione
- la mantissa varia tra \(1\) e \(2\) (escluso);
- l’esponente varia tra \(-1022\) e \(1023\)
e i numeri normali si estendono da \[ 1 \times 2^{-1022} \approx 2.2 \times 10^{-308} \quad\text{fino a }\quad (2 - 2^{-52}) \times 2^{1023} \approx 2 \times 2^{1023} \approx 1.8 \times 10^{308} \]
Nota
Per completezza: quando l’esponente assume il valore minimo possibile, che corrisponderebbe a \(-1023\), lo standard prevede di passare alla rappresentazione sub-normale, in cui l’esponente vale ancora \(-1022\), ma il bit più significativo della mantissa si assume essere \(0\) anziché \(1\). Si tratta di un regime completamente diverso, in cui la precisione relative cessa di essere costante, ma il range si estende verso lo zero fino a \[ 2^{-52} \times 2^{-1022} = 2^{-1074} \approx 4.9 \times 10^{-324} \]
Guardiamo le limitazioni sul range in opera:
print(1.e-324)
print(5.e-324)
print(1.5e308)
print(2.e308)0.0
5e-324
1.5e+308
inf(Occhio alle operazioni matematiche—non servono numeri stratosferici per triggerare overflow: \(171!\) è troppo grande per essere convertito in virgola mobile, così come \(144^{144}\), oppure \(e^{710}\).)
Come abbiamo capito, l’aritmetica in virgola mobile è progettata per operare a precisione relativa costante. Quest’ultima è determinata essenzialmente dai bit allocati per la mantissa: con \(53\) bit effettivi (i \(52\) fisici e quello più significativo, che si assume essere \(1\)) abbiamo \[ 2^{53} \approx 9 \times 10^{15} \] valori possibili nell’intervallo \([1, 2[\), per cui la rappresentazione in doppia precisione permette di avere circa \(16\) cifre significative—\(2^{-53}\) corrisponde all’errore (relativo) di arrotondamento massimo, che a volte viene anche chiamato half spacing.
Nota
Fermatevi un attimo a riflettere su quest’ultima cosa. Siccome nella propagazione degli errori su moltiplicazioni e divisioni si sommano le incertezze relative, l’aritmetica in virgola mobile, by design, permette di moltiplicare o dividere tra loro senza problemi numeri che differiscono di molti ordini di grandezza. Ma la storia è completamente diversa per la somma e la differenza, come vedremo in dettaglio tra un attimo!
Il fatto che la precisione relativa sia costante, implica che la spaziatura fisica tra valori consecutivi per i numeri in virgola mobile non è costante, ma dipende dal valore del numero che vogliamo rappresentare—più specificamente è approssimativamente proporzionale a quest’ultimo. Così, se per i numeri tra \(1\) e \(2\) la spaziatura tra due numeri in virgola mobile consecutivi è \(2^{-52}\), tra \(2\) e \(4\) diventa \(2^{-51}\) e così via: raddoppia ogni volta che si attraversa una potenza di \(2\). Il modulo math della libreria standard di Python fornisce la funzione ulp (unit in the last place) che, per un numero in virgola mobile dato \(x\) restituisce la distanza dal numero immediatamente più alto che si può rappresentare
import math
for x in (1., 1.e10, 1.e20, 1.e50, 1.e100):
print(math.ulp(x))2.220446049250313e-16
1.9073486328125e-06
16384.0
2.076918743413931e+34
1.942668892225729e+84Quando si arriva a valori dell’ordine di \(2^{52} \approx 4.5 \times 10^{15}\) la spaziatura diventa più grande di \(1\) e possono succedere cose curiose. Se volessimo trattate il numero di Avogadro come un vero intero, ad esempio, incapperemmo nella situazione paradossale in cui aggiungere una molecola non ha effetto:
import math
N = 6.02e23
print(math.ulp(N))
print(N == N + 1)67108864.0
True(Intorno al numero di Avogadro la spaziatura della rappresentazione in virgola mobile, ovvero il nostro quanto di molecole, è di circa \(70,000,000\).)
A questo punto possiamo chiederci: come mai chi ha progettato lo standard IEEE 754 ha deciso di allocare i 64 disponibili proprio in questo modo—1 per il segno, 52 per la mantissa e 11 per l’esponente? Abbiamo visto che il numero di bit \(n_e\) che riserviamo per l’esponente determina il range dinamico del nostro standard, ovvero sia il numero di ordini di grandezza (o decadi) tra il minimo ed il massimo numero rappresentabile \[ \text{dex}(x_{\text{min}}, x_{\text{max}}) = \log_{10} \left( \frac{x_{\text{max}}}{x_{\text{min}}} \right) \approx \log_{10} \left( 2^{2^{n_e}}\right) = 2^{n_e} \log_{10} 2 \] Il numero di bit \(n_m\) che riserviamo per la mantissa, d’altra parte, determina la precisione relativa con cui possiamo rappresentare i numeri, o, equivalentemente, il numero di cifre significative nell’ordinaria rappresentazione decimale \[ \log_{10} \left( 2^{n_m + 1} \right) = (n_m + 1) \log_{10} 2 \]
La tabella che segue mostra il range dinamico ed il numero di cifre significative che otterremmo con diverse scelte (per 11 bit di esponente e 52 + 1 di mantissa torniamo al nostro range dinamico di circa 616 decadi, con una precisione di circa 16 cifre significative).
| Bit di esponente | Bit di mantissa | Range dinamico | Cifre significative |
|---|---|---|---|
| 5 | 58 + 1 | 9.6 | 17.8 |
| 6 | 57 + 1 | 19.3 | 17.5 |
| 7 | 56 + 1 | 38.5 | 17.2 |
| 8 | 55 + 1 | 77.1 | 16.9 |
| 9 | 54 + 1 | 154.1 | 16.6 |
| 10 | 53 + 1 | 308.3 | 16.3 |
| 11 | 52 + 1 | 616.5 | 16.0 |
| 12 | 51 + 1 | 1233.0 | 15.7 |
| 13 | 50 + 1 | 2466.0 | 15.4 |
| 14 | 49 + 1 | 4932.1 | 15.1 |
| 15 | 48 + 1 | 9864.2 | 14.8 |
Cosa possiamo dire, come fisiche? Beh, se pensiamo alle scale di lunghezza possiamo dire che ci sono circa 42 ordini di grandezza tra il raggio del protone ed il raggio dell’universo osservabile—61 se partiamo dalla scala di Planck. Ci sono circa 90 ordini di grandezza tra la massa del neutrino più pesante (che ancora non conosciamo) e la massa che stimiamo comporre l’universo. In altre parole: non è facile trovare un range dinamico rilevante in Fisica che ecceda di molto i 100 ordini di grandezza.
E per la precisione? La misura più recente del momento magnetico dell’elettrone (Fan et al. 2023), con una precisione di 0.1 ppb (part per billion), ovvero con 13 cifre significative, detiene il record della proprietà meglio misurata di una particelle. Ovviamente le nostre amiche metrologhe non smettono di fare progressi, ed il record più recente di precisione sulla misura di una frequenza (Marshall et al. 2025) sfiora le 19 cifre significative.
Ma nel complesso possiamo dire che i 64 bit dello standard IEEE sono ben allocati: ci permettono di rappresentare comodamente tutti gl intervalli dinamici che potrebbero essere rilevanti per la Fisica, con una precisione che supera abbondantemente quella delle misure più all’avanguardia, con l’eccezione di alcune applicazioni molto specifiche.
11.3.2 Un esempio di accortezza numerica
Quindi siamo salve? Non proprio, perché a queste condizioni la tragedia è potenzialmente dietro l’angolo attraverso underflow/overflow nei calcoli intermedi e/o cancellazioni catastrofiche non volute. Prendiamo ad esempio la nostra amica distribuzione di Poisson \[ P(k; \mu) = \frac{\mu^k e^{-\mu}}{k!} \] Da un punto di vista numerico, l’espressione è un coacerbo di tutte le funzioni matematiche che sono più prone all’overflow, tanto è vero che non bisogna arrivare a numeri stratosferici per far fallire l’implementazione più naive che una potrebbe immaginare—provare per credere!
import math
def poisson(k, mu):
return mu**k * math.exp(-mu) / math.factorial(k)
mu = 10.
for k in (100, 150, 200):
print(f"P({k}, {mu}) = {poisson(k, mu)}")P(100, 10.0) = 4.864649182067611e-63
P(150, 10.0) = 7.94624168593895e-118--------------------------------------------------------------------------- OverflowError Traceback (most recent call last) Cell In[10], line 8 4 return mu**k * math.exp(-mu) / math.factorial(k) 5 6 mu = 10. 7 for k in (100, 150, 200): ----> 8 print(f"P({k}, {mu}) = {poisson(k, mu)}") Cell In[10], line 4, in poisson(k, mu) 3 def poisson(k, mu): ----> 4 return mu**k * math.exp(-mu) / math.factorial(k) OverflowError: int too large to convert to float
(Avete capito bene. Se non ci mettiamo un pizzico di intelligenza, calcolare la probabilità che una variabile Poissoniana con media \(10\) valga \(200\) diventa un problema intrattabile.)
L’espediente standard, qui, è passare ai logaritmi, e passare attraverso la funzione \(\Gamma\) per calcolare il logaritmo del fattoriale \[
\log P(k; \mu) = k \log \mu - \mu - \log k! = k \log \mu - \mu - \log \Gamma(k + 1)
\] Curiosamente, il modulo math della libreria standard di Python fornisce proprio una funzione che restituisce il logaritmo della funzione \(\Gamma\), per cui una possibile implementazione di questo approccio potrebbe essere
import math
def poisson(k, mu):
logp = k * math.log(mu) - mu - math.lgamma(k + 1)
return math.exp(logp)
mu = 10.
for k in (100, 150, 200, 400):
print(f"P({k}, {mu}) = {poisson(k, mu)}")P(100, 10.0) = 4.864649182067763e-63
P(150, 10.0) = 7.946241685938612e-118
P(200, 10.0) = 5.7566064628489175e-180
P(400, 10.0) = 0.0Come vedete, adesso il caso con \(k = 200\) è trattato correttamente—il problema non era il valore finale, ma i numeri nei calcoli intermedi, che passando ai logaritmi diventano molto più piccoli. Ovviamente ad un certo punto ci scontriamo con il fatto che il risultato è troppo piccolo per essere rappresentato, ed in quel caso otteniamo semplicemente \(0\) come risposta, senza alcun overflow.
Nota
Resistete alla tentazione di re-implementare la Poisson nella vita reale. Il modulo scipy.stats offre, tra le altre cose, un’implementazione della distribuzione di Poisson robusta e ben testata (anche se l’implementazione non è molto diversa da quella che abbiamo appena visto) e non c’è ragione di re-inventare la ruota!
from scipy.stats import poisson
mu = 10.
for k in (100, 150, 200, 400):
print(f"P({k}, {mu}) = {poisson(mu).pmf(k)}")P(100, 10.0) = 4.864649182067763e-63
P(150, 10.0) = 7.946241685938612e-118
P(200, 10.0) = 5.7566064628489175e-180
P(400, 10.0) = 0.011.3.3 Proprietà dell’aritmetica elementare
Noi tutte diamo per scontato le proprietà dell’aritmetica cha abbiamo studiato alle scuole elementari: la proprietà commutativa della somma e del prodotto, la proprietà associativa e la proprietà distributiva della moltplicazione rispetto alla somma. Ebbene, il carattere intrinsecamente inesatto dell’artitmetica in virgola mobile fa sì che queste proprietà, in generale, non possano essere preservate tutte insieme in un calcolatore digitale. (Giusto per chiarezza: questo è vero solo in virgola mobile—l’aritmetica intera, da questo punto di vista, si comporta come ci aspetteremmo.)
a = 0.2
b = 0.1
c = 0.3
print((a + b) - c == a + (b - c))
print(a * (b + c) == a * b + a * c)False
FalseAll’inizio questa realizzazione può apparire sorprendente, ma se ci pensate bene è una conseguenza abbastanza diretta di ciò che abbiamo detto fino ad ora. Se ad ogni operazione in virgola mobile, potenzialmente, facciamo approssimazioni, allora è chiaro che il risultato finale può dipendere dall’ordine in cui facciamo le cose. Nell’esempio che abbiamo appena visto il risultato dei due confronti è falso perché le due cose che confrontiamo differiscono di una quantità piccola—al limite della precisione della macchina—ma a volte la cosa può diventare distrastosa:
a = 1.e16
b = -1.e16
c = 1.
print((a + b) + c)
print(a + (b + c))1.0
0.0Questo è un caso interessante, perché chiaramente la prima risposta è quella che ci aspetteremmo ingenuamente, ma a seconda di come mettiamo le parentesi possiamo ottenere un nuemero completamente diverso. Il problema, chiaramente, è quello che abbiamo cercato di mettere in evidenza con l’esempio del numero di Avogadro nella sezione precedente, e la soluzione (se così possiamo dire) è cercare, ove possibile, di non sommare o sottrarre direttamente numeri che differiscono tra loro di molti ordini di grandezza.
11.3.4 Confronti tra numeri in virgola mobile
Ormai abbiamo capito che confrontare tra loro numeri in virgola mobile è potenzialmente pericoloso, eppure vi assicuro che, andando avanti, vi imbatterete di frequente in programmi del tipo
a = 0.1
b = 0.2
c = 0.3
if a + b <= c:
print("win!")
else:
print("lose...")lose...Se ci mettiamo per un attimo in un caso realistico, in cui a + b rappresenta una generica operazione matematica e a, b e c sono distribuiti in intervalli grandi rispetto alla precisione relativa dell’aritmetica, la cosa non è così disastrosa, perché il risultato della disuguaglianza sarà quello atteso nella maggior parte dei casi. Resta il fatto che, in senso stretto, non possiamo far affidamento sul confronto tra numeri in virgola mobile.
La libreria standard, oltre che pacchetti esterni come numpy, offrono funzioni specifiche per gestire il confronto tra numeri in virgola mobile, che possono essere utili in alcune situazioni
import math
import numpy as np
a = 0.1
b = 0.2
c = 0.3
print(math.isclose(a + b, c))
print(np.allclose(a + b, c))True
TrueAllora l’esempio iniziale di questa sezione si potrebbe scrivere in modo leggermente più robusto come
import math
a = 0.1
b = 0.2
c = 0.3
if a + b < c or math.isclose(a + b ,c):
print("win!")
else:
print("lose...")win!11.4 E quindi?
Semplice: l’aritmetica in virgola mobile è intrinsecamente inesatta e dovete sempre fare attenzione. Lo standard IEEE 754 è progettato in modo che le cose funzionino in modo trasparente nella maggior parte dei casi, ma non è magico, per cui non dovete assumere che i numeri in un calcolatore digitale si comportino come quelli sul libro di analisi. In particolare dovete fare attenzione a:
- addizioni e sottrazioni di numeri che differiscono tra loro di molti ordini di grandezza;
- sottrazioni di numeri molto vicini tra loro;
- moltiplicazioni e divisioni in catena in cui i risultati intermedi sono potenzialmente soggetti a underflow o overflow;
- accumulo coerente di piccoli errori dentro loop molto lunghi (se avete ancora tempo da dedicare all’argomento e vi suggerisco la lettura di questo report su un avvenimento drammatico in un teatro di guerra).
A volte ci sono espedienti che permettono di aggirare il problema—ad esempio passare ai logaritmi come abbiamo visto nel caso della distribuzione di Poisson. Altre volte le potenziali soluzione sono più sottili, come nel caso dell’algoritmo di Welford per il calcolo della varianza (Welford 1962). Ove possibile, usare librerie popolari e ben testate come numpy e scipy è sicuaramente da preferire rispetto al reimplementare tutto from scratch.
La lettura di questo capitolo non vi rende esperte di calcolo numerico, ma siete state avvertite!
Fan, X., T. G. Myers, B. A. D. Sukra, e G. Gabrielse. 2023. «Measurement of the Electron Magnetic Moment». Phys. Rev. Lett. 130 (febbraio): 071801. https://doi.org/10.1103/PhysRevLett.130.071801.
Marshall, Mason C., Daniel A. Rodriguez Castillo, Willa J. Arthur-Dworschack, et al. 2025. «High-Stability Single-Ion Clock with \(5.5\times10^{-19}\) Systematic Uncertainty». Phys. Rev. Lett. 135 (luglio): 033201. https://doi.org/10.1103/hb3c-dk28.
Welford, B. P. 1962. «Note on a Method for Calculating Corrected Sums of Squares and Products». Technometrics 4 (3): 419–20. https://doi.org/10.1080/00401706.1962.10490022.