Floating Point Error, quando i computer commettono errori di calcolo.

Floating Point Error, quando i computer commettono errori di calcolo

Comprendere il Floating Point error

L’errore di floating point origina dalle limitazioni nella rappresentazione dei numeri reali nella memoria del computer, portando a imprecisioni nei calcoli che possono accumularsi e influenzare i risultati.

La maggior parte delle frazioni decimali non può essere rappresentata con precisione come frazioni binarie.

Rappresentazione binaria:

Molti numeri decimali non possono essere rappresentati esattamente in binario. Ad esempio, 0,1 non può essere espresso con precisione in forma binaria, causando piccoli errori di rappresentazione. La rappresentazione binaria più vicina potrebbe essere simile a 0,1000000000000000055511151231257827021181583404541015625.

Limiti di precisione:

La doppia precisione fornisce circa 15-17 cifre decimali significative. Se i calcoli danno come risultato numeri che richiedono una precisione maggiore, si verifica un arrotondamento che può portare a imprecisioni.

Accumulo di errori:

Nei calcoli che coinvolgono più operazioni in virgola mobile, possono accumularsi piccoli errori. Ad esempio, sottraendo due numeri quasi uguali si può ottenere un errore sostanziale a causa della perdita di precisione.

Arrotondamento:

Quando il risultato di un calcolo deve essere memorizzato come doppio e non può essere rappresentato esattamente, si verifica un arrotondamento. Questo arrotondamento porta a discrepanze. Ad esempio, la somma di diversi numeri piccoli potrebbe non produrre il risultato esatto previsto, ma un arrotondamento molto vicino.

Di seguito è possibile osservare come attraverso la scrittura di un semplice codice Delphi, utilizzando numeri di tipo Real di varie dimensioni e compilato su architetture diverse (32 e 64 bit), possa produrre alcuni errori evidenti. Per dimostrare l’entita di questi errori sono stati utilizzati tutti i tipi numerici in virgola mobile sia a precisione singola che precisione doppia.

Nota: il tipo real48 non è nativo dell’architettura del processore Intel. In Delphi è stato mantenuto per retrocompatibilità ed è alias del tipo Real.

Codice Delphi, compilato a 32 bit su Sistema Operativo Windows 10 x64 eseguito su processori Intel Xeon E5 2678 v3 (64 bit).

Real48 =>    2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 2,16840434497101E-19
Real48 =>    0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Real48 =>   -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = -2,71050543121376E-20

Single =>    2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 2,16840434497101E-19
Single =>    0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Single =>   -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = -2,71050543121376E-20

Double =>    2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 2,16840434497101E-19
Double =>    0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Double =>   -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = -2,71050543121376E-20

Real =>      2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 2,16840434497101E-19
Real =>      0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Real =>     -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = -2,71050543121376E-20

Extended =>  2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 2,16840434497101E-19
Extended =>  0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Extended => -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = -2,71050543121376E-20

Comp =>      2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 0
Comp =>      0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Comp =>     -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = 0

Currency =>  2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 0
Currency =>  0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Currency => -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = 0

Real48 =>    6.92736 - 6.92735 = 9,99999999999612E-6

Single =>    6.92736 - 6.92735 = 9,99999974737875E-6

Double =>    6.92736 - 6.92735 = 1E-5

Real =>      6.92736 - 6.92735 = 1E-5

Extended =>  6.92736 - 6.92735 = 1E-5

Comp =>      6.92736 - 6.92735 = 0

Currency =>  6.92736 - 6.92735 = 0

Double =>                             1.01^(365) = 37,7834343328873
Double =>    1.01 * 1.01 * 1.01 * 1.01 * 1.01... = 37,783434332523

Double =>                             0.99^(365) = 0,0255179644522911
Double =>    0.99 * 0.99 * 0.99 * 0.99 * 0.99... = 0,0255179644522912

Codice Delphi, compilato a 64 bit su Sistema Operativo Windows 10 x64 eseguito su processori Intel Xeon E5 2678 v3 (64 bit).

Real48 =>    2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 8,88178419700125E-16
Real48 =>    0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Real48 =>   -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = -5,55111512312578E-17

Single =>    2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 8,88178419700125E-16
Single =>    0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Single =>   -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = -5,55111512312578E-17

Double =>    2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 8,88178419700125E-16
Double =>    0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Double =>   -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = -5,55111512312578E-17

Real =>      2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 8,88178419700125E-16
Real =>      0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Real =>     -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = -5,55111512312578E-17

Extended =>  2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 8,88178419700125E-16
Extended =>  0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Extended => -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = -5,55111512312578E-17

Comp =>      2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 0
Comp =>      0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Comp =>     -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = 0

Currency =>  2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 - 3 = 0
Currency =>  0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 - 3 = 0
Currency => -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = 0

Real48 =>    6.92736 - 6.92735 = 1,00000000005096E-5

Single =>    6.92736 - 6.92735 = 9,99999974737875E-6

