Inny wynik zmiennoprzecinkowy z włączoną optymalizacją - błąd kompilatora?

Poniższy kod działa w programie Visual Studio 2008 z optymalizacją i bez niej. Ale działa tylko na g ++ bez optymalizacji (O0).

#include <cstdlib>
#include <iostream>
#include <cmath>

double round(double v, double digit)
{
    double pow = std::pow(10.0, digit);
    double t = v * pow;
    //std::cout << "t:" << t << std::endl;
    double r = std::floor(t + 0.5);
    //std::cout << "r:" << r << std::endl;
    return r / pow;
}

int main(int argc, char *argv[])
{
    std::cout << round(4.45, 1) << std::endl;
    std::cout << round(4.55, 1) << std::endl;
}

Wynik powinien być:

4.5
4.6

Ale g ++ z optymalizacją ( O1- O3) wyświetli:

4.5
4.5

Jeśli dodam volatilesłowo kluczowe przed t, to działa, więc czy może wystąpić jakiś błąd optymalizacji?

Test na g ++ 4.1.2 i 4.4.4.

Oto wynik na ideone: http://ideone.com/Rz937

Opcja, którą testuję na g ++, jest prosta:

g++ -O2 round.cpp

Ciekawszy wynik, nawet jeśli włączę /fp:fastopcję na Visual Studio 2008, wynik nadal jest poprawny.

Kolejne pytanie:

Zastanawiałem się, czy zawsze powinienem włączać -ffloat-storeopcję?

Ponieważ testowana przeze mnie wersja g ++ jest dostarczana z CentOS / Red Hat Linux 5 i CentOS / Redhat 6 .

Skompilowałem wiele moich programów na tych platformach i obawiam się, że spowoduje to nieoczekiwane błędy w moich programach. Wydaje się, że trochę trudno jest zbadać cały mój kod C ++ i używane biblioteki, czy mają takie problemy. Jakieś sugestie?

Czy ktoś jest zainteresowany tym, dlaczego nawet /fp:fastwłączony, Visual Studio 2008 nadal działa? Wygląda na to, że Visual Studio 2008 jest bardziej niezawodny w tym problemie niż g ++?

Procesory Intel x86 używają wewnętrznie 80-bitowej rozszerzonej precyzji, podczas gdy doublezwykle mają 64-bitową szerokość. Różne poziomy optymalizacji wpływają na to, jak często wartości zmiennoprzecinkowe z procesora są zapisywane w pamięci, a tym samym zaokrąglane z precyzji 80-bitowej do precyzji 64-bitowej.

Użyj -ffloat-storeopcji gcc, aby uzyskać te same wyniki zmiennoprzecinkowe z różnymi poziomami optymalizacji.

Alternatywnie, użyj long doubletypu, który zwykle ma 80-bitową szerokość w gcc, aby uniknąć zaokrąglania z 80-bitowej do 64-bitowej precyzji.

man gcc mówi wszystko:

   -ffloat-store
       Do not store floating point variables in registers, and inhibit
       other options that might change whether a floating point value is
       taken from a register or memory.

       This option prevents undesirable excess precision on machines such
       as the 68000 where the floating registers (of the 68881) keep more
       precision than a "double" is supposed to have.  Similarly for the
       x86 architecture.  For most programs, the excess precision does
       only good, but a few programs rely on the precise definition of
       IEEE floating point.  Use -ffloat-store for such programs, after
       modifying them to store all pertinent intermediate computations
       into variables.

W kompilacjach x86_64 kompilatory używają rejestrów SSE dla floati doubledomyślnie, dzięki czemu nie jest używana rozszerzona precyzja i ten problem nie występuje.

gccopcja kompilatora-mfpmath kontroluje to.

Wynik powinien wyglądać następująco: 4.5 4.6 Tak wyglądałyby dane wyjściowe, gdybyś miał nieskończoną precyzję lub gdybyś pracował z urządzeniem, które używało dziesiętnej zamiast binarnej reprezentacji zmiennoprzecinkowej. Ale tak nie jest. Większość komputerów używa binarnego standardu zmiennoprzecinkowego IEEE.

