11 - Numeri in virgola mobile: perché 0.1 + 0.2 non fa 0.3

Python per i dati, parte 4: i numeri decimali in un computer (IEEE 754) e perché 0.1 + 0.2 != 0.3

Corso Machine Learning Engineer · Modulo 01 · Python per i dati · Lezione 4

Scrivi 0.1 + 0.2 == 0.3 e Python risponde False. Non è un difetto del linguaggio e non lo si risolve cambiando libreria: è la conseguenza di come ogni computer scrive i numeri con la virgola, e nel lavoro sui dati si manifesta in forme molto meno innocue di un confronto che fallisce. Questa lezione segue liste, tuple, dizionari e insiemi e funzioni, file, eccezioni e oggetti: dai bit ricavi perché 0.1 non esiste in binario, calcoli la distanza fra due numeri vicini e incontri tre casi del mestiere in cui l'aritmetica in virgola mobile rompe un algoritmo corretto sulla carta, cioè la varianza con la formula «da libro», la somma di molti termini e la softmax.

Il luogo comune dice che i float sono approssimati ma l'errore è trascurabile. Nella varianza calcolata con la formula comoda l'errore produce un risultato negativo, una grandezza che per definizione non può esserlo: lo vedrai eseguito più sotto, con l'uscita vera.

Cos'è, perché esiste, quando NON basta

Il problema

Un computer ha memoria finita, i numeri reali sono infiniti. Bisogna scegliere un modo di rappresentare un sottoinsieme di reali in un numero fisso di bit. Lo standard IEEE 754 (1985, rivisto nel 2008 e nel 2019) è la risposta universale: quasi ogni processore lo implementa, quindi le operazioni di base (somma, prodotto, divisione) danno lo stesso risultato su macchine diverse.

L'analogia: la notazione scientifica con spazio limitato

Immagina di poter scrivere i numeri solo con 3 cifre significative e un esponente: 1.23 × 10⁴. Rappresenti bene 12300, ma non 12345 (diventa 1.23 × 10⁴ = 12300) e non puoi distinguere 12300 da 12349. Più il numero è grande, più è grande il salto tra due numeri rappresentabili. Un float è questo, in binario, con 53 cifre significative.

Quando i float NON bastano

• Denaro: usa interi (centesimi) o decimal.Decimal.
• Uguaglianza esatta tra risultati di calcoli: usa tolleranze.
• Interi molto grandi (oltre 2⁵³): un float li perde; usa int, che ha precisione arbitraria.
• Calcolo simbolico esatto (frazioni): fractions.Fraction.

Il filo comune è che il float è una buona scelta per grandezze misurate, dove l'incertezza della misura è già più grande dell'errore di arrotondamento, e una scelta cattiva per quantità contate o confrontate con esattezza. Per il denaro il problema è doppio: i centesimi non sono rappresentabili in binario e l'errore cresce con l'importo, come mostra l'esercizio 2 più sotto.

Preparazione

Serve solo NumPy; il resto è libreria standard (struct, math, decimal, fractions). Gli esempi sono stati eseguiti due volte ciascuno, con la stessa uscita, su Python 3.14.4 e NumPy 2.3.5. Il comando di installazione sotto è uno schema da adattare al tuo ambiente (non eseguito qui, dove NumPy era già presente).

python -m pip install numpy

Primo esempio: il famoso 0.1 + 0.2

import math

print(0.1 + 0.2)
print(0.1 + 0.2 == 0.3)
print(math.isclose(0.1 + 0.2, 0.3))

Output: 0.30000000000000004, False, True.

Python non sbaglia: stampa la stringa più corta che identifica il float ottenuto dal calcolo, e quel float non è 0.3. Il terzo confronto, con math.isclose, è il rimedio che sviluppiamo più avanti. Per vedere perché il risultato è proprio quello, guardiamo i bit.

Teoria: il formato float64, ricavato

Struttura (IEEE 754 binary64)

Un float di Python occupa 64 bit, divisi in tre campi:

Per i numeri «normali» il valore è

x = (−1)^s × (1 + f/2⁵²) × 2^(e − 1023)

Il «1 +» davanti alla frazione è il bit nascosto: in forma normalizzata la mantissa ha sempre la forma 1.xxxx… in binario, quindi quel primo 1 non si memorizza. Questo regala una cifra gratis: 52 bit memorizzati, 53 di precisione.

