Niespójność obliczeń sinusów w Visual C++ 2012?

Question

Rozważmy następujący kod:

// 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; 
} 

Kiedy skompilowany z Visual C++ 2012 (cl fputest.cpp), a program jest wykonywany, wyjście jest następujące:

3FEDF640D8D36174 
3FEDF640D8D36175 

Pytania:

• Dlaczego te dwie wartości są różne?
• Czy jest możliwe wydanie niektórych opcji kompilatora, aby obliczone wartości sinusoidalne były dokładnie takie same?

Answer 1

Ten problem nie jest spowodowany konwersją z długich podwójnych na podwójne. Przyczyną może być niedokładność procedury sin w bibliotece matematycznej.

Instrukcja jest podana w celu uzyskania wyniku w zakresie 1 ULP (w długim podwójnym formacie) dla operandów w swoim zakresie (według Intel 64 i IA-32 Architectures Software Developer's Manual, październik 2011, tom 1, 8.3.10), w trybie od okrągłego do najbliższego. Na Intel Core i7, fsin wartości pytający użytkownika, -5.07121364272633190495298549649305641651153564453125 lub -0x1.448ec3aaa278dp + 2, produkuje 0xe.fb206c69b0ba402p-4. Z tej heksadecymalnej liczby można łatwo zauważyć, że ostatnie 11 bitów to 100 0000 0010. Są to bity, które zostaną zaokrąglone podczas konwersji z długiego podwójnego. Jeśli są większe niż 100 0000 0000, liczba zostanie zaokrąglona w górę. Są więksi. Dlatego wynikiem konwersji tej długiej podwójnej wartości na double jest 0xe.fb206c69b0ba8p-4, co równa się 0x1.df640d8d36175p-1 i 0.93631021832247418590355891865328885614871978759765625. Zauważ również, że nawet jeśli wynik byłby o jeden ULP niższy, ostatnie 11 bitów wciąż byłby większy niż 100 0000 0000 i nadal zaokrąglałby. Dlatego wynik ten nie powinien się różnić w przypadku procesorów Intela zgodnych z powyższą dokumentacją.

Porównaj to do obliczania podwójnej precyzji sinus bezpośrednio, wykorzystując idealne sin rutynę, która produkuje poprawnie zaokrąglone wyniki. Sinus wartości wynosi około 0,93631021832247413051857150785044253634581268961333520518023697738674775240815140702992025520721336793516756640679315765619707343171517531053811196321335899848286682535203710849065933755262347468763562 (obliczony dla Maple 10). Podwójne najbliżej tego jest 0x1.df640d8d36175p-1. Jest to ta sama wartość, którą uzyskaliśmy, przekształcając wynik na 22,000.

Dlatego rozbieżność nie jest spowodowana przez konwersję długich podwójnych na podwójne; przekształcenie podwójnego wyniku o podwójnej wartości daje wynik dokładnie taki sam, jak w przypadku idealnej rutyny podwójnej precyzji o dokładności sin.

Nie mamy specyfikację dokładności rutyny sin używanego przez pakiet Visual Studio pytający użytkownika. W bibliotekach komercyjnych powszechne jest zezwalanie na błędy rzędu 1 ULP lub kilku ULP. Obserwuj, jak blisko jest sinus do punktu, w którym wartość podwójnej precyzji jest zaokrąglona: Jest to 0,498864 ULP (podwójna precyzja ULP) z dala od podwójnego, czyli wynosi 0,001136 ULP od punktu, w którym zaokrąglenia się zmieniają. Dlatego nawet bardzo niewielka niedokładność w procedurze sin spowoduje, że zwróci 0x1.df640d8d36174p-1 zamiast bliższego 0x1.df640d8d36175p-1.

Dlatego przypuszczam, że źródłem rozbieżności jest bardzo mała niedokładność w procedurze sin.

Answer 2

(Uwaga:..! Jak wspomniano w komentarzach, to nie działa na VC2012 Zostawiłem ją tutaj informacje ogólne nie polecam powołując się na wszystko, co zależy od poziomu optymalizacji tak)

Nie mam VS2012, ale na kompilatorze VS2010 można podać /fp:fast w wierszu poleceń, a następnie otrzymam te same wyniki. Powoduje to, że kompilator generuje "szybki" kod, który niekoniecznie musi dokładnie odpowiadać wymaganemu zaokrągleniu i regułom w C++, ale który pasuje do obliczania języka asemblerowego.

Nie mogę tego wypróbować w VS2012, ale wyobrażam sobie, że ma tę samą opcję.

To tylko wydaje się działać w zoptymalizowanej kompilacji również z /Ox jako opcja.

Answer 3

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

Jak David mentioned in the comment, rozbieżność pochodzi od ruchu danych w FP zarejestrować do lokalizacji pamięci (zarejestrować/RAM) o różnej wielkości. I nie zawsze jest to zadanie; nawet inna operacja w trybie zmiennoprzecinkowym może być wystarczająca do przepłukania rejestru FP, czyniąc wszelkie próby zagwarantowania, że ​​dana wartość będzie daremna. Jeśli trzeba zrobić porównanie, może być w stanie złagodzić niektóre z tych zmuszając wszystkich wyników do miejsca pamięci w następujący sposób:

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

Jednak nawet to może nie działać. Najlepszym podejściem jest zaakceptowanie, że każda operacja zmiennoprzecinkowa będzie jedynie "w większości dokładna" i że z punktu widzenia programistów ten błąd będzie faktycznie nieprzewidywalny i losowy, nawet jeśli ta sama operacja jest wykonywana wielokrotnie. Zamiast

if (F1 == F2) 

należy użyć coś do skutku

if (isArbitrarilyClose(F1, F2)) 

lub

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

gdzie n to mała liczba.

Answer 4

Generator kodu w VS2012 został znacznie zmieniony w celu wspierania automatic vectorization. Częścią tej zmiany jest to, że matematyka zmiennoprzecinkowa x86 jest teraz wykonywana w SSE2 i nie używa już FPU, ponieważ kod FPU nie może być wektoryzowany. SSE2 oblicza z 64-bitową precyzją zamiast 80-bitową precyzją, dając dobre kursy, że wyniki są wyłączone o jeden bit z powodu zaokrąglenia. Również powód, dla którego @ J99 może uzyskać spójny wynik z/fp: fast w VS2010, jego kompilator nadal używa FPU i/fp: fast używa bezpośrednio wyniku FSIN.

Istnieje wiele hukiem w tej funkcji, sprawdź wideo Jim Hogg pod adresem URL połączonego aby dowiedzieć się, jak z niej skorzystać.