Jak Maxim Yegorushkin już zauważył w swojej odpowiedzi, część problemu polega na tym, że wewnętrznie twój komputer używa 80-bitowej reprezentacji zmiennoprzecinkowej. To tylko część problemu. Podstawą problemu jest to, że żadna liczba w postaci n.nn5 nie ma dokładnej binarnej reprezentacji zmiennoprzecinkowej. Te przypadki narożne są zawsze niedokładnymi liczbami.

Jeśli naprawdę chcesz, aby Twoje zaokrąglanie mogło niezawodnie zaokrąglić te przypadki narożne, potrzebujesz algorytmu zaokrąglania, który uwzględnia fakt, że n.n5, n.nn5 lub n.nnn5 itd. (Ale nie n.5) jest zawsze niedokładny. Znajdź przypadek narożny, który określa, czy jakaś wartość wejściowa jest zaokrąglana w górę, czy w dół, i zwraca wartość zaokrągloną w górę lub w dół na podstawie porównania z tym przypadkiem narożnym. I musisz uważać, aby optymalizujący kompilator nie umieścił tego znalezionego narożnika w rejestrze o rozszerzonej precyzji.

Zobacz, w jaki sposób program Excel pomyślnie zaokrągla liczby zmienne, mimo że są one nieprecyzyjne? dla takiego algorytmu.

Lub możesz po prostu żyć z faktem, że narożniki czasami będą błędnie zaokrąglane.

Różne kompilatory mają różne ustawienia optymalizacji. Niektóre z szybszych ustawień optymalizacji nie zachowują ścisłych reguł zmiennoprzecinkowych zgodnie z IEEE 754 . Visual Studio ma specyficzne ustawienie, /fp:strict, /fp:precise, /fp:fast, gdzie /fp:fastjest niezgodny ze standardem na to, co można zrobić. Może się okazać, że ta flaga steruje optymalizacją w takich ustawieniach. Możesz również znaleźć podobne ustawienie w GCC, które zmienia zachowanie.

Jeśli tak jest, jedyną różnicą między kompilatorami jest to, że GCC domyślnie szukałby najszybszego zachowania zmiennoprzecinkowego przy wyższych optymalizacjach, podczas gdy Visual Studio nie zmienia zachowania zmiennoprzecinkowego przy wyższych poziomach optymalizacji. Dlatego niekoniecznie musi to być rzeczywisty błąd, ale zamierzone zachowanie opcji, o której nie wiedziałeś, że włączasz.

Do tych, którzy nie mogą odtworzyć błędu: nie odkomentuj zakomentowanych instrukcji debugowania, mają one wpływ na wynik.

Oznacza to, że problem jest związany z instrukcjami debugowania. Wygląda na to, że wystąpił błąd zaokrąglania spowodowany ładowaniem wartości do rejestrów podczas instrukcji wyjściowych, dlatego inni odkryli, że można to naprawić za pomocą-ffloat-store

Kolejne pytanie:

Zastanawiałem się, czy zawsze powinienem włączać -ffloat-storeopcję?

Aby być niepoważny, to musi być jakiś powód, że niektórzy programiści nie włączyć -ffloat-store, w przeciwnym razie opcja nie istnieje (podobnie, tam musi być jakiś powód, że niektórzy programiści nie włącza się -ffloat-store). Nie radziłbym zawsze go włączać lub wyłączać. Włączenie go zapobiega niektórym optymalizacjom, ale wyłączenie go pozwala na zachowanie, które otrzymujesz.

Ale ogólnie rzecz biorąc, istnieje pewne niedopasowanie między binarnymi liczbami zmiennoprzecinkowymi (używanymi przez komputer) a dziesiętnymi liczbami zmiennoprzecinkowymi (które ludzie są zaznajomieni), a ta niezgodność może powodować podobne zachowanie do tego, co otrzymujesz (aby było jasne, zachowanie które otrzymujesz, nie jest spowodowane niedopasowaniem, ale podobne zachowanie może być). Rzecz w tym, że skoro masz już pewne niejasności, gdy masz do czynienia z zmiennoprzecinkowymi, nie mogę powiedzieć, -ffloat-storeże to poprawi lub pogorszy.

