Obliczenia zmiennoprzecinkowe IEEE-754, równość i zwężenie

Question

W poniższym kodzie funkcje foo1, foo2 i foo3 mają być równoważne. Jednak, gdy uruchamiamy foo3 nie kończy się z pętli, czy istnieje powód, dlaczego tak jest?

template <typename T> 
T foo1() 
{ 
    T x = T(1); 
    T y = T(0); 
    for (;;) 
    { 
     if (x == y) break; 
     y = x; 
     ++x; 
    } 
    return x; 
} 

template <typename T> 
T foo2() 
{ 
    T x = T(0); 
    for (;;) 
    { 
     T y = x + T(1); 
     if (!(x != y)) break; 
     ++x; 
    } 
    return x; 
} 

template <typename T> 
T foo3() 
{ 
    T x = T(0); 
    while (x != (x + T(1))) ++x; 
    return x; 
} 

int main() 
{ 
    printf("1 float: %20.5f\n", foo1<float>()); 
    printf("2 float: %20.5f\n", foo2<float>()); 
    printf("3 float: %20.5f\n", foo3<float>()); 
    return 0; 
} 

Uwaga: ta została opracowana przy użyciu VS2010 z/fp precyzyjne w trybie zwolnienia. Nie wiem, jak GCC itp. Traktowałby ten kod, każda informacja byłaby świetna. Czy to może być problem, w przypadku gdy w foo3 wartości x i x + stają się w jakiś sposób NaN?

Answer 1

Co się dzieje jest najprawdopodobniej następujące. Na łuku x86, pośrednie obliczenia można wykonać z 80 bitami dokładności (długie podwójne to odpowiedni typ C/C++). Kompilator używa wszystkich 80 bitów do operacji (+1) i operacji (! =), Ale obcina wyniki przed zapisaniem.

Więc co kompilator naprawdę jest taka:

while ((long double)(x) != ((long double)(x) + (long double)(1))) { 
    x = (float)((long double)(x) + (long double)(1)); 
} 

To jest absolutnie nie-zgodne ze standardem IEEE-a powoduje niekończące się bóle głowy dla wszystkich, ale jest to ustawienie domyślne dla MSVC. Użyj flagi kompilatora /fp:strict, aby wyłączyć to zachowanie.

To jest moje wspomnienie o problemie od około 10 lat temu, więc proszę mi wybaczyć, jeśli nie jest to jakoś całkowicie poprawne. See this for the official Microsoft documentation.

EDIT Byłem bardzo zaskoczony tym, że g ++ domyślnie eksponatów dokładnie takie samo zachowanie (na i386 linux, ale nie z np -mfpmath = sse).