Guardiamo i bit di 0.1

import struct

def campi_float64(x):
    """Restituisce (segno, esponente, frazione) come stringhe di bit."""
    n = struct.unpack(">Q", struct.pack(">d", x))[0]
    b = f"{n:064b}"
    return b[0], b[1:12], b[12:]

s, e, f = campi_float64(0.1)
print(s, e, f)
print("esponente:", int(e, 2), "-> 2^", int(e, 2) - 1023)
mantissa = 1 + int(f, 2) / 2**52
print(mantissa * 2 ** (int(e, 2) - 1023))

Output: segno 0, esponente 01111111011, frazione 1001100110011001100110011001100110011001100110011010; poi esponente: 1019 -> 2^ -4; poi 0.1.

Come leggere il codice. struct.pack(">d", x) converte il float in 8 byte secondo lo standard IEEE 754 (il segno > fissa l'ordine dei byte, d significa «double», cioè 64 bit); struct.unpack(">Q", ...) reinterpreta gli stessi 8 byte come un intero senza segno a 64 bit. La stringa f"{n:064b}" scrive quell'intero in binario con 64 cifre. Poi il programma ritaglia i tre campi della figura: il primo bit è il segno, i successivi 11 sono l'esponente, gli ultimi 52 la frazione. L'ultima riga ricostruisce il valore con la formula del formato (mantissa × 2^(e − 1023)), e dà 0.1: la formula è davvero ciò che il computer fa.

Leggi: e = 1019, quindi il fattore è 2^(1019−1023) = 2⁻⁴ = 1/16. La mantissa è circa 1.6, e 1.6/16 = 0.1. Ma osserva la frazione: 1001 1001 1001 … si ripete. Perché?

Derivazione: 0.1 in binario non finisce mai

Per convertire una frazione decimale in binario si moltiplica ripetutamente per 2 e si prende la parte intera come cifra successiva.

Al passo 6 il resto è di nuovo 0.2, già visto al passo 2: il ciclo 0011 si ripete per sempre. Quindi

0.1₁₀ = 0.0001 1001 1001 1001 …₂

(i bit 1001 1001 … che hai visto nella frazione sono questo stesso ciclo, letto a partire dal primo 1 dopo aver spostato la virgola di quattro posti per l'esponente −4; l'ultimo gruppo è 1010 invece di 1001 perché il taglio arrotonda per eccesso). Un numero ha un'espansione binaria finita se e solo se, scritto come frazione irriducibile p/q, il denominatore q è una potenza di 2. 0.1 = 1/10 = 1/(2·5): contiene un fattore 5, quindi l'espansione è infinita. Esattamente come 1/3 = 0.3333… in decimale, perché 3 non divide 10.

Il computer deve tagliare dopo 53 cifre significative e arrotondare. Il numero che ottiene è leggermente diverso da 0.1. Quanto?

from decimal import Decimal
from fractions import Fraction

print(Decimal(0.1))
print(Fraction(0.1))

Output: 0.1000000000000000055511151231257827021181583404541015625 e 3602879701896397/36028797018963968. Il float 0.1 è in realtà questa frazione, il cui denominatore è 36028797018963968 = 2⁵⁵: una potenza di 2, come deve essere. Ecco perché 0.1 + 0.2 non dà esattamente 0.3: sommi due numeri che sono già stati arrotondati, e la somma viene arrotondata di nuovo.

Non associativa: l'ordine delle somme conta

Un fatto sorprendente conseguente dell'arrotondamento: l'addizione tra float non è associativa. In matematica (a + b) + c = a + (b + c); in virgola mobile no, perché ogni somma arrotonda il risultato intermedio.

import math
import numpy as np

a, b, c = 1e16, -1e16, 1.0
print("(a + b) + c =", (a + b) + c)
print("a + (b + c) =", a + (b + c))

rng = np.random.default_rng(0)
x = rng.random(1_000_000)
ordinamenti = {
    "ordine originale": x,
    "crescente": np.sort(x),
    "decrescente": np.sort(x)[::-1],
}
esatta = math.fsum(x)                                    # somma accurata, con riferimento affidabile
for nome, v in ordinamenti.items():
    s = float(np.sum(v))
    print(f"{nome:<18}: {s!r}  differenza dalla somma accurata: {s - esatta:+.2e}")

Output ottenuto:

(a + b) + c = 1.0
a + (b + c) = 0.0
ordine originale  : 500159.2564636844  differenza dalla somma accurata: -5.82e-11
crescente         : 500159.25646368443  differenza dalla somma accurata: +0.00e+00
decrescente       : 500159.25646368443  differenza dalla somma accurata: +0.00e+00

Come leggere il codice. La prima parte mostra il caso estremo: 1e16 + (-1e16) vale 0 e poi si somma 1 (risultato 1), ma sommando prima -1e16 + 1 si ottiene di nuovo -1e16 (l'uno viene perso perché i float attorno a 10¹⁶ sono distanti 2) e il risultato finale è 0. Stessi numeri, stesso insieme di operazioni, risultati diversi.

La seconda parte somma un milione di numeri casuali nei tre ordini; math.fsum fornisce un riferimento ottenuto con una procedura che tiene traccia degli errori parziali (restituisce il float più vicino alla somma esatta). Le differenze sono minuscole: l'ordine originale si discosta dalla somma accurata di circa 6·10⁻¹¹, mentre i due ordinamenti coincidono con essa; con un altro array, o con un altro seme, la situazione sarebbe diversa. Ciò che conta è che il risultato dipende dall'ordine, anche se di pochissimo. Questo è il motivo per cui due esecuzioni dello stesso addestramento su GPU, che sommano in ordini diversi, possono dare risultati con differenze nelle ultime cifre, e per cui un test di uguaglianza esatta fra due calcoli «equivalenti» non è mai affidabile.

Epsilon macchina e spaziatura

Tra due potenze consecutive di 2, 2^k e 2^(k+1), ci sono esattamente 2⁵² numeri float equispaziati (i valori possibili della frazione). Quindi la distanza tra due float vicini in quell'intervallo è

ulp(x) = 2^(k−52) per 2^k ≤ x < 2^(k+1)

(ulp = unit in the last place). In particolare, sopra 1.0: ulp(1) = 2⁻⁵² ≈ 2.22 × 10⁻¹⁶. Questo numero si chiama epsilon macchina.

Arrotondando al più vicino, l'errore massimo è metà della spaziatura, quindi l'errore relativo di una singola operazione è al massimo u = 2⁻⁵³ ≈ 1.11 × 10⁻¹⁶: circa 16 cifre decimali di precisione.

import math
import sys

print(sys.float_info.epsilon, 2**-52, math.ulp(1.0))
print(math.ulp(1e16), math.ulp(0.0))
print(2.0**53 == 2.0**53 + 1)

Output: tre volte 2.220446049250313e-16; poi 2.0 e 5e-324; poi True.

Interpreta le ultime tre righe:

• math.ulp(1e16) vale 2.0: attorno a 10¹⁶ i float sono distanti 2. Quindi 1e16 + 1 non è rappresentabile e viene arrotondato.
• math.ulp(0.0) è 5e-324, il più piccolo float positivo (denormalizzato).
• 2.0**53 == 2.0**53 + 1 è True: da 2⁵³ in poi i float non rappresentano più tutti gli interi. Per questo i grandi interi (ID, timestamp in nanosecondi) non vanno mai convertiti in float.

La spaziatura relativa resta circa costante (tra 2⁻⁵³ e 2⁻⁵²); quella assoluta cresce con il numero. I float sono molto precisi vicino allo zero e sempre meno precisi lontano.

Un secondo esempio, contrastante: float32 e float16

Il float64 dà circa 16 cifre decimali. Le reti neurali spesso usano float32 (circa 7 cifre) o float16 (circa 3). Vediamo che cosa si perde, con due esperimenti.

import numpy as np

# 1) accumulare 20 milioni di uni, in float32
uni = np.ones(20_000_000, dtype=np.float32)
print("cumsum in float32 (sequenziale):", float(np.cumsum(uni)[-1]))
print("np.sum in float32 (a coppie)   :", float(uni.sum()))
print("np.sum con accumulatore float64:", float(uni.sum(dtype=np.float64)))

# 2) i limiti di float16
with np.errstate(over="ignore"):
    print("float16: 60000 * 2 =", np.float16(60000) * np.float16(2))
print("float16: 0.1 diventa", float(np.float16(0.1)), "| float32: 0.1 diventa", float(np.float32(0.1)))
print("float32: 16777216 + 1 =", float(np.float32(16777216) + np.float32(1)))

Output ottenuto:

cumsum in float32 (sequenziale): 16777216.0
np.sum in float32 (a coppie)   : 20000000.0
np.sum con accumulatore float64: 20000000.0
float16: 60000 * 2 = inf
float16: 0.1 diventa 0.0999755859375 | float32: 0.1 diventa 0.10000000149011612
float32: 16777216 + 1 = 16777216.0

Come leggere il codice. Nella prima parte, np.cumsum somma in sequenza: quando il totale parziale raggiunge 2²⁴ = 16 777 216, in float32 la spaziatura tra numeri consecutivi diventa 2, quindi aggiungere 1 non cambia più il totale, che resta bloccato a 16 777 216 invece di arrivare a 20 milioni. np.sum, invece, usa internamente una somma a coppie (somma gruppi piccoli e poi li combina), che mantiene gli addendi di dimensione simile e produce il risultato corretto di 20 milioni: l'algoritmo di somma conta quanto la precisione del tipo. Con dtype=np.float64 l'accumulatore ha un margine enorme. Nella seconda parte, float16 non sa rappresentare 120000 (il massimo è 65504) e restituisce inf; rappresenta 0.1 come 0.09997... (errore dello 0.02%), mentre float32 come 0.10000000149.... La lezione pratica per il deep learning: le precisioni ridotte accelerano e risparmiano memoria ma vanno usate con accumulatori più precisi e con attenzione ai limiti dell'intervallo rappresentabile (da cui la scalatura della perdita, di cui trovi la definizione nell'elenco più sotto).

Tre casi in cui l'algoritmo si rompe

Caso 1: la cancellazione catastrofica

Sottrarre due numeri quasi uguali azzera le cifre significative in comune e lascia solo il «rumore» dei bit meno significativi.

print((1e16 + 1) - 1e16)
print((1e16 + 2) - 1e16)

Output: 0.0 e 2.0. Il primo dovrebbe valere 1: 1e16 + 1 è stato arrotondato a 1e16, e la differenza è zero. Il secondo vale 2 perché 1e16 + 2 è rappresentabile.

Applicazione: la varianza «da libro»

La formula comoda Var(X) = E[X²] − E[X]² è algebricamente corretta (la ricavi dalle proprietà del valore atteso). In virgola mobile, se i dati hanno una media grande rispetto alla dispersione, i due termini sono enormi e quasi uguali, e la differenza è dominata dal rumore.

import numpy as np

rng = np.random.default_rng(0)
x = rng.normal(0, 1, 1000) + 1e9          # dispersione ~1, media ~1 miliardo

var_ingenua = np.mean(x**2) - np.mean(x) ** 2
var_due_passi = np.mean((x - x.mean()) ** 2)
print(var_ingenua, var_due_passi, np.var(x))

Output: -128.0, 0.9540463497672959, 0.9540463497672959. La formula ingenua restituisce una varianza negativa: E[X²] vale circa 10¹⁸ e a quella scala i float sono distanti 128 (math.ulp(1e18)), quindi la differenza è dominata dall'arrotondamento. Il metodo a due passi (prima sottrai la media, poi quadra) sottrae numeri piccoli e funziona.

Lezione di ingegneria: due formule uguali sulla carta non sono uguali sul computer. NumPy usa il metodo a due passi.

Caso 2: sommare tanti numeri

import math

import numpy as np

def somma_a_mano(valori):
    s = 0.0
    for v in valori:
        s += v
    return s

print(somma_a_mano([0.1] * 10))
print(somma_a_mano([0.1] * 1_000_000))
print(sum([0.1] * 1_000_000), math.fsum([0.1] * 1_000_000))
print(np.sum(np.full(1_000_000, 0.1)))

Output: 0.9999999999999999, 100000.00000133288, 100000.0 100000.0, 100000.00000000003.

Il ciclo naïf accumula errore: ogni somma arrotonda, e gli errori si sommano (100000.00000133288 invece di 100000). Le alternative sono tre. math.fsum somma con precisione estesa e restituisce il risultato correttamente arrotondato: costa di più, ma è l'arbitro quando serve un riferimento. La sum() di Python dalla versione 3.12 usa una sommatoria compensata (algoritmo di Neumaier) per i float, come dicono le note di rilascio della 3.12, e da qui viene il 100000.0; prima si comporta come il ciclo a mano.

np.sum usa una sommatoria a coppie parziale: secondo la documentazione di numpy.sum è sempre applicata quando non indichi un asse, mentre lungo un asse dipende da quale. Somma a blocchi in modo ricorsivo e riduce la crescita dell'errore, nel caso peggiore, da proporzionale a n a proporzionale a log n. Qui non dà esattamente 100000, ma l'errore (circa 3e-11) è oltre quattro ordini di grandezza più piccolo di quello del ciclo (circa 1.3e-6).

L'algoritmo di Kahan (1965) tiene una variabile di compensazione che raccoglie l'errore di arrotondamento di ogni passo e lo reintroduce al successivo:

def kahan(valori):
    somma, compensazione = 0.0, 0.0
    for v in valori:
        y = v - compensazione
        t = somma + y
        compensazione = (t - somma) - y     # ciò che è andato perso nell'addizione
        somma = t
    return somma

print(kahan([0.1] * 1_000_000))

Output: 100000.0. Attenzione: un compilatore troppo «intelligente» può semplificare algebricamente (t − somma) − y in zero e disattivare la compensazione; in Python non accade.

Trade-off. Fsum e Kahan sono più lenti della somma semplice. Per np.sum su array normali la sommatoria a coppie è già un ottimo compromesso. Serve cura speciale solo con moltissimi termini, valori di scala molto diversa, o dove il risultato esatto conta.

Caso 3: softmax

La funzione softmax trasforma un vettore di punteggi z in probabilità, ed è all'uscita di quasi ogni classificatore a più classi:

softmax(z)ᵢ = e^(zᵢ) / Σⱼ e^(zⱼ)

Con punteggi grandi, e^z va in overflow: np.exp(710) vale inf (il massimo float64 è circa 1.8 × 10³⁰⁸, e e^709 ≈ 8.2 × 10³⁰⁷ è ancora rappresentabile). Allora inf/inf = nan.

Rimedio e dimostrazione di correttezza. Per qualunque costante c:

e^(zᵢ − c) / Σⱼ e^(zⱼ − c) = ( e^(zᵢ) · e^(−c) ) / ( e^(−c) · Σⱼ e^(zⱼ) ) = e^(zᵢ) / Σⱼ e^(zⱼ)

perché e^(−c) è un fattore comune e si semplifica. Il risultato non cambia. Scegliendo c = max(z) l'esponente più grande diventa e⁰ = 1 e nessun termine supera 1: niente overflow.

import numpy as np

def softmax_ingenuo(z):
    e = np.exp(z)
    return e / e.sum()

def softmax_stabile(z):
    e = np.exp(z - z.max())
    return e / e.sum()

z = np.array([1000.0, 1001.0, 1002.0])
with np.errstate(all="ignore"):
    print(softmax_ingenuo(z), softmax_stabile(z))
piccolo = np.array([1.0, 2.0, 3.0])
print(softmax_ingenuo(piccolo), softmax_stabile(piccolo))

Output: [nan nan nan] contro [0.09003057 0.24472847 0.66524096]; sull'esempio piccolo i due coincidono. Questo trucco (e il suo parente per il logaritmo, log-sum-exp) è nel codice di ogni libreria di deep learning. Se un giorno una tua rete restituisce nan, sospetta prima di tutto un overflow o una divisione per zero.

Come ci si difende

Confrontare con tolleranza

math.isclose(a, b, rel_tol=1e-9, abs_tol=0.0) restituisce vero se

|a − b| ≤ max( rel_tol · max(|a|, |b|), abs_tol )

La tolleranza relativa va bene per numeri lontani da zero; vicino a zero è inutile (rel_tol · 0 = 0), quindi serve una tolleranza assoluta:

import math
import numpy as np

print(math.isclose(1e-12, 2e-12))                  # False: relativamente sono il doppio
print(math.isclose(1e-12, 2e-12, abs_tol=1e-9))    # True: assolutamente sono vicinissimi
print(np.allclose([1.0, 2.0], [1.0 + 1e-9, 2.0]))  # True

Output: False, True, True. Per gli array: np.isclose e np.allclose. Nei test: np.testing.assert_allclose.

Denaro e decimali esatti

from decimal import Decimal

print(Decimal("0.1") + Decimal("0.2"), Decimal("0.1") + Decimal("0.2") == Decimal("0.3"))

Output: 0.3 True. Decimal ragiona in base 10, quindi 0.1 è esatto. Attenzione: costruisci da stringa, non da float (Decimal(0.1) eredita l'errore, come visto sopra). In pratica, per il denaro spesso basta usare interi in centesimi.

Arrotondamento

round() di Python arrotonda .5 al pari più vicino (round half to even), non sempre per eccesso: round(0.5) = 0, round(1.5) = 2, round(2.5) = 2, round(3.5) = 4. È la scelta dello standard, che evita una deriva sistematica nelle somme di molti arrotondamenti. Inoltre round(2.675, 2) vale 2.67: perché il float vicino a 2.675 è in realtà leggermente sotto 2.675.

Float a 32 e 16 bit: il machine learning vive qui

(Valori ottenibili da np.finfo(tipo); np.finfo(np.float16).max vale 65504.)

Le reti neurali si addestrano spesso in float32 o in formati a 16 bit, perché la memoria e la velocità contano più della precisione: i pesi sono rumorosi comunque. Conseguenze da conoscere:

• Con float32, 16777216 + 1 == 16777216 (a 2²⁴ i float32 non distinguono più gli interi consecutivi). Gli accumulatori vengono spesso tenuti a precisione più alta.
• In float16 il massimo è 65504: un valore più grande va in inf. Nasce da qui la tecnica della scalatura della perdita (loss scaling) nell'addestramento a precisione mista.
• Il comportamento cambia con l'ordine delle operazioni: addestrare due volte la stessa rete su GPU può dare risultati leggermente diversi perché le somme parallele sono in ordini diversi (e la somma in virgola mobile non è associativa: (a + b) + c ≠ a + (b + c) in generale).

Le tre conseguenze hanno la stessa radice: meno bit significano uno spazio fra numeri più largo e un massimo più basso. Il rimedio non è «usa sempre float64», che costa memoria e velocità, ma scegliere dove tenere la precisione (nell'accumulatore e non nei pesi) e dove proteggere l'intervallo dei valori.

Valori speciali. inf (overflow), -inf, nan (risultato non definito: 0.0 * inf). nan non è uguale a niente, nemmeno a se stesso: float("nan") == float("nan") è False. Per testarlo usa math.isnan o np.isnan. Un nan si propaga: ogni operazione con un nan dà nan, quindi un solo valore corrotto può inquinare un'intera metrica.

Un piccolo protocollo di debug numerico

Quando un calcolo dà un risultato strano (un nan, una varianza negativa, una metrica che cambia a ogni esecuzione), conviene un ordine di indagine, dal più economico al più costoso. Primo: individua dove compare il primo valore anomalo (np.isnan(x).argmax() dà la posizione del primo nan in un array); l'errore nasce lì, tutto il resto è propagazione. Secondo: stampa gli ordini di grandezza delle quantità coinvolte (minimi, massimi, np.log10(np.abs(x).max())): una sottrazione fra due numeri con ordini di grandezza simili e molto più grandi del risultato è un candidato alla cancellazione.

Terzo: rifai lo stesso calcolo in un tipo più preciso (float64 al posto di float32, o Decimal/Fraction per piccoli esempi): se il risultato cambia in modo visibile, il problema è numerico e non logico. Quarto: riscrivi la formula in una forma stabile (due passi, sottrazione del massimo, log1p, fsum) e verifica che coincida con la versione ingenua su dati piccoli. Quinto: sostituisci ogni == fra float con una tolleranza dichiarata. Seguire sempre lo stesso ordine evita di perdere ore a cercare un bug logico che non c'è.

Errori tipici

Best practice

• Mai == tra float calcolati. Sempre una tolleranza dichiarata, con giustificazione.
• Normalizza prima di sottrarre valori grandi; preferisci le formule numericamente stabili (due passi, log-sum-exp).
• Usa le funzioni delle librerie (np.var, np.sum, np.logaddexp) invece di riscrivere formule: sono state studiate per la stabilità.
• Dichiara il tipo (dtype) quando importa; non fidarti dei default per dati sensibili.
• Controlla i nan (np.isnan(x).any()) all'ingresso e dopo i passaggi critici.
• Denaro: interi o Decimal.

Le sei regole discendono dallo stesso fatto: il float non è il numero reale, è il suo rappresentante più vicino. Funzioni di libreria, tolleranze dichiarate e controllo dei nan sono tre modi di non fidarti di quel rappresentante senza averlo guardato. Decidi la politica una volta per progetto e scrivila nel README, così non resta nella testa di una persona sola.

Mini-progetto

Costruisci un piccolo laboratorio numerico lab_float.py:

1. Una funzione ulp_attorno(x) che confronta math.ulp(x) con la formula 2**(floor(log2(x)) - 52) per x in [1e-5, 1e20] e verifica che coincidano.
2. Una funzione somma_confronto(n) che restituisce gli errori assoluti (rispetto a math.fsum) di ciclo naïf, np.sum, Kahan per n valori casuali.
3. Un confronto della varianza ingenua e stabile al variare dell'offset 10**k, k = 0..12, con grafico dell'errore relativo in funzione di k.
4. Un test che mostra che (a + b) + c != a + (b + c) per tre valori scelti da te.

Il primo punto verifica la formula dell'ulp contro la funzione di libreria, il secondo e il terzo ripetono i due casi della lezione su dati tuoi, il quarto ti obbliga a trovare tre numeri che violano l'associatività. Per il grafico del terzo punto serve Matplotlib, trattato in visualizzazione con Matplotlib. Non c'è una soluzione ufficiale: il criterio è che ogni funzione stampi un numero che avevi previsto prima di eseguirla.

Challenge

• Easy. Trova a mano la rappresentazione binaria di 0.375 e verifica con struct che la frazione termini.
• Medium. Scrivi log_somma_exp(z) stabile: log Σ e^(zᵢ) calcolata come c + log Σ e^(zᵢ − c) con c = max(z). Dimostra l'identità.
• Hard. Implementa l'algoritmo di Welford per la varianza in streaming e confrontalo con l'ingenuo sulla stessa sfida dei dati con offset 10⁹.
• Expert. Mostra che la somma in virgola mobile non è associativa con un esempio minimo trovato per ricerca casuale. Poi calcola, per un array di n numeri, il limite superiore dell'errore di sommatoria ricorsiva e di sommatoria a coppie (Higham, cap. 4).
• Enterprise. Definisci per il tuo team una politica numerica: quando si usa float32 e quando float64, quali tolleranze nei test, come si rileva un nan in produzione.

I cinque livelli salgono da un controllo a mano a una decisione di squadra. Easy e Medium si risolvono con ciò che hai visto qui; Hard e Expert chiedono di consultare la letteratura citata in fondo; Enterprise non ha una risposta di codice ma un documento di una pagina che i colleghi possano applicare.

Che cosa sai fare adesso

Sai leggere i tre campi di un float64 e ricostruire il valore con la formula del formato, sai spiegare perché 0.1 non ha espansione binaria finita (il denominatore ha un fattore 5) e stimare con l'ulp la distanza fra due float vicini. Riconosci i tre punti in cui un algoritmo corretto sulla carta si rompe, cioè cancellazione, somma di molti termini e overflow nell'esponenziale, e sai quale rimedio usare in ciascuno. Il passo successivo è NumPy: array, vettorizzazione e broadcasting, che porta tutto questo su milioni di numeri alla volta.

Nella stessa serie trovi anche: 13 - Pandas: tabelle di dati senza sorprese, 21 - Analisi matematica 1: limiti e derivate, 25 - Calcolo matriciale e differenziazione automatica: dal Jacobiano a un autograd scritto da zero, 27 - Entropia, entropia incrociata e divergenza KL: perché le funzioni di costo sono tutte la stessa cosa (che richiamano questa lezione).

Fonti

Riferimenti bibliografici citati nel testo, non consultati online in questa revisione: IEEE Std 754-2019; D. Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, ACM Computing Surveys 23(1), 1991; N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2ª ed., SIAM, 2002; W. Kahan, Communications of the ACM 8(1), 1965. Le misure in questo articolo (ulp, epsilon, uscite) sono state ricalcolate con Python 3.14.4 e NumPy 2.3.5.