Incoerenza di calcolo sinusoidale in Visual C++ 2012?

Question

Si consideri il seguente codice:

// Filename fputest.cpp 

#include <cmath> 
#include <cstdio> 

int main() 
{ 
    double x; 
    *(__int64 *) &x = 0xc01448ec3aaa278di64; // -5.0712136427263319 
    double sine1 = sin(x); 
    printf("%016llX\n", sine1); 
    double sine2; 
    __asm { 
    fld x 
    fsin 
    fstp sine2 
    } 
    printf("%016llX\n", sine2); 
    return 0; 
} 

Quando compilato con Visual C++ 2012 (cl fputest.cpp) e il programma viene eseguito, l'uscita è la seguente:

3FEDF640D8D36174 
3FEDF640D8D36175 

Domande:

• Perché questi due valori sono diversi?
• È possibile emettere alcune opzioni del compilatore in modo che i valori sinusoidali calcolati siano esattamente gli stessi?

Answer 1

Questo problema non è causato dalla conversione da lunga doppia a doppia. Potrebbe essere dovuto a inaccuratezza nella routine sin nella libreria matematica.

L'istruzione fsin viene specificata per produrre un risultato entro 1 ULP (nel formato doppio lungo) per gli operandi all'interno del suo intervallo (in base a Intel 64 e IA-32 Architectures Software Developer's Manual, ottobre 2011, volume 1, 8.3.10), in modalità round to to nearest. In un Intel Core i7, fsin del valore dell'interrogante, -5.07121364272633190495298549649305641651153564453125 o -0x1.448ec3aaa278dp + 2, produce 0xe.fb206c69b0ba402p-4. Possiamo facilmente vedere da questo esadecimale che gli ultimi 11 bit sono 100 0000 0010. Questi sono i bit che verranno arrotondati durante la conversione dal doppio lungo. Se sono superiori a 100 0000 0000, il numero verrà arrotondato per eccesso. Sono più grandi Pertanto, il risultato della conversione di questo doppio valore lungo in double è 0xe.fb206c69b0ba8p-4, che equivale a 0x1.df640d8d36175p-1 e 0.93631021832247418590355891865328885614871978759765625. Si noti inoltre che, anche se il risultato fosse inferiore di un ULP, gli ultimi 11 bit sarebbero comunque superiori a 100 0000 0000 e andrebbero comunque arrotondati. Pertanto questo risultato non dovrebbe variare su CPU Intel conformi alla documentazione di cui sopra.

Confronta questo per calcolare direttamente un seno a doppia precisione, utilizzando una routine ideale che produce risultati correttamente arrotondati. Il seno del valore è approssimativamente 0.93631021832247413051857150785044253634581268961333520518023697738674775240815140702992025520721336793516756640679315765619707343171517531053811196321335899848286682535203710849065933755262347468763562 (calcolato con Maple 10). Il doppio più vicino a questo è 0x1.df640d8d36175p-1. È lo stesso valore che abbiamo ottenuto convertendo il risultato fsin in double.

Pertanto, la discrepanza non è causata dalla conversione del lungo da doppio a doppio; la conversione del doppio risultato lungo fsin in raddoppio produce esattamente lo stesso risultato di una routine sin a precisione doppia ideale.

Non abbiamo una specifica per la precisione della routine sin utilizzata dal pacchetto Visual Studio del questionario. Nelle biblioteche commerciali è comune la possibilità di errori di 1 ULP o diversi ULP. Osservare quanto vicino il seno è in un punto in cui il valore di precisione doppia è arrotondato: È .498864 ULP (doppia precisione ULP) lontano da un doppio, quindi è 0,100366 ULP lontano dal punto in cui l'arrotondamento cambia. Pertanto, anche una minima inesattezza nella routine sin causerà il ritorno 0x1.df640d8d36174p-1 anziché il più vicino 0x1.df640d8d36175p-1.

Pertanto, suppongo che l'origine della discrepanza sia un'imprecisione molto piccola nella routine sin.

Answer 2

(Nota:..! Come accennato nei commenti, questo non funziona sul VC2012 che ho lasciato qui per informazioni di carattere generale non mi consiglia fare affidamento su tutto ciò che dipende dal livello di ottimizzazione comunque)

Non ho VS2012, ma sul compilatore VS2010 è possibile specificare /fp:fast sulla riga di comando e quindi ottengo gli stessi risultati. Questo fa sì che il compilatore generi codice "veloce" che non segue necessariamente l'arrotondamento e le regole richiesti in C++, ma che corrisponde al calcolo del linguaggio assembly.

Non riesco a provare questo in VS2012 ma immagino che abbia la stessa opzione.

Questo sembra funzionare in una build ottimizzata anche con /Ox come opzione.

Answer 3

Vedi Why is cos(x) != cos(y) even though x == y?

Come David mentioned in the comment, la discrepanza viene da spostare i dati in un registro FP una posizione di memoria (registro/ram) di dimensioni diverse. E non è sempre un compito neanche; anche un'altra operazione in virgola mobile vicina può essere sufficiente per svuotare un registro FP, rendendo inutile qualsiasi tentativo di garantire un particolare valore. Se avete bisogno di fare un paragone, si potrebbe essere in grado di mitigare alcune di queste forzando tutti i risultati in una posizione di memoria come segue:

float F1 = sin(a); float F2 = sin(b); if (F1 == F2) 

Tuttavia, anche questo non potrebbe funzionare. L'approccio migliore consiste nell'accettare che qualsiasi operazione in virgola mobile sarà solo "per lo più accurata" e che dal punto di vista dei programmatori questo errore sarà effettivamente imprevedibile e casuale, anche se la stessa operazione viene eseguita ripetutamente. Invece di

if (F1 == F2) 

si dovrebbe usare qualcosa per l'effetto di

if (isArbitrarilyClose(F1, F2)) 

o

if (absf(F1 - F2) <= n) 

dove n è un numero molto piccolo.

Answer 4

Il generatore di codice in VS2012 è stato sostanzialmente modificato per supportare automatic vectorization. Parte di questo cambiamento è che la matematica in virgola mobile x86 viene ora eseguita in SSE2 e non utilizza più la FPU, necessaria perché il codice FPU non può essere vettorizzato. SSE2 calcola con precisione a 64 bit invece di precisione a 80 bit, dando buone probabilità che i risultati siano spenti di un bit a causa del round-off. Anche il motivo per cui @ J99 può ottenere un risultato consistente con/fp: veloce in VS2010, il suo compilatore usa ancora FPU e/fp: fast usa direttamente il risultato di FSIN.

C'è molto bang in questa funzione, controlla il video di Jim Hogg nell'URL collegato per scoprire come sfruttarlo.