Strange uint32_t unosić konwersję tablicy

Question

Mam następujący fragment kodu:

#include <cstdio> 
#include <cstdint> 

static const size_t ARR_SIZE = 129; 

int main() 
{ 
    uint32_t value = 2570980487; 

    uint32_t arr[ARR_SIZE]; 
    for (int x = 0; x < ARR_SIZE; ++x) 
    arr[x] = value; 

    float arr_dst[ARR_SIZE]; 
    for (int x = 0; x < ARR_SIZE; ++x) 
    { 
    arr_dst[x] = static_cast<float>(arr[x]); 
    } 

    printf("%s\n", arr_dst[ARR_SIZE - 1] == arr_dst[ARR_SIZE - 2] ? "OK" : "WTF??!!"); 

    printf("magic = %0.10f\n", arr_dst[ARR_SIZE - 2]); 
    printf("magic = %0.10f\n", arr_dst[ARR_SIZE - 1]); 
    return 0; 
} 

Gdybym go skompilować pod MS Visual Studio 2015 Widzę, że wyjście jest:

WTF??!! 
magic = 2570980352.0000000000 
magic = 2570980608.0000000000 

Więc ostatnia arr_dst element różni się od poprzedniego, ale te dwie wartości zostały uzyskane przez konwersję tej samej wartości, która zapełnia tablicę arr! Czy to błąd?

Zauważyłem, że jeśli zmodyfikować pętlę konwersji w następujący sposób, pojawia się „OK”, wynik:

for (int x = 0; x < ARR_SIZE; ++x) 
{ 
    if (x == 0) 
    x = 0; 
    arr_dst[x] = static_cast<float>(arr[x]); 
} 

Więc to prawdopodobnie jest jakiś problem z optymalizacją Wektoryzacja.

To zachowanie nie jest odtwarzane w gcc 4.8. Jakieś pomysły?

Answer 1

Zrobiłem śledztwo w sprawie impregnacji PowerPC (Freescale MCP7450), ponieważ IMHO są znacznie lepiej udokumentowane niż jakikolwiek inny Intel voodoo.

Okazuje się, że jednostka zmiennoprzecinkowa, jednostka FPU i wektor mogą mieć różne zaokrąglenia dla operacji zmiennoprzecinkowych. Jednostka FPU może być skonfigurowana do używania jednego z czterech trybów zaokrąglania; od okrągłego do najbliższego (domyślnie), obcięte, w kierunku dodatniej nieskończoności iw kierunku ujemnej nieskończoności. Jednostka wektorowa jest jednak w stanie zaokrąglić do najbliższej, z kilkoma wybranymi instrukcjami mającymi określone zasady zaokrąglania. Wewnętrzna precyzja FPU jest 106-bitowa. Jednostka wektorowa spełnia IEEE-754, ale dokumentacja nie zawiera dużo więcej informacji.

Patrząc na twój wynik, konwersja 2570980608 jest bliższa pierwotnej liczbie całkowitej, co sugeruje, że FPU ma lepszą dokładność wewnętrzną niż jednostka wektorowa LUB różne tryby zaokrąglania.

Answer 2

32-bitowy binarny zmiennoprzecinkowy IEEE-754, taki jak MSVC++, zapewnia jedynie 6-7 cyfr dziesiętnych precyzji. Twoja wartość początkowa mieści się w zakresie zakresu tego typu, ale wydaje się, że nie jest dokładnie reprezentowana przez ten typ, tak jak w przypadku większości wartości typu uint32_t.

W tym samym czasie jednostka zmiennoprzecinkowa procesora x86 lub x86_64 używa szerszej reprezentacji nawet niż 64-bitowy kod MSVC++ double. Wydaje się prawdopodobne, że po zakończeniu pętli ostatnio obliczony element tablicy pozostaje w rejestrze FPU, w jego rozszerzonej formie precyzji. Program może następnie użyć tej wartości bezpośrednio z rejestru zamiast odczytywać go z pamięci, co jest obowiązkowe w przypadku poprzednich elementów.