Zamiast tego możesz przyjrzeć się innym rozwiązaniom problemu, który próbujesz rozwiązać (niestety, Koenig nie wskazuje na rzeczywisty artykuł, a ja nie mogę znaleźć dla niego oczywistego „kanonicznego” miejsca, więc będę musiał wysłać Cię do Google ).

Jeśli nie zaokrąglasz w celach wyjściowych, prawdopodobnie spojrzałbym na std::modf()(in cmath) i std::numeric_limits<double>::epsilon()(in limits). Zastanawiając się nad oryginalną round()funkcją, uważam, że czystsze byłoby zastąpienie wywołania std::floor(d + .5)funkcji wywołaniem tej funkcji:

// this still has the same problems as the original rounding function
int round_up(double d)
{
    // return value will be coerced to int, and truncated as expected
    // you can then assign the int to a double, if desired
    return d + 0.5;
}

Myślę, że sugeruje to następującą poprawę:

// this won't work for negative d ...
// this may still round some numbers up when they should be rounded down
int round_up(double d)
{
    double floor;
    d = std::modf(d, &floor);
    return floor + (d + .5 + std::numeric_limits<double>::epsilon());
}

Prosta uwaga: std::numeric_limits<T>::epsilon()jest definiowana jako „najmniejsza liczba dodana do 1, która tworzy liczbę różną od 1”. Zwykle musisz użyć względnego epsilon (tj. Skalować epsilon w jakiś sposób, aby uwzględnić fakt, że pracujesz z liczbami innymi niż „1”). Suma d, .5i std::numeric_limits<double>::epsilon()powinien być bliski 1, więc grupowanie Oznacza to ponadto, że std::numeric_limits<double>::epsilon()będzie o odpowiedniej wielkości za to, co robimy. Jeśli już, std::numeric_limits<double>::epsilon()będzie za duża (gdy suma wszystkich trzech jest mniejsza niż jeden) i może spowodować, że będziemy zaokrąglać niektóre liczby w górę, gdy nie powinniśmy.

W dzisiejszych czasach powinieneś to rozważyć std::nearbyint().

Zaakceptowana odpowiedź jest poprawna, jeśli kompilujesz do celu x86, który nie zawiera SSE2. Wszystkie nowoczesne procesory x86 obsługują SSE2, więc jeśli możesz z tego skorzystać, powinieneś:

-mfpmath=sse -msse2 -ffp-contract=off

Rozbijmy to.

-mfpmath=sse -msse2. To wykonuje zaokrąglanie przy użyciu rejestrów SSE2, co jest znacznie szybsze niż przechowywanie każdego pośredniego wyniku w pamięci. Zauważ, że jest to już domyślne ustawienie w GCC dla x86-64. Z wiki GCC :

Na bardziej nowoczesnych procesorach x86, które obsługują SSE2, określenie opcji kompilatora -mfpmath=sse -msse2zapewnia, że ​​wszystkie operacje zmiennoprzecinkowe i podwójne są wykonywane w rejestrach SSE i prawidłowo zaokrąglane. Opcje te nie mają wpływu na ABI i dlatego powinny być stosowane, gdy tylko jest to możliwe, w celu uzyskania przewidywalnych wyników liczbowych.

-ffp-contract=off. Jednak kontrolowanie zaokrąglania nie wystarczy, aby uzyskać dokładne dopasowanie. Instrukcje FMA (łączone mnożenie i dodawanie) mogą zmienić zachowanie zaokrąglania w porównaniu z ich nieskondensowanymi odpowiednikami, więc musimy je wyłączyć. Jest to domyślne ustawienie w Clang, a nie GCC. Jak wyjaśnia ta odpowiedź :

FMA ma tylko jedno zaokrąglenie (efektywnie utrzymuje nieskończoną precyzję dla wewnętrznego tymczasowego wyniku mnożenia), podczas gdy ADD + MUL ma dwa.

Wyłączając FMA, otrzymujemy wyniki, które dokładnie pasują do debugowania i wydania, kosztem pewnej wydajności (i dokładności). Nadal możemy korzystać z innych zalet wydajnościowych SSE i AVX.

