risultati inattesi con specifica modalità di arrotondamento

Question

Prima di tutto, mi spiace per il titolo cattivo:/

sto cercando di riprodurre i risultati di una carta sul calcolo gli autovalori di una matrice simmetrica tridiagonale. Sto determinando alcuni valori 'limiti superiori e inferiori' arrotondando rispettivamente a più e meno infinito.

Invece di modificare la modalità di arrotondamento ogni volta, utilizzo solo il "trucco": fl⁻ (y) = -fl⁺ (-y), dove fl⁻ (y) è il valore di y quando si utilizza il segno meno la modalità di arrotondamento infinito e fl⁺ (y) è il valore di y quando si utilizza la modalità di arrotondamento a più infinito. Così, ho il seguente pezzo di codice in C:

fesetround(FE_UPWARD); 
first = - (-d[i] + x); 
second = (- ((e[i-1]*e[i-1])/a_inf)); 
a_inf = first + second; 

first = d[i] - x; 
second = - ((e[i-1]*e[i-1])/a_sup); 
a_sup = first + second; 

e funziona bene tranne che per un esempio in cui A_INF mi dà il risultato giusto, ma a_sup dà il risultato sbagliato, anche se entrambe le prime e seconde variabili sembrano avere gli stessi valori.

Tuttavia, se faccio così:

fesetround(FE_UPWARD); 
    first = - (-d[i] + x); 
    second = (- ((e[i-1]*e[i-1])/a_inf)); 

    fesetround(FE_DOWNWARD); 
    first = - (-d[i] + x); 
    second = (- ((e[i-1]*e[i-1])/a_sup)); 

ottengo i risultati giusti. Quindi, se uso il trucco fl⁻ (y) = -fl⁺ (-y), ottengo i risultati corretti, se cambio la modalità di arrotondamento e utilizzo l'espressione originale ottengo risultati errati. Qualche idea del perché?

In entrambi i casi, le variabili primi e secondi valori sono i seguenti:

first 1.031250000000000e+07, second -1.031250000000000e+07 
first 1.031250000000000e+07, second -1.031250000000000e+07 

E i valori corretti per A_INF e a_sup sono -1.862645149230957e-09 e + 1.862645149230957e-09, rispettivamente, ma in il primo caso a_sup = 0, che è sbagliato

Quello che sto cercando di indovinare che sta accadendo è una sorta di cancellazione catastrofica, ma non ho idea su come risolverlo in questo caso ...

Grazie a avanzare!

Answer 1

Il primo problema riscontrato è che si sta utilizzando solo il "trucco" per ottenere l'arrotondamento di first verso -inf, e non second, quindi è ancora arrotondato verso + inf.

Il secondo problema è che C non offre alcun tipo di garanzia su come esegue calcoli in virgola mobile. In particolare, il compilatore è libero di riordinare e riassociare le operazioni come ritiene opportuno per le prestazioni, anche se tali riarrangiamenti possono modificare il comportamento di arrotondamento del programma. Quindi, quando si dice:

first = - (-d[i] + x); 

il compilatore può riorganizzare quello e fare un unico sottrarre, invece di negare, aggiungere, negare, che inverte la direzione di arrotondamenti. Ora a volte puoi far funzionare le cose come ti aspetti disabilitando tutte le ottimizzazioni, ma anche con quelle non c'è garanzia.