Double =>    6.92736 - 6.92735 = 1,00000000005096E-5

Real =>      6.92736 - 6.92735 = 1,00000000005096E-5

Extended =>  6.92736 - 6.92735 = 1,00000000005096E-5

Comp =>      6.92736 - 6.92735 = 0

Currency =>  6.92736 - 6.92735 = 0

Double =>                             1.01^(365) = 37,7834343328873
Double =>    1.01 * 1.01 * 1.01 * 1.01 * 1.01... = 37,783434332523

Double =>                             0.99^(365) = 0,025517964452291
Double =>    0.99 * 0.99 * 0.99 * 0.99 * 0.99... = 0,025517964452348

Gli esempi sopra, tre tipi di operazioni aritmentiche;

Addizione, (real 32 bit) :

• 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 – 3 = 2,16840434497101E-19 oppure 0,000000000000000000216840434497101
• 0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 – 3 = 0
• -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = -2,71050543121376E-20 oppure -0,0000000000000000000271050543121376

Addizione, (real 64 bit) :

• 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 – 3 = 8,88178419700125E-16 oppure 0,000000000000000888178419700125
• 0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 2 – 3 = 0
• -3 + 2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 = -5,55111512312578E-17 oppure -0,0000000000000000555111512312578

Sottrazione, (real 32) :

• 6.92736 – 6.92735 = 1E-5 oppure 0,00001

Sottrazione, (real 64 bit) :

• 6.92736 – 6.92735 = 1,00000000005096E-5 oppure 0,0000100000000005096

Elevamento a potenza, (double 32 bit) :

• 1.01^(365) = 37,7834343328873
• 0.99^(365) = 0,0255179644522911

Elevamento a potenza, (double 64 bit) :

• 1.01^(365) = 37,7834343328873
• 0.99^(365) = 0,025517964452291

Moltiplicazione ricorsiva, (double 32 bit) :

• 1.01 * 1.01 * 1.01 * 1.01 * 1.01… = 37,783434332523
• 0.99 * 0.99 * 0.99 * 0.99 * 0.99… = 0,0255179644522912

Moltiplicazione ricorsiva, (double 64 bit) :

• 1.01 * 1.01 * 1.01 * 1.01 * 1.01… = 37,783434332523
• 0.99 * 0.99 * 0.99 * 0.99 * 0.99… = 0,025517964452348

Piccoli errori.

Appare ovvio che l’ entità di errore è molto piccola, sicuramente non in grado di influenzare la maggior parte di calcoli comuni, ma nel caso di calcolo scientifico, ed esempio chimica, fisica, astronomia sicuramente questo tipo di errore può avere impatto notevole.

Metodi per ridurre l’errore.

• Utilizzare tipi di dati con maggiore precisione
  • Doppia precisione (64 bit):
    utilizzare la doppia precisione invece della precisione singola per aumentare l’accuratezza. Ciò consente di ottenere cifre più significative e un intervallo più ampio, riducendo gli errori di arrotondamento.
  • Precisione estesa:
    alcuni sistemi supportano tipi con precisione ancora maggiore (ad esempio, long double in C o librerie specializzate), che possono fornire una maggiore accuratezza per i calcoli che coinvolgono numeri grandi.

• Algoritmi migliori
  • Somma di Kahan:
    utilizzare algoritmi progettati per ridurre al minimo l’accumulo di errori, come la Somma di Kahan, che compensa la perdita di precisione durante l’addizione.
  • Evitare la cancellazione catastrofica:
    riorganizzare i calcoli per evitare di sottrarre numeri quasi uguali, che possono portare a una significativa perdita di precisione.

• Librerie numeriche
  • Utilizzare librerie specializzate:
    Sfruttare librerie che puntano sull’alta precisione e sulla stabilità numerica, come la libreria GNU MPFR (Multiple Precision Floating-Point Reliable) per C o la libreria Arbitrary Precision Arithmetic.

• Calcolo simbolico:
  • Valutare l’utilizzo di librerie matematiche simboliche che manipolano le espressioni algebricamente anziché numericamente, mantenendo l’esattezza fino a quando non è necessaria un’approssimazione numerica.

• Ordine accurato delle operazioni:
  • Pianificare l’ordine delle operazioni per ridurre gli errori di arrotondamento, specialmente nelle serie di addizioni e moltiplicazioni. Raggruppare i termini in modo strategico per mantenere la precisione.

Calcoli Astronomici (meccanica orbitale)

Ciò che ha ispirato la scrittura di questo articolo è proprio il calcolo astronomico, ridurre gli errori di virgola mobile nei calcoli della meccanica orbitale è fondamentale a causa dell’ampia gamma di valori e della precisione richiesta.

Ecco una delle soluzioni per ridurre al minimo gli errori in questo tipo di calcolo tramite l’utilizzo della Libreria MPFR:

Esempio: Calcolo dell’Anomalia Media e Anomalia Vera del sole:

Il calcolo dell’Anomalia Media e Anomalia Vera del Sole richiede l’uso dei principi della meccanica orbitale, in particolare della formula:

Mean Anomaly = E−e⋅sin(E)
True Anomaly = 2⋅tan−1(1−e1+e​​⋅tan(2M​))

Dove:

M è l’anomalia media.
E è l’anomalia eccentrica.
e è l’eccentricità dell’orbita.

Supponiamo di avere:

Eccentricità e=0,0167
Anomalia eccentrica E=1,0 radianti.

Calcolo dell’Anomalia Media e Anomalia Vera utilizzando la libreria standard C

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

double calculate_mean_anomaly(double e, double E) {
    return E - e * sin(E);
}

double calculate_true_anomaly(double e, double M) {
    double sqrt_term = sqrt((1 + e) / (1 - e));
    return 2 * atan(sqrt_term * tan(M / 2.0));
}

int main() {
    double e = 0.0167; // Eccentricity
    double E = 1.0;    // Eccentric anomaly

    // Calculate Mean Anomaly
    double M = calculate_mean_anomaly(e, E);
    // Calculate True Anomaly
    double true_anomaly = calculate_true_anomaly(e, M);

    // Output results
    printf("Mean Anomaly (Standard C): %.15f\n", M);
    printf("True Anomaly (Standard C): %.15f\n", true_anomaly);

    return 0;
}

Mean Anomaly (Standard C): 0.996989144563351
True Anomaly (Standard C): 1.199230506769381

Calcolo dell’Anomalia Media e Anomalia Vera utilizzando la libreria standard mpfr:

#include <stdio.h>
#include <mpfr.h>

void calculate_mean_anomaly(mpfr_t M, mpfr_t e, mpfr_t E) {
    mpfr_t sin_E;
    mpfr_init2(sin_E, 256);
    mpfr_sin(sin_E, E, MPFR_RNDN);
    mpfr_mul(M, e, sin_E, MPFR_RNDN);
    mpfr_sub(M, E, M, MPFR_RNDN);
    mpfr_clear(sin_E);
}

void calculate_true_anomaly(mpfr_t true_anomaly, mpfr_t e, mpfr_t M) {
    mpfr_t sqrt_term, tan_half_M;
    mpfr_init2(sqrt_term, 256);
    mpfr_init2(tan_half_M, 256);

    mpfr_add_d(sqrt_term, e, 1.0, MPFR_RNDN); // 1 + e
    mpfr_sub_d(tan_half_M, 1.0, e, MPFR_RNDN); // 1 - e
    mpfr_sqrt(tan_half_M, tan_half_M, MPFR_RNDN); // sqrt(1 - e)
    mpfr_div(sqrt_term, sqrt_term, tan_half_M, MPFR_RNDN); // (1 + e)/(1 - e)
    mpfr_tan(tan_half_M, M, MPFR_RNDN); // tan(M/2)
    mpfr_mul(tan_half_M, tan_half_M, sqrt_term, MPFR_RNDN); // Final computation
    mpfr_atan(true_anomaly, tan_half_M, MPFR_RNDN); // atan(...)

    mpfr_mul_d(true_anomaly, true_anomaly, 2.0, MPFR_RNDN); // Multiply by 2
    mpfr_clear(sqrt_term);
    mpfr_clear(tan_half_M);
}

int main() {
    mpfr_t e, E, M, true_anomaly;
    mpfr_init2(e, 256);
    mpfr_init2(E, 256);
    mpfr_init2(M, 256);
    mpfr_init2(true_anomaly, 256);

    // Set values
    mpfr_set_d(e, 0.0167, MPFR_RNDN); // Eccentricity
    mpfr_set_d(E, 1.0, MPFR_RNDN);    // Eccentric anomaly

    // Calculate Mean Anomaly
    calculate_mean_anomaly(M, e, E);
    // Calculate True Anomaly
    calculate_true_anomaly(true_anomaly, e, M);

    // Output results
    mpfr_printf("Mean Anomaly (MPFR): %.15Ff\n", M);
    mpfr_printf("True Anomaly (MPFR): %.15Ff\n", true_anomaly);

    // Clear MPFR variables
    mpfr_clear(e);
    mpfr_clear(E);
    mpfr_clear(M);
    mpfr_clear(true_anomaly);

    return 0;
}

Mean Anomaly (MPFR): 0.9969891445633509000000000
True Anomaly (MPFR): 1.199230506769381183000000000

Standard e Pratiche

Lo standard IEEE 754 regola l’aritmetica in virgola mobile, garantendo un comportamento coerente su diversi sistemi informatici. Una profonda comprensione di questi principi è fondamentale per gli sviluppatori al fine di gestire e mitigare efficacemente gli errori di punto mobile. La maggior parte dei processori moderni aderisce allo standard IEEE 754 per l’aritmetica in virgola mobile. Questo standard stabilisce come vengono rappresentati i valori in virgola mobile e come devono essere eseguite le operazioni aritmetiche, ma non ne elimina gli errori.