MinGW GCC 4.9.1 e determinismo in virgola mobile

Question

Ho scritto un piccolo programma per calcolare la norma euclidea di un vettore a 3 coordinate. Eccolo:

#include <array> 
#include <cmath> 
#include <iostream> 

template<typename T, std::size_t N> 
auto norm(const std::array<T, N>& arr) 
    -> T 
{ 
    T res{}; 
    for (auto value: arr) 
    { 
     res += value * value; 
    } 
    return std::sqrt(res); 
} 

int main() 
{ 
    std::array<double, 3u> arr = { 4.0, -2.0, 6.0 }; 
    std::cout << norm(arr) - norm(arr) << '\n'; 
} 

Sul mio computer, stampa -1.12323e-016.

---

So che i tipi in virgola mobile devono essere gestiti con cura. Tuttavia, pensavo che le operazioni in virgola mobile fossero almeno in qualche modo deterministiche. This article merito a virgola mobile determinismo precisa che:

Alcune delle cose che sono garantiti sono i risultati di addizione, sottrazione, moltiplicazione, divisione e radice quadrata. I risultati di queste operazioni sono garantiti per essere il risultato esatto arrotondato correttamente (ne riparleremo più avanti) quindi se fornisci gli stessi valori di input, con le stesse impostazioni globali e la stessa precisione di destinazione, ti viene garantito lo stesso risultato.

Come si può vedere, le sole operazioni eseguite da questo programma su valori in virgola mobile sono addizione, sottrazione, moltiplicazione e radice quadrata. Se mi fido dell'articolo che ho citato sopra, considerando che gira in un singolo thread e che non cambio le modalità di arrotondamento o qualsiasi altra cosa legata in virgola mobile, ho pensato che norm(arr) - norm(arr) sarebbe 0 poiché faccio esattamente le stesse operazioni sul stessi valori due volte.

Le mie supposizioni sono errate o si tratta di un caso di compilatore non strettamente conforme alla matematica in virgola mobile IEEE? Attualmente sto usando MinGW-W64 GCC 4.9.1 32 bit (ho provato ogni livello di ottimizzazione da -O0 a -O3). Lo stesso programma con MinGW-W64 GCC 4.8.x visualizzato 0, che è quello che mi sarei aspettato.

---

EDIT: ho smontato il codice. Non invierò l'intero assembly generato perché è troppo grande. Tuttavia, credo che la parte rilevante è qui:

call ___main 
fldl LC0 
fstpl -32(%ebp) 
fldl LC1 
fstpl -24(%ebp) 
fldl LC2 
fstpl -16(%ebp) 
leal -32(%ebp), %eax 
movl %eax, (%esp) 
call __Z4normIdLj3EET_RKSt5arrayIS0_XT0_EE 
fstpl -48(%ebp) 
leal -32(%ebp), %eax 
movl %eax, (%esp) 
call __Z4normIdLj3EET_RKSt5arrayIS0_XT0_EE 
fsubrl -48(%ebp) 
fstpl (%esp) 
movl $__ZSt4cout, %ecx 
call __ZNSolsEd 
subl $8, %esp 
movl $10, 4(%esp) 
movl %eax, (%esp) 
call __ZStlsISt11char_traitsIcEERSt13basic_ostreamIcT_ES5_c 
movl $0, %eax 
movl -4(%ebp), %ecx 
.cfi_def_cfa 1, 0 
leave 

Come si può vedere, __Z4normIdLj3EET_RKSt5arrayIS0_XT0_EE è chiamato due volte, quindi, non è inline. Tuttavia, non capisco tutto e non posso dire quale sia il problema.

Answer 1

Come @MatthiasB ha sottolineato, questo sembra essere un problema di gcc memorizzare temporaneamente un valore in virgola mobile a 80 bit in un registro a 64 bit/posizione di memoria. Si consideri il seguente programma semplificato che riproduce ancora il problema:

#include <cmath> 
#include <iostream> 

double norm() { 
    double res = 4.0 * 4.0 + (-2.0 * -2.0) + (6.0 * 6.0); 
    return std::sqrt(res); 
} 

int main() { 
    std::cout << norm() - norm() << '\n'; 
    return 0; 
} 

Il codice montaggio della parte essenziale norm() - norm() si presenta così (usando 32 bit mingw 4.8.0 compilatore)

... 
call __Z4normv  ; call norm() 
fstpl -16(%ebp)  ; store result (80 bit) in temporary (64 bit!) 
call __Z4normv  ; call norm() again 
fsubrl -16(%ebp)  ; subtract result (80 bit) from temporary (64 bit!) 
... 

In sostanza, vorrei prendere in considerazione questo un bug gcc, ma sembra essere un complicated topic ...

Answer 2

C'è una differenza di precisione di un numero in virgola mobile a seconda di dove è memorizzato. Se il compilatore mantiene una variabile in un registro, ha maggiore precisione come variabile memorizzata. Si può cercare di forzare le variabili da memorizzare nella memoria prima, come:

int main() 
{ 
    std::array<double, 3u> arr = { 4.0, -2.0, 6.0 }; 
    volatile double v1 = norm(arr); 
    volatile double v2 = norm(arr); 
    std::cout << v1 - v2 << '\n'; 
} 

Questo ti dà il risultato atteso di 0.