Zagłębiłem się bardziej w ten problem i mogę wprowadzić więcej szczegółów. Po pierwsze, dokładne reprezentacje 4,45 i 4,55 zgodnie z gcc na x84_64 są następujące (z libquadmath wypisuje ostatnią precyzję):

float 32:   4.44999980926513671875
double 64:  4.45000000000000017763568394002504646778106689453125
doublex 80: 4.449999999999999999826527652402319290558807551860809326171875
quad 128:   4.45000000000000000000000000000000015407439555097886824447823540679418548304813185723105561919510364532470703125

float 32:   4.55000019073486328125
double 64:  4.54999999999999982236431605997495353221893310546875
doublex 80: 4.550000000000000000173472347597680709441192448139190673828125
quad 128:   4.54999999999999999999999999999999984592560444902113175552176459320581451695186814276894438080489635467529296875

Jak powiedział powyżej Maxim , problem wynika z rozmiaru 80 bitów rejestrów FPU.

Ale dlaczego problem nigdy nie występuje w systemie Windows? na IA-32 jednostka FPU x87 została skonfigurowana do używania wewnętrznej precyzji mantysy wynoszącej 53 bity (co odpowiada całkowitemu rozmiarowi 64 bitów:) double. W systemach Linux i Mac OS zastosowano domyślną precyzję 64 bitów (odpowiednik całkowitego rozmiaru 80 bitów:) long double. Zatem problem powinien być możliwy lub nie na tych różnych platformach poprzez zmianę słowa kontrolnego FPU (zakładając, że sekwencja instrukcji wywoła błąd). Problem został zgłoszony do gcc jako błąd 323 (przeczytaj przynajmniej komentarz 92!).

Aby pokazać precyzję mantysy w systemie Windows, możesz skompilować to w 32 bitach za pomocą VC ++:

#include "stdafx.h"
#include <stdio.h>  
#include <float.h>  

int main(void)
{
    char t[] = { 64, 53, 24, -1 };
    unsigned int cw = _control87(0, 0);
    printf("mantissa is %d bits\n", t[(cw >> 16) & 3]);
}

oraz w systemie Linux / Cygwin:

#include <stdio.h>

int main(int argc, char **argv)
{
    char t[] = { 24, -1, 53, 64 };
    unsigned int cw = 0;
    __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw));
    printf("mantissa is %d bits\n", t[(cw >> 8) & 3]);
}

Zauważ, że za pomocą gcc możesz ustawić precyzję FPU -mpc32/64/80, chociaż jest ona ignorowana w Cygwin. Ale pamiętaj, że zmieni on rozmiar mantysy, ale nie wykładnika, otwierając drzwi dla innych rodzajów różnych zachowań.

W architekturze x86_64, SSE jest używane, jak powiedział tmandry , więc problem nie wystąpi, chyba że wymusisz stary FPU x87 do obliczeń FP -mfpmath=387lub jeśli nie będziesz kompilować w trybie 32-bitowym z -m32(będziesz potrzebować pakietu multilib). Mogę odtworzyć problem w systemie Linux z różnymi kombinacjami flag i wersji gcc:

g++-5 -m32 floating.cpp -O1
g++-8 -mfpmath=387 floating.cpp -O1

Wypróbowałem kilka kombinacji w systemie Windows lub Cygwin z VC ++ / gcc / tcc, ale błąd nigdy się nie pojawił. Przypuszczam, że sekwencja generowanych instrukcji nie jest taka sama.

Na koniec zwróć uwagę, że egzotycznym sposobem uniknięcia tego problemu z 4.45 lub 4.55 byłoby użycie _Decimal32/64/128, ale wsparcie jest naprawdę rzadkie ... Spędziłem dużo czasu tylko po to, aby móc wykonać printf z libdfp!

Osobiście napotkałem ten sam problem idąc w drugą stronę - od gcc do VS. W większości przypadków uważam, że lepiej unikać optymalizacji. Jedyny przypadek, w którym jest to warte zachodu, dotyczy metod numerycznych obejmujących duże tablice danych zmiennoprzecinkowych. Nawet po demontażu często jestem rozczarowany wyborami kompilatorów. Bardzo często po prostu łatwiej jest użyć wewnętrznych funkcji kompilatora lub po prostu napisać zestaw samodzielnie.