Numeri macchina^{1}

A. Sommariva^{2}

Keywords: Errore di troncamento e di arrotondamento, rappresentazione floating-point dei reali, precisione di macchina; operazioni aritmetiche con numeri approssimati, potenziale instabilita’ della sottrazione, stabilizzazione della formula risolutiva per equazioni di secondo grado, condizionamento di funzioni, propagazione degli errori in algoritmi iterativi (successione di Archimede per il calcolo di pigreco),esempi di costo computazionale degli algoritmi numerici: schema di Hoerner per i polinomi, potenza rapida tramite codifica binaria dell’esponente, calcolo della funzione exp(x) con scalatura della variabile e formula di Taylor, calcolo di determinanti.

Revisione: 24 dicembre 2018

1. Numeri macchina e loro proprietà

Fissato un numero naturale β > 1 è possibile vedere che ogni numero reale x ≠ 0 ha una unica rappresentazione del tipo

x = sign(x)β^{e} ∑_{k=1}^{+∞} d_{k}β^{−k}, 0 < d_{1}, 0 ≤ d_{k} ≤ β − 1, e ∈ ℤ  (1)

dove

sign(y) = \{ \begin{matrix}1, y > 0 \\ −1, y < 0 \\ 0, y = 0. \end{matrix}

è la funzione segno.

Dall’avvio dello standard IEEE 754-1985, in cui venne realizzato lo standard per l’aritmetica binari floating point, la base più utilizzata risulta β = 2, seppure col nuovo standard 854 si stiano avviando studi anche relativamente a β = 10. [? ].

Vediamo quale esempio come rappresentare il numero π al variare di due classiche scelte di β:

• base 10: 3, 141592653589793238462643383279502884197169399 … secondo la rappresentazione sopra descritta e’

+1 · 10^{1} · (3 · 10^{−1} + 1 · 10^{−2} + 4 · 10^{−2} + …)

che viene spesso denotato come

π = (+1)·(0.3141592653589793238462643383279502884197169399 …)_{10}·10^{1};

• base 2: 11, 0010010000111111011010101000100010000101 … secondo la rappresentazione sopra descritta e’

+1·2^{2}·(1·2^{−1}+1·2^{−2}+0·2^{−3} …+0·2^{−4}+1·2^{−5} …) = +4(1/2+1/4+1/32+…),

che viene spesso denotato come

π = (+1) · (0.110010010000111111011010101000100010000101 …)_{2} · 2^{2}.

La particolarità della sommatoria in (??), ovvero

x = sign(x)β^{e} ∑_{k=1}^{+∞} d_{k}β^{−k}, 0 < d_{1}, 0 ≤ d_{k} ≤ β − 1, e ∈ ℤ  (2)

è ovviamente tale che x potrebbe essere non rappresentabile esattamente dal calcolatore se esistono infiniti d_{k} > 0, in quanto il calcolatore è capace di immagazzinare solo un numero finito di tali d_{k}.

Ad esempio, non è difficile vedere che:

• i numeri irrazionali come π, indipendentemente dalla base β, hanno sempre un numero infinito di cifre d_{k};

• in base 10 alcuni numeri razionali hanno comunque un numero infinito di cifre d_{k} (si pensi ai numeri periodici come 1/3 = 0.333 …, ma comunque in certe basi possono avere una rappresentazione finita (nel nostro esempio 1/3 = 1·3^{−1} per β = 3).

Nota. 1.1. Osserviamo che storicamente, la scelta della base β è risultata molto eterogenea:

• la base utilizzata dalla maggior parte dei computer attuali è β = 2,

• il computer russo SETUN (≈ 1957) usava la base β = 3,

• il PDP-10 (≈ 1965), Burroughs 570 e 6700 utilizzavano la base ottale, cioè β = 8,

• la base che usiamo nella vita di tutti i giorni è β = 10, ed è stata una scelta comune anche per computer del passato come l’HP-11C, HP-15C (≈ 1980) o Eniac (≈ 1947),

• esistono popoli in Nigeria, India, Nepal, la cui è base β = 12,

• molti computers, come l’IBM 360, l’IBM 3033 (≈ 1977), usavano la base β = 16,

• i Maya usavano la base β = 20.

• in tavolette babilonesi si utilizzava la base β = 60.

1.1. Numeri floating-point

In virtú di di considerazioni sulla finitezza delle cifre immagazzibili, ipotesi necessaria perché il numero sia utilizzabile da un computer, risulta naturale la seguente

Definizione. L’insieme dei numeri macchina F (β, t, L, U) di ℝ è costituito da quei numeri y che sono 0 oppure sono rappresentati con il sistema a virgola mobile normalizzata o floating point normalizzato da

y = sign(y) · (0.d_{1} … d_{t})_{β}β^{e} = sign(y) · β^{e} ∑_{k=1}^{t} d_{k}β^{−k}, d_{1} > 0, 0 ≤ d_{k} ≤ β − 1

con

• β > 1 numero naturale detto base, t numero prefissato di cifre di mantissa,

• e ∈ ℤ l’esponente intero tale che L ≤ e ≤ U.

La parola normalizzato sottolinea che d_{1} > 0. Si può provare che ∑_{k=1}^{t} d_{k}β^{−k} < 1 e quindi e = floor(1 + log_{β}(x)).

Ogni elemento di F (β, t, L, U) è detto numero macchina, ed evidentemente è un numero razionale.

Nota. 1.2. Osserviamo che se e > 0 allora

y = sign(y) · (0.d_{1} … d_{t})_{β}β^{e} = sign(y) · β^{e} ∑_{k=1}^{t} d_{k}β^{−k}

= sign(y) · (\underbrace{∑_{k=1}^{e} d_{k}β^{e−k}}_{parte intera}+ \underbrace{∑_{k=e+1}^{t} d_{k}β^{e−k})}_{parte frazionaria}

dove l’ultima somma è nulla per definizione se t < e + 1.

Nota. 1.3. Si osservi che il sistema floating-point è solo una delle scelte possibili, in quanto una volta si utilizzava un sistema a virgola fissa, ovvero con un numero fissato di cifre prima della virgola e un certo numero di cifre dopo la virgola, abbandonato (nonostante presentasse alcuni vantaggi) perchè l’intervallo di valori rappresentati è modesto e la precisione dei numeri frazionari scarsa.

Di ogni numero macchina si immagazzinano • una cifra per il segno (ovvero 0 per positivo e 1 per negativo),

• le cifre dell’esponente,

• le cifre della mantissa.

segno esponente mantissa
1 cifra m cifre  t cifre

Ognuna di queste cifre è detta bit. Ogni numero è registrato in un vettore detto parola di un certo numero di M bits.

Nota. 1.4. Quello che proponiamo è un modello dei numeri floating-point. Le implementazioni sono molte più complesse e posso avere qualche trascurabile scostamento da quanto affermato.

Ad esempio, nella notazione binaria normalizzata, immagazzinate le cifre (1, d_{1}, …, d_{t−1})_{β} e l’esponente e, a volte si intende

y = sign(y) · (1.d_{2} … d_{t})_{β}β^{e} = sign(y) · β^{e}(1 + ∑_{k=2}^{t} d_{k}β^{−k})  (3)

e non come in precedenza

y = sign(y) · (1.d_{1} … d_{t−1})_{β}β^{e} = sign(y) · β^{e}(β^{−1} + ∑_{k=2}^{t} d_{k−1}β^{−k}). (4)

Questo potenzialmente può creare qualche confusione. Di seguito scegliamo la notazione più utilizzata nei manuali di analisi numerica, ovvero (??).

1.2. Alcune proprietá di F (β, t, L, U)

Vediamo di seguito, alcune proprietá di F (β, t, L, U).

Proposizione 1.5. Il piú piccolo numero macchina positivo (normalizzato) é β^{L−1}.

Dimostrazione. Ponendo d_{1} = 1, d_{2} = … = d_{t} = 0 ed e = L otteniamo che il più piccolo numero macchina (e quindi normalizzato!) positivo rappresentabile è β^{L−1}.

Proposizione 1.6. Il piú grande numero macchina positivo (normalizzato) é β^{U}(1 − β^{−t}), dove t é il numero di cifre della mantissa.

Dimostrazione. Ponendo d_{k} = β − 1, per k = 1, …, t ed e = U otteniamo che il massimo numero macchina rappresentabile è

β^{U}(1 − β^{−t}).

Proposizione 1.7. La somma relativa alla mantissa, cioè

∑_{k=1}^{t} d_{k}β^{−k},

appartiene all’intervallo (0, 1 − β^{−t}] e quindi positiva e strettamente minore di 1.

Dimostrazione. La dimostrazione è come la precedente, prendendo 0 invece di U.

Inoltre

Proposizione 1.8. L’insieme dei numeri macchina F (β, t, L, U) ha un numero finito di elementi e la sua cardinalità è

2(β − 1)β^{t−1}(U − L + 1) + 1.

Dimostrazione. Per una fissata mantissa, a t cifre, il numero piú piccolo che possiamo scrivere è

(1 0 … 0 0)_{β}

mentre il più grande è

(\underbrace{β − 1 β − 1 … β − 1 β − 1}_{t})_{β}.

Così non è difficile vedere che i numeri scrivibili con t cifre di mantissa sono quelli in cui la prima cifra va da 1 a β − 1 e le altre da 0 a β − a, e quindi la loro cardinalità risulta

(β − 1) · \underbrace{β … ββ}_{t−1} = (β − 1)β^{t−1}.

Quindi, al variare dell’esponente e in L, L + 1, …, U − 1, U, e del segno + o −, i numeri floating-point sono

2(β − 1)β^{t−1}(U − L + 1).

Proposizione 1.9. I numeri macchina sono numeri razionali. In particolare quelli che hanno esponente e ≥ t sono interi.

Dimostrazione. La scrittura dei numeri floating-point come somma di un numero finito di numeri frazionari implica che i numeri risultanti siamo razionali.

Per vedere se e ≥ t allora sono interi, da

x = sgn(x)β^{e} ∑_{i=1}^{t} d_{i}β^{−i} = sgn(x) ∑_{i=1}^{t} d_{i}β^{e−i}

visto che per i = 1, …, t ≥ e che β^{e−i} è intero, come d’altra parte i d_{i}, lo è anche x.

Proposizione 1.10. L’insieme dei numeri macchina F (β, t, L, U) non consiste di punti equispaziati.

Dimostrazione. Per capirlo, il numero positivo più piccolo, e quindi più vicino a 0 è β^{1−L} ≪ 1 mentre i numeri floating-point che hanno esponente e ≥ t sono interi e quindi la distanza di due tali numeri successivi è sicuramente maggiore o uguale 1. In definitiva, due numeri successivi vicino a 0 distano molto poco, e due numeri successivi con esponente molto grande distano molto, e quindi l’insieme dei floating-point non consiste di punti equispaziati.

Esempio. In Matlab, in precisione doppia, il numero successivo a 0 è circa 2.2251 · 10^{−308} mentre i due numeri piú grandi distano l’uno dall’altro ≈ 1.9958 · 10^{+292}. Di conseguenza l’insieme numerico del Matlab non consiste di punti equispaziati perchè 1.9958 · 10^{292} ≫ 2.2251 · 10^{−308}.

Proposizione 1.11. L’insieme dei numeri macchina F (β, t, L, U) non consiste di punti equispaziati.

Dimostrazione. Per capirlo, il numero positivo più piccolo, e quindi più vicino a 0 è β^{1−L} ≪ 1 mentre i numeri floating-point che hanno esponente e ≥ t sono interi e quindi la distanza di due tali numeri successivi è sicuramente maggiore o uguale 1. In definitiva, due numeri successivi vicino a 0 distano molto poco, e due numeri successivi con esponente molto grande distano molto, e quindi l’insieme dei floating-point non consiste di punti equispaziati.

Nota. 1.12. Nel caso della notazione

y = sign(y) · (1.d_{2} … d_{t})_{β}β^{e} = sign(y) · β^{e}(1 + ∑_{k=2}^{t} d_{k}β^{−k})  (5)

abbiamo i seguenti risultati per F (β, t, L, U), che si possono dimostrare con tecniche simili a quelle viste sopra:

• la somma degli elementi della mantissa (1, d_{1} … d_{t−1} è un numero nell’intervallo [1, 2 − β^{−(t−1)}] ⊂ [1, 2);

• Il più piccolo numero positivo è β^{L} (usare quale mantissa (1.0 … 0));

• Il più grande numero positivo è (2 − β^{−(t−1)}) · β^{U};

Gli altri asserti sono uguali a quelli visti per la notazione classica.

Matlab utilizza questa notazione e quindi per quanto riguarda ad esempio la precisione doppia F (2, 53, −1022, 1023):

• la somma degli elementi della mantissa (1, d_{1} … d_{52} è un numero nell’intervallo [1, 2 − β^{−(52)}] ≈ [1, 2 − 2.22 · 10^{−16};

• Il più piccolo numero positivo è β^{L} = 2^{−1022} ≈ 2.225073858507201e − 308; • Il più grande numero positivo è (2 − β^{−(t−1)}) · β^{U} = (2 − β^{−(52)}) · 2^{1023} =

1.797693134862316e + 308.

Di seguito faremo esempi di queste quantitá relativamente agli insiemi F (2, 24, −126, 127) e F (2, 53, −1022, 1023).

Nota. 1.13. Per quanto riguarda numeri speciali

• per rappresentare +inf si utilizza il segno pari a 0, tutti gli esponenti uguali a 1 e tutta la mantissa uguale a 0;

• per rappresentare -inf si utilizza il segno pari a 1, tutti gli esponenti uguali a 1 e tutta la mantissa uguale a 0;

• per rappresentare lo 0 si distingue +0 da −0 (generalmente utilizzati per denotare un numero troppo piccolo in modulo che viene arrotondato a 0, pur mantenendo il segno); a parte il segno (positivo di default), tutte le altre cifre di esponente e mantissa sono uguali a 0;

• per rappresentare con NaN il risultato di un’operazione (numerica) eseguita su operandi non validi (specialmente in calcoli in virgola mobile). Secondo IEEE 754 i NaN sono rappresentati con il campo dell’esponenziale riempito di ”1” e un numero diverso da zero nel campo della mantissa.

1.3. Precisione singola e doppia

Descriviamo di seguito i più comuni insiemi di numeri floating-point. Usualmente se M è il numero di bits necessari per un numero in precisione singola, nel caso di precisione doppia questo è usualmente composto da 2 parole di M bits.

Classici esempi sono per quanto riguarda i processori che seguono lo standard IEEE754 (inizialmente proposto da INTEL, studiato dal 1977 e introdotto nel 1985), e i successivi aggiornamenti fino all’attuale versione 2008 e alle proposte di IEEE854,

• per la precisione singola F (2, 24, −126, 127) in cui ogni numero macchina occupa 32 bit (o 4 byte essendo un byte pari a 8 bit)

• per la doppia F (2, 53, −1022, 1023) in cui ogni numero macchina occupa 64 bit.

Nota. 1.14. Ricordiamo che nella implementazione effettiva di IEEE754 ci sono molti dettagli tecnici che la rendono leggermente diversa da quanto analizzato, anche se essenzialmente i risultati sono gli stessi del nostro modello scelto.

1.3.1. Precisione singola Nel caso della precisione singola F (2, 24, −126, 127),

• Si usa la base 2.

• Si usa un bit per il segno.

• Si usa una mantissa di 24 cifre binarie, ma visto che si utilizza la notazione normalizzata, solo 23 sono necessarie.

• Visto che si compone di 32 bits, di cui 24 per mantissa e segno, 8 possono essere dedicate all’esponente; quindi si possono rappresentare numeri interi da 0 a 255; in questo ambito i numeri 0 e 255 vengono utilizzati per significati speciali e quindi utilizzati esclusivamente gli interi N_{e} da 1 a 254; fissato bias=127 = 2^{7} − 1 per ottenere gli esponenti e li si codifica con N_{e} tale che

e = N_{e} − bias .

Evidentemente il piú piccolo esponente e_{min} = 1 − 127 = 126, mentre il piú grande é e_{min} = 254 − 127 = 127.

Nell’implementazione Matlab,

• Il minimo numero positivo in precisione singola é x_{min} = 1.1754944e − 38;

• Il massimo numero positivo in precisione singola é x_{max} = 3.4028235e + 38;

• la cardinalità dell’insieme F (2, 24, −126, 127) è

≈ 4.2950e + 09.

1.3.2. Precisione doppia Nel caso della precisione singola F (2, 53, −1022, 1023),

• Si usa la base 2.

• Si usa un bit per il segno.

• Si usa una mantissa di 53 cifre binarie, ma visto che si utilizza la notazione normalizzata, solo 52 sono necessarie.

• Visto che si compone di 64 bits, di cui 53 per mantissa e segno, 11 possono essere dedicate all’esponente; quindi si possono rappresentare numeri interi da 0 a 2047; in questo ambito i numeri 0 e 2047 vengono utilizzati per significati speciali e quindi utilizzati esclusivamente gli interi N_{e} da 1 a 2046; fissato bias=1023 = 2^{10} − 1 per ottenere gli esponenti e li si codifica con N_{e} tale che

e = N_{e} − bias .

Evidentemente il piú piccolo esponente e_{min} = −(1 − 1023) = −1022, mentre il piú grande é e_{max} = 2046 − 1023 = 1023.

Nell’implementazione Matlab,

• Il minimo numero positivo in precisione doppia é

x_{min} = 2.225073858507201e − 308;

• Il massimo numero positivo in precisione doppia é

x_{max} = 1.797693134862316e + 308;

• la cardinalità dell’insieme F (2, 53, −1022, 1023) è ≈ 3.7038e + 19.

Nota. 1.15. Si può credere che i numeri utilizzati per le precisioni singole e doppie siano sempre stati questi, ma non è cosí. In altri sistemi, differentemente da IEEE754, il numero di cifre per esponente e mantissa sono risultati molteplici, ad esempio su

• Zuse Z1 (≈ 1936), la mantissa corrispondeva a 22 cifre (14 per la mantissa e 8 per l’esponente),

• IBM 3033 in precisione doppia si aveva, oltre alla base β = 16, una mantissa di 14 cifre,

• PRIME 850, che lavorava in base β = 2, si aveva in precisione singola 23 cifre di mantissa, mentre in precisione doppia 47 cifre di mantissa, [? , p.13].

Nota. 1.16. In generale la cardinalità dei numeri macchina è molto alta, spesso ordini di grandezza maggiore di quanto sia necessaria nelle scienze applicate [? ].

On a computer, the interval [1, 2], for example, is approximated by about 10^{16} numbers. It is interesting to compare the fineness of this discretization with that of the discretizations of physics. In a handful of solid or liquid or a balloonful of gas, the number of atoms or molecules in a line from one point to another is on the order of 108 (the cube root of Avogadro’s number).

Such a system behaves enough like a continuum to justify our definitions of physical quantities such as density, pressure, stress, strain, and temperature. Computer arithmetic, however, is more than a million times finer than this.

Another comparison with physics concerns the precision to which fundamental constants are known, such as (roughly) 4 digits for the gravitational constant G, 7 digits for Planck’s constant h and the elementary charge e, and 12 digits for the ratio μ_{e}/μ_{B} of the magnetic moment of the electron to the Bohr magneton. At present, almost nothing in physics is known to more than 12 or 13 digits of accuracy.

Thus IEEE numbers are orders of magnitude more precise than any number in science. (Of course, purely mathematical quantities like π are another matter.)

1.3.3. Troncamento e arrotondamento

Se x è il numero reale

x = sign(x)(0.d_{1}, … d_{t}, …)_{β}β^{e} = sign(x)β^{e} ∑_{k=1}^{+∞} d_{k}β^{−k}, d_{1} > 0

allora con fl(x) denotiamo

fl(x) = sign(x)(0.d_{1}, …, d_{t−1}, \overline{d_{t}})β^{e} = sign(x)β^{e}(∑_{k=1}^{t−1} d_{k}β^{−k} + \overline{d_{t}}β^{−t})

in cui

• se si effettua il troncamento

\overline{d_{t}} = d_{t};

• se viene effettuato l’arrotondamento.

\overline{d_{t}} = \{ \begin{matrix}d_{t}, d_{t+1} < β/2 \\ d_{t} + 1, d_{t+1} ≥ β/2 \end{matrix}

La scelta tra troncamento o arrotondamento è eseguita a priori dal sistema numerico utilizzato dal calcolatore. L’ultima, è la più diffusa.

Curiosa la descrizione in [? , p.5]

A machine register acts much like the infamous Procrustes bed in Greek mythology. Procrustes was the innkeeper whose inn had only beds of one size. If a fellow came along who was too tall to fit into his beds, he cut off his feet. If the fellow was too short, he stretched him. In the same way, if a real number comes along that is too long, its tail end (not the head) is cutoff; if it is too short, it is padded by zeros at the end.

Esempio. Si consideri

π = (+1) · (0.31415926535 . . . . . .)_{10} · 10^{1}.

Se lo tronchiamo alla quarta cifra decimale otteniamo

fl(π) = +1 · (0.3141)_{10} · 10^{1},

mentre se lo arrotondiamo

fl(π) = +1 · (0.3142)_{10} · 10^{1}.

Evidentemente, per effettuare un arrotondamento, il calcolatore deve immagazzinare in un extra-bit la cifra d_{t+1} (rifletterci un po’).

Se l’esponente e del numero x è

• minore di L si commette un’errore di underflow,

• maggiore di U si commette un’errore di overflow.

Se invece e ∈ [L, U], x = ±β^{e} ∑_{k=1}^{+∞} d_{k}β^{−k} e per qualche t^{∗} > t si ha d_{t^{∗}} > 0 allora si commette un errore di troncamento o arrotondamento a seconda si effettui l’uno o l’altro.

1.4. La precisione di macchina

Una delle costanti rilevanti all’interno del sistema floating point è la precisione di macchina.

Definizione. Si dice precisione di macchina, quel numero eps che rappresenta la distanza tra 1 ed il successivo numero in virgola mobile.

Calcoliamo tale valore. Necessariamente il numero successivo a 1 ha per mantissa d_{1} = 1, …, d_{t−1} = 0, d_{t} = 1 e esponente e = 1 e quindiè x = 1+β^{1}·β^{−t} = 1+β^{1−t}, per cui

eps = β^{1−t}.

Si sottolinea che la precisione di macchina eps non coincide in generale con il più piccolo numero positivo rappresentabile dal calcolatore, in quanto se t ≠ −L, come usualmente accade, allora

β^{1−t} ≠ β^{L−1}.

In Matlab ad esempio, x_{min} è molto inferiore di eps, in quanto

• il piú piccolo numero positivo in modulo è

x_{min} = 2^{−1022} = 2.225073858507201 · 10^{−308},

• la precisione di macchina è

eps = 2^{−52} = 2.220446049250313 · 10^{−16}.

Nota. 1.17. Usualmente

• un numero macchina il cui esponente è più piccolo di L produce il cosidetto underflow,

• un numero macchina il cui esponente è più grande di U produce il cosidetto overflow.

Osserviamo che esistono numeri denormalizzati, il cui range estende quello dei normalizzati. II piu piccolo numero positivo denormalizzato e 0.494·10^{−323}, permettendo di passare gradualmente dall’underflow allo 0.

1.5. Errori relativi e assoluti

Fissato un numero da approssimare

x^{∗} ∈ ℝ,

si definisce

1. errore assoluto tra x e x^{∗} la quantità

|x − x^{∗}|

2. errore relativo tra x e x^{∗} ≠ 0 la quantità

\frac{|x − x^{∗}|}{|x^{∗}|}.

Si noti che se |x^{∗}| = 0, cioè x^{∗} = 0, allora non ha senso la scrittura propria dell’errore relativo.

Nota. 1.18. Fissato un vettore da approssimare

x^{∗} = (x_{1}^{∗}, …, x_{n}^{∗}) ∈ ℝ^{n},

si definisce 1. errore assoluto tra x e x^{∗} la quantità

∥x − x^{∗}∥

dove

∥y∥ = \sqrt{∑_{i=1}^{n} y_{i}^{2}}, y = (y_{1}, …, y_{n});

2. errore relativo tra x e x^{∗} ≠ 0 la quantità

\frac{∥x − x^{∗}∥}{∥x^{∗}∥}.

Si noti che se ∥x^{∗}∥ = 0, cioè x^{∗} = 0 per la proprietà delle norme [? , p.481], allora non ha senso la scrittura propria dell’errore relativo. Inoltre se il vettore ha un solo componente, allora la norma ∥ · ∥ coincide con l’usuale valore assoluto | · |.

Sia

x = sign(x)β^{e} ∑_{k=1}^{+∞} d_{k}β^{−k} ≠ 0, d_{1} > 0

un numero reale e

fl(x) = sign(x)(0.d_{1}, …, d_{t−1}, \overline{d_{t}})β^{e} = sign(x)β^{e}(∑_{k=1}^{t−1} d_{k}β^{−k} + \overline{d_{t}}β^{−t}), d_{1} > 0

la sua approssimazione in virgola mobile.

L’errore relativo compiuto nell’approssimare x con fl(x), visto che |x| ≥ |d_{1}β^{e−1}| ≥ |β^{e−1}| poichè d_{1} > 0, verifica

\frac{|x − fl(x)|}{|x|} ≤ \frac{|sign(x)β^{e} ∑_{k=1}^{+∞} d_{k}β^{−k} − (sign(x)β^{e}(∑_{k=1}^{t−1} d_{k}β^{−k} + \overline{d_{t}}β^{−t}))|}{|β^{e−1}|}

= β · |∑_{k=t}^{+∞} d_{k}β^{−k} − \overline{d_{t}}β^{−t}|  (6)

Nel caso si effettui il troncamento abbiamo d_{t} = \overline{d_{t}} e quindi, essendo eps = β^{1−t} la precisione di macchina

β · |∑_{k=t}^{+∞} d_{k}β^{−k} − \overline{d_{t}}β^{−t}| = β · | ∑_{k=t+1}^{+∞} d_{k}β^{−k}| < β · (β · β^{−t−1}) = β^{1−t} = eps

mentre se si effettua l’arrotondamento abbiamo dopo alcuni calcoli che

\frac{|x − fl(x)|}{|x|} < \frac{β^{1−t}}{2} = \frac{eps}{2} .

La quantità \frac{eps}{2} è detta unità di arrotondamento.

Il minor errore relativo compiuto nell’approssimare x con fl(x) via arrotondamento, suggerisce perchè quest’ultimo sistema sia più utilizzato.

Nota. 1.19. Da questo si deduce che i numeri in virgola mobile non sono dunque un insieme denso, come sono invece i numeri reali, bensì un insieme discreto di valori sull’asse reale, in posizioni non equispaziate (se così fosse, visto che i numeri più grandi sono interi, pure eps lo sarebbe).

2. Operazioni con i numeri macchina

Per quanto concerne le quattro operazioni aritmetiche, si può dimostrare che la divisione e la moltiplicazione di tipo floating-point producono un errore relativo piccolo rispetto al risultato esatto, cosa che non succede sempre per somma e sottrazione.

Più precisamente se fl(x) è il numero macchina che corrisponde a x, denotiamo con

x ⊕ y := fl(fl(x) + fl(y))  (7)
x ⊖ y := fl(fl(x) + fl(y))  (8)
x ⊗ y := fl(fl(x) · fl(y))  (9)
x ⊘ y := fl(fl(x) : fl(y))  (10)

e per ⊙ una delle operazioni precedentemente introdotte, cioè una tra ⊕, ⊖, ⊗, ⊘ corrispondenti a op cioè una tra +, −, ·, :. Inoltre porremo

ϵ_{x} := \frac{|x − fl(x)|}{|x|}  (11)
ϵ_{x,y}^{⊙} := \frac{|(x opy) − (x ⊙ y)|}{|x opy|}  (12)

Nota. 2.1. Osserviamo, che a volte nei manuali si trova qualcosa di diverso. Ad esempio,

x ⊕ y = fl(x + y).

In realtà il sistema floating point, utilizzando ulteriori bits nei calcoli, garantisce che sia

x ⊕ y = fl(fl(x) + fl(y)) = fl(x + y).

In effetti, da [? ]

The IEEE standard requires that the result of addition, subtraction, multiplication and division be exactly rounded. That is, the result must be computed exactly and then rounded to the nearest floating-point number (using round to even).

The section Guard Digits pointed out that computing the exact difference or sum of two floating-point numbers can be very expensive when their exponents are substantially different. That section introduced guard digits, which provide a practical way of computing differences while guaranteeing that the relative error is small.

However, computing with a single guard digit will not always give the same answer as computing the exact result and then rounding. By introducing a second guard digit and a third sticky bit, differences can be computed at only a little more cost than with a single guard digit, but the result is the same as if the difference were computed exactly and then rounded . Thus the standard can be implemented efficiently.

Altre volte si trova equivalentemente [? ]

If ∗ is one of these four operations in its ideal form and ⊙ is the same operation as realized on the computer, then for any floating-point numbers x and y , assuming that there is no underflow or overflow,

x ⊙ y = (x ∗ y)(1 + ϵ).

Here ϵ is a very small quantity, no greater in absolute value than a number known as machine epsilon, denoted by ϵ_{mach}, that measures the accuracy of the computer. In the IEEE system, ϵ_{mach} = 2^{−53} ≈ 1.1 × 10^{−16}.

I calcoli che faremo successivamente per studiare le quattro operazioni in floatingpoint, possono essere eseguiti anche in questo modo.

2.1. Proprietá commutativa, associativa e distributiva delle operazioni floating point

Per prima cosa osserviamo che alcune proprietà caratteristiche di +, −, · e : non valgono per i corrispettivi ⊕, ⊖, ⊗, ⊘. Infatti per la somma ⊕ e il prodotto ⊗ vale la proprietà commutativa ma non vale quella associativa, e nemmeno quella distributiva.

Quindi l’ordine in cui sono eseguite le operazioni può variare il risultato. Inoltre esistono l’elemento neutro della moltiplicazione, divisione, addizione e

sottrazione, come pure l’opposto di ogni numero floating-point.

Esempio. Mostriamo come non valga la proprietà associativa della somma. Consideriamo F (10, 4, −10, +10) con fl basato sul troncamento. Sia a = 2000, b = 2.5, c = 7.8. Allora:

(a ⊕ b) ⊕ c = 0.2002 · 10^{4} ⊕ 0.7800 · 10^{1} = 0.2009 · 10^{4}
a ⊕ (b ⊕ c) = 0.2000 · 10^{4} ⊕ 0.1030 · 10^{2} = 0.2010 · 10^{4}

e quindi (a ⊕ b) ⊕ c ≠ a ⊕ (b ⊕ c).

Per quanto riguarda la proprietà distributiva, in F (10, 4, −10, +10), posti a = 0.800000 · 10^{9}, b = 0.500009 · 10^{4}, c = 0.500008 · 10^{4}

a ⊗ (b ⊖ c) = 0.800000 · 10^{9} ⊗ 0.1 · 10^{−2} = 0.800000 · 10^{9}
(a ⊗ b) ⊖ (a ⊗ c) = (0.800000 · 10^{9} ⊗ 0.500009 · 10^{4}) ⊖

⊖(0.800000 · 10^{9} ⊗ b = 0.500008 · 10^{4}) → overflow(13)

(l’overflow è dovuto alla sequenza di operazioni da compiere).

2.2. Errori nelle operazioni e loro propagazione

Consideriamo ora come si propagano gli errori in ⊕, ⊖, ⊗, ⊘. Si mostra che (cf. [? , p.78])

ϵ_{x,y}^{⊕} ≤ |\frac{x}{x + y} | ϵ_{x} + |\frac{y}{x + y} | ϵ_{y}, x + y ≠ 0  (14)
ϵ_{x,y}^{⊖} ≤ |\frac{x}{x − y} | ϵ_{x} + |\frac{y}{x − y} | ϵ_{y}, x.y ≠ 0  (15)
ϵ_{x,y}^{⊗} ⪅ ϵ_{x} + ϵ_{y}  (16)
ϵ_{x,y}^{⊘} ≤ |ϵ_{x} − ϵ_{y}|  (17)

dove per ogni operazione floating point ⊙ relativamente all’operazione op (somma, sottrazione, moltiplicazione o divisione)

ϵ_{x} = |x − fl(x)|/|x|,
ϵ_{y} = |y − fl(y)|/|y|,
ϵ_{x,y}^{⊙} = |fl(x op y) − (x op y)|/|x op y|

Proviamo di seguito sotto opportune ipotesi il caso della somma, sottrazione e della moltiplicazione, lasciando al lettore quello della divisione. Con a ⪅ b si intende a inferiore o uguale circa a b.

Proposizione 2.2. Se,

• x + y ≠ 0, x ≠ 0, y ≠ 0,

• fl(fl(x) + fl(y)) = fl(x) + fl(y),

• ϵ_{x} = |x − fl(x)|/|x|,

• ϵ_{y} = |y − fl(y)|/|y|,

• ϵ_{x,y}^{⊕} = |fl(x ⊕ y) − (x + y)|/|x + y|

allora

ϵ_{x,y}^{⊕} ≤ |\frac{x}{x + y} | ϵ_{x} + |\frac{y}{x + y} | ϵ_{y}.

Dimostrazione. Nelle ipotesi precedenti abbiamo

ϵ_{x,y}^{⊕} = \frac{|(x + y) − (x ⊕ y)|}{|x + y|} = \frac{|(x + y) − fl(fl(x) + fl(y))|}{|x + y|} = \frac{|(x − fl(x)) + (y − fl(y))|}{|x + y|}
≤ \frac{|x|}{|x + y|} · \frac{|x − fl(x)|}{|x|} + \frac{|y|}{|x + y|} · \frac{|y − fl(y)|}{|y|} = |\frac{x}{x + y} | ϵ_{x} + |\frac{y}{x + y} | ϵ_{y}.  (18)

Il punto chiave è che per x + y ≈ 0 abbiamo che per piccoli errori sui dati ϵ_{x}, ϵ_{y} possono aversi errori rilevanti in ϵ_{x,y}^{⊕} (tale fenomeno è noto come cancellazione).

Proposizione 2.3. Se,

• x − y ≠ 0, x ≠ 0, y ≠ 0,

• fl(fl(x) − fl(y)) = fl(x) − fl(y)

• ϵ_{x} = |x − fl(x)|/|x|,

• ϵ_{y} = |y − fl(y)|/|y|,

• ϵ_{x,y}^{⊖} = |fl(x ⊖ y) − (x + y)|/|x + y|

allora

ϵ_{x,y}^{⊖} ≤ |\frac{x}{x − y} | ϵ_{x} + |\frac{y}{x − y} | ϵ_{y}.

Dimostrazione. Nelle ipotesi precedenti abbiamo

ϵ_{x,y}^{⊖} = \frac{|(x − y) − (x ⊖ y)|}{|x − y|} = \frac{|(x − y) − fl(fl(x) − fl(y))|}{|x − y|} = \frac{|(x − fl(x)) − (y − fl(y))|}{|x − y|}
≤ \frac{|x|}{|x − y|} · \frac{|x − fl(x)|}{|x|} + \frac{|y|}{|x − y|} · \frac{|y − fl(y)|}{|y|} = |\frac{x}{x − y} | ϵ_{x} + |\frac{y}{x − y} | ϵ_{y}  (19)

Nuovamente, per x − y ≈ 0 abbiamo che per piccoli errori sui dati ϵ_{x}, ϵ_{y} possono aversi errori rilevanti in ϵ_{x,y}^{⊖}. Qualora ciò accada, si parla nuovamente di fenomeno di cancellazione.

Cancellazione Consideriamo il sistema F (10, 5, −5, +5) (con arrotondamento) e i due numeri

a = 0.73415776, b = 0.73402350.

Naturalmente,

fl(a) = 0.73416, fl(b) = 0.73402

e quindi

a − b = 0.00013426, a ⊖ b = 0.00014 = 10^{−3} · 0.14000

con un errore relativo pari a

|(a − b) − (a ⊖ b)|/|a − b| ≈ 4.3/100

cioè molto grande (oltre il 4 per cento!).

Proposizione 2.4. Supposti,

• x · y ≠ 0, x ≠ 0, y ≠ 0, |fl(y)|/y ≈ 1,

• fl(fl(x) · fl(y)) = fl(x) · fl(y) (cioè fl(x) · fl(y) numero macchina),

posto ϵ_{x} = \frac{|x−fl(x)|}{|x|}ed ϵ_{y} = \frac{|y−fl(y)|}{|y|}, ϵ_{x,y}^{⊗} = |x ⊗ y − x · y|/|x · y| allora

ϵ_{x,y}^{⊗} ⪅ ϵ_{x} + ϵ_{y}.

Dimostrazione. Nelle ipotesi precedenti abbiamo

ϵ_{x,y}^{⊗} = \frac{|x · y − x ⊗ y|}{|x · y|} = \frac{|x · y − fl(fl(x) · fl(y))|}{|x · y|} = \frac{|x · y − fl(x) · fl(y)|}{|x · y|}
= \frac{|x · y − x · fl(y) + x · fl(y) − fl(x) · fl(y)|}{|x · y|}
≤ \frac{|x · (y − fl(y))| + |fl(y) · (x − fl(x))|}{|x · y|}
= \frac{|y − fl(y)|}{|y|} + \frac{|fl(y) · (x − fl(x))|}{|x · y|}
≈ ϵ_{x} + ϵ_{y}  (20)

dove si ricorda che

\frac{|fl(y)|}{|x · y|} ≈ \frac{1}{|x|}.

Per quanto riguarda la divisione, forniamo una traccia: basta studiare il caso 1/y e poi ricordare i risultati del prodotto.

Nota. 2.5. Per stabilità di una generica operazione macchina op (si pensi ad esempio ∗ come alla somma oppure al prodotto), si intende che se ⊙ è l’omologa versione floating-point allora esiste C > 0 tale che

ϵ_{x,y}^{⊙} ≤ C(ϵ_{y} op ϵ_{y})

per ogni x, y per cui la scrittura abbia significato.

Ne consegue che la somma e la sottrazione non sono stabili mentre lo sono il prodotto e la divisione [? , p.10].

Nota. 2.6 (Difficile, necessario aver seguito il corso di Analisi II). Osserviamo che dalla formula di Taylor in due variabili, posto per brevità di notazione f′_{x} = ∂f/∂x, f′_{y} = ∂f/∂y,

f(x, y) ≈ f(\overline{x}, \overline{y}) + f′_{x}(\overline{x}, \overline{y})(x − \overline{x}) + f′_{y}(\overline{x}, \overline{y})(y − \overline{y})  (21)

necessariamente, se x, y, f(x, y) ≠ 0 allora posti ϵ_{x} = |x − \overline{x}|/|x|, ϵ_{y} = |y − \overline{y}|/|y| gli errori relativi sui dati

\frac{|f(x, y) − f(\overline{x}, \overline{y})|}{|f(x, y)|} ⪅ \frac{|f′_{x}(\overline{x}, \overline{y})||x − \overline{x}||x|}{|f(x, y)||x|} + \frac{|f′_{y}(\overline{x}, \overline{y})||y − \overline{y}||y|}{|f(x, y)||y|}
= \frac{|f′_{x}(\overline{x}, \overline{y})||x|}{|f(x, y)|} · ϵ_{x} + \frac{|f′_{y}(\overline{x}, \overline{y})||y|}{|f(x, y)|} · ϵ_{y}  (22)

Da (??), si possono ottenere risultati simili ai precedenti per ⊕, ⊗. Uno dei vantaggi è che si possono studiare operazioni più generali come ad esempio cosa succeda se si esegue exp(x + y).

3. Alcune problematiche numeriche

Si consideri un problema scientifico che abbia a che fare con una procedura di calcolo. Sussitono varie tipologie di errore nei calcoli:

• errori di modelizzazione: vengono utilizzate equazioni per rappresentare un evento fisico, e il modello matematico è solo una semplificazione di quello reale;

• errori di tabulazione: vengono immessi dal programmatore o dal programma dati errati (cosa che ad esempio era comune una volta quando bisognava tabulare manualmente numeri con molte cifre); un errore di questo tipo ha causato nel 1996 la perdita del razzo Ariane 5 per una cattiva conversione di un numero;

Il 4 giugno 1996, è avvenuta un esplosione a soli 40 secondi dal lancio dovuta a un overflow per una conversione di un numero reale memorizzato usando 64 bit in un numero intero con soli 16 bit. Il numero che rappresentava la velocitá orizzontale del razzo era piú grande del massimo numero rappresentabile con soli 16 bit. La missione era costata 7.5 miliardi di dollari.

Errori di questo tipo hanno causato pure il crollo di una piattaforma petrolifera nel mare del Nord (Norvegia) nel 1991 e la distruzione del veicolo spaziale Mars Climate orbiter nel 1999 (cf. [? ], [? ]).

• errori dovuti a misure fisiche: sono dovuti a una imprecisa osservazione sperimentale, spesso dovuta agli strumenti utilizzati, ad esempio la velocitá della luce nel vuoto è

c = (2.997925 + ϵ) · 10^{10}cm/s, |ϵ| ≤ 0.000003

e quindi nei calcoli che hanno a che vedere con c necessariamente si compie un errore; • errori di rappresentazione e di aritmetica di macchina: sono errori dovuti all’ar-

rotondamento o troncamento di un numero, e sono inevitabili quando si usa l’aritmetica floating-point; sono la sorgente principale di errore di alcuni problemi come ad esempio la soluzione di sistemi lineari.

Di seguito vediamo alcune problematiche dovute agli errori di rappresentazione e di aritmetica di macchina. Per una buona serie di esempi si veda [? , p.46].

4. Valutazione di una funzione

Nel valutare una funzione continua f : Ω ⊆ ℝ → ℝ in un punto x ∈ Ω bisogna tener conto che il dato x può essere per varie ragioni affetto da errori, tipicamente di misura o di rappresentazione al calcolatore, così da essere approssimato con un certo \overline{x}. Di conseguenza, invece di f(x) si valuta f(\overline{x}) e ci si domanda se ci siano considerevoli differenze.

Una funzione f risulta difficile da valutare al calcolatore nel punto x ≠ 0 in cui f(x) ≠ 0 qualora a piccoli errori (relativi) sui dati

\frac{|x − x_{c}|}{|x|}

corrispondano grandi errori (relativi) sui risultati

\frac{|f(x) − f(x_{c})|}{|f(x)|}.

Quindi è importante discutere la quantità

K(f, x, x_{c}) = \frac{|f(x) − f(x_{c})|/|f(x)|}{|x − x_{c}|/|x|}.

La valutazione di f in un punto x si dice bencondizionata se a piccole variazioni relative sui dati corrispondono piccole variazioni sui risultati, malcondizionata se ciò non avviene.

Se f è derivabile con continuità nel più piccolo intervallo I aperto contenente x ed x_{c}, per il teorema della media, essendo

f(x) − f(x_{c}) = f′(ξ) · (x − x_{c}), ξ ∈ I,

e quindi da |f(x) − f(x_{c})|/|x − x_{c}| = f′(ξ) ≈ f′(x) ricaviamo facilmente che

K(f, x, x_{c}) ≈ K(f, x) := \frac{|x · f′(x)|}{|f(x)|}.

La quantità K(f, x) si chiama condizionamento di f nel punto x. Più piccola risulta K(f, x) e meno la funzione amplifica un piccolo errore sul dato x. Naturalmente la funzione può essere bencondizionata su un certo dato x_{1} ∈ Ω e malcondizionata su un certo altro x_{2} ∈ Ω.

Osserviamo che il condizionamento di una funzione non dipende dall’algoritmo con cui viene valutata, ma è inerente alla funzione stessa f e al dato x.