Jeżeli program wykonuje porównanie == poprzez promowanie węższy reprezentację szerszej zamiast na odwrót, to dwie wartości może rzeczywiście porównać nierówne, jako podróż w obie strony z rozszerzonym dokładnością do float iz powrotem traci precyzji. W każdym przypadku obie wartości są konwertowane na typ double po przekazaniu do printf(); jeśli rzeczywiście porównują nierówność, to jest prawdopodobne, że wyniki tych konwersji również się różnią.

Nie mam opcji kompilacji MSVC++, ale najprawdopodobniej istnieje taka, która może przerwać to zachowanie. Takie opcje czasami mają takie nazwy, jak "ścisła matematyka" lub "ścisłe fp". Należy jednak pamiętać, że włączenie takiej opcji (lub wyłączenie jej przeciwnej strony) może być bardzo kosztowne w ciężkim programie FP.

Answer 3

Konwersja między unsigned i float nie jest prosta na x86; nie ma jednej instrukcji dla niego (do AVX512). Typową techniką jest konwersja jako podpisana, a następnie poprawienie wyniku. Istnieje wiele sposobów na zrobienie tego. (Zobacz this Q&A for some manually-vectorized methods with C intrinsics, z których nie wszystkie mają idealnie zaokrąglone wyniki.)

MSVC wektoryzuje pierwsze 128 za pomocą jednej strategii, a następnie stosuje inną strategię (która nie będzie wektoryzować) dla ostatniego elementu skalarnego, który wymaga konwersji do double, a następnie od double do float.

---

gcc szczęk wytworzenia 2570980608.0 wyniku ich wektorowy i metody skalarnych. 2570980608 - 2570980487 = 121 i 2570980487 - 2570980352 = 135 (bez zaokrąglania wejść/wyjść), więc gcc i clang dają w rezultacie prawidłowo zaokrąglony wynik (mniej niż 0,5ulp błędu). IDK, jeśli to prawda dla każdego możliwego uint32_t (ale są tylko 2^32 z nich, we could exhaustively check). Wynik końcowy MSVC dla wektoryzowanej pętli ma nieco ponad 0,5ulp błędu, ale metoda skalarna jest poprawnie zaokrąglona dla tego wejścia.

IEEE matematyka żąda +-*/ i sqrt produkować poprawnie zaokrąglone wyniki (mniej niż 0.5ulp błędu), ale pozostałe funkcje (jak log) nie mają takiej ścisłe wymagania. IDK, jakie są wymagania dotyczące zaokrąglania konwersji int-> float, więc IDK, jeśli to, co robi MSVC, jest ściśle legalne (jeśli nie używasz /fp:fast lub cokolwiek innego).

Zobacz także Bruce Dawson's Floating-Point Determinism blog post (część jego doskonałej serii o matematyce FP), chociaż nie wspomina o liczbach całkowitych < -> Konwersje FP.

---

Widzimy w asm połączonego przez OP co MSVC zrobił (uproszczoną tylko ciekawych wskazówek i skomentował ręcznie):

; Function compile flags: /Ogtp 
# assembler macro constants 
_arr_dst$ = -1040     ; size = 516 
_arr$ = -520      ; size = 516 
_main PROC      ; COMDAT 

    00013  mov  edx, 129 
    00018  mov  eax, -1723986809 ; this is your unsigned 2570980487 
    0001d  mov  ecx, edx 
    00023  lea  edi, DWORD PTR _arr$[esp+1088] ; edi=arr 
    0002a  rep stosd    ; memset in chunks of 4B 
    # arr[0..128] = 2570980487 at this point 

    0002c  xor  ecx, ecx  ; i = 0 
    # xmm2 = 0.0 in each element (i.e. all-zero) 
    # xmm3 = __xmm@4f8000004f8000004f8000004f800000 (a constant repeated in each of 4 float elements) 


    ####### The vectorized unsigned->float conversion strategy: 
    $LL7@main:          ; do{ 
    00030  movups xmm0, XMMWORD PTR _arr$[esp+ecx*4+1088] ; load 4 uint32_t 
    00038  cvtdq2ps xmm1, xmm0     ; SIGNED int to Single-precision float 
    0003b  movaps xmm0, xmm1 
    0003e  cmpltps xmm0, xmm2     ; xmm0 = (xmm0 < 0.0) 
    00042  andps xmm0, xmm3     ; mask the magic constant 
    00045  addps xmm0, xmm1     ; x += (x<0.0) ? magic_constant : 0.0f; 
    # There's no instruction for converting from unsigned to float, so compilers use inconvenient techniques like this to correct the result of converting as signed. 
    00048  movups XMMWORD PTR _arr_dst$[esp+ecx*4+1088], xmm0 ; store 4 floats to arr_dst 
    ; and repeat the same thing again, with addresses that are 16B higher (+1104) 
    ; i.e. this loop is unrolled by two 

    0006a  add  ecx, 8   ; i+=8 (two vectors of 4 elements) 
    0006d  cmp  ecx, 128 
    00073  jb SHORT $LL7@main ; }while(i<128) 

#### End of vectorized loop 
# and then IDK what MSVC smoking; both these values are known at compile time. Is /Ogtp not full optimization? 
# I don't see a branch target that would let execution reach this code 
# other than by falling out of the loop that ends with ecx=128 
    00075  cmp  ecx, edx 
    00077  jae  $LN21@main  ; if(i>=129): always false 

    0007d  sub  edx, ecx  ; edx = 129-128 = 1 

... niektóre bardziej znane śmieszne -w-czasie kompilacji skoki później ...

######## The scalar unsigned->float conversion strategy for the last element 
$LC15@main: 
    00140  mov  eax, DWORD PTR _arr$[esp+ecx*4+1088] 
    00147  movd xmm0, eax 
    # eax = xmm0[0] = arr[128] 
    0014b  cvtdq2pd xmm0, xmm0  ; convert the last element TO DOUBLE 
    0014f  shr  eax, 31   ; shift the sign bit to bit 1, so eax = 0 or 1 
    ; then eax indexes a 16B constant, selecting either 0 or 0x41f0... (as whatever double that represents) 
    00152  addsd xmm0, QWORD PTR __xmm@41f00000000000000000000000000000[eax*8] 
    0015b  cvtpd2ps xmm0, xmm0  ; double -> float 
    0015f  movss DWORD PTR _arr_dst$[esp+ecx*4+1088], xmm0 ; and store it 

    00165  inc  ecx   ; ++i; 
    00166  cmp  ecx, 129  ; } while(i<129) 
    0016c  jb SHORT $LC15@main 
    # Yes, this is a loop, which always runs exactly once for the last element 

---

drodze O f porównania, clang i gcc również nie optymalizują całości podczas kompilacji, ale zdają sobie sprawę, że nie potrzebują czyszczenia pętli, i po prostu zrobić jeden sklep skalarny lub przekonwertować po odpowiednich pętlach. (Klang faktycznie w pełni rozwija wszystko, chyba że powiesz mu, żeby tego nie robił.)

Zobacz kod na Godbolt compiler explorer.

gcc po prostu konwertuje połówki górną i dolną 16b w celu unoszenia się osobno i łączy je z pomnożeniem przez 65536 i dodaniem.

Strategia konwersji Clang'a unsigned ->float jest interesująca: w ogóle nie korzysta z instrukcji cvt. Myślę, że to wypycha dwie 16-bitowe połówki niepodpisanej liczby całkowitej bezpośrednio do mantys z dwóch pływaków (z kilkoma sztuczkami, aby ustawić wykładniki (bitowe boolean stuff i ADDPS), a następnie dodaje niską i wysoką połowę razem jak gcc.

Oczywiście, jeśli kompilujesz do kodu 64-bitowego, konwersja skalarna może zerowe rozszerzenie uint32_t na 64-bitowe i przekonwertować ją jako sygnaturę int64_t do float. Podpisany int64_t może reprezentować każdą wartość uint32_t, a x86 może przekonwertować 64-bitową int podpisaną na float. Ale to nie wektoryzuje.