Strange uint32_t для преобразования массива с плавающей запятой

У меня есть следующий фрагмент кода:

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

Если я скомпилирую его в MS Visual Studio 2015, я вижу, что вывод:

WTF??!!
magic = 2570980352.0000000000
magic = 2570980608.0000000000

Итак, последний элемент arr_dst отличается от предыдущего, но эти два значения были получены путем преобразования того же значения, которое заполняет массив arr! Это ошибка?

Я заметил, что если я модифицирую цикл преобразования следующим образом, я получаю результат "ОК":

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

Итак, это, вероятно, проблема с оптимизацией векторизации.

Это поведение не воспроизводится на gcc 4.8. Любые идеи?

Ответ 1

Я провел исследование по внедрению PowerPC (Freescale MCP7450), поскольку они IMHO гораздо лучше документированы, чем любой вуду Intel.

Как оказалось, блок с плавающей точкой, FPU и векторный блок могут иметь различное округление для операций с плавающей запятой. FPU может быть настроен на использование одного из четырех режимов округления; округление до ближайшего (по умолчанию), усечение, к положительной бесконечности и к отрицательной бесконечности. Однако векторная единица может только округлять до ближайшего, с несколькими инструкциями выбора, имеющими конкретные правила округления. Внутренняя точность FPU составляет 106 бит. Единица вектора соответствует IEEE-754, но документация не указала намного больше.

Глядя на ваш результат, преобразование 2570980608 ближе к исходному целому числу, предполагая, что FPU имеет лучшую внутреннюю точность, чем векторная единица ИЛИ различные режимы округления.

Ответ 2

32-битный двоичный float IEEE-754, такой как MSVС++, использует только 6-7 десятичных цифр точности. Ваше начальное значение находится в пределах диапазона этого типа, но похоже, что он не может быть точно представлен этим типом, как это действительно имеет место для большинства значений типа uint32_t.

В то же время блок с плавающей точкой процессора x86 или x86_64 использует более широкое представление, отличное от MSVС++ 64-бит double. Кажется вероятным, что после выхода цикла, последний вычисленный элемент массива остается в регистре FPU в форме расширенной точности. Затем программа может использовать это значение непосредственно из регистра вместо того, чтобы читать его обратно из памяти, которое оно обязано делать с предыдущими элементами.

Если программа выполняет сравнение ==, продвигая более узкое представление к более широкому, а не наоборот, то эти два значения могут действительно сравнивать неравные, так как округление от расширенной точности до float и обратно теряет точность. В любом случае оба значения преобразуются в тип double при передаче на printf(); если действительно они сравнивают неравные, то, скорее всего, результаты этих преобразований тоже отличаются.

Я не использую параметры компиляции MSVС++, но, скорее всего, есть тот, который отменяет это поведение. Такие варианты иногда бывают такими именами, как "строгая математика" или "строгое fp". Однако имейте в виду, что включение такой опции (или отключение ее противоположности) может оказаться очень дорогостоящим в программе FP-heavy.

Ответ 3

Преобразование между unsigned и float не просто на x86; для него нет единой инструкции (до AVX512). Общим методом является преобразование как подписанное, а затем исправление результата. Существует несколько способов сделать это. (См. этот Q & A для некоторых рутинных векторизованных методов с использованием C intrinsics, не все из которых имеют идеально округленные результаты.)

MSVC векторизовывает первые 128 с одной стратегией и затем использует другую стратегию (которая не будет векторизовать) для последнего скалярного элемента, которая включает преобразование в double, а затем от double до float.

gcc и clang производят результат 2570980608.0 из их векторизованных и скалярных методов. 2570980608 - 2570980487 = 121 и 2570980487 - 2570980352 = 135 (без округления входов/выходов), поэтому gcc и clang дают корректный округленный результат в этом случае (менее 0,5 п.п. ошибки). IDK, если это верно для всех возможных uint32_t (но их всего 2 ^ 32, мы могли бы полностью проверить). Конечный результат MSVC для векторизованного цикла имеет чуть более 0,5 п.п. ошибки, но скалярный метод правильно округлен для этого ввода.

Математика IEEE требует, чтобы + - * / и sqrt производили корректно округленные результаты (менее 0,5 п.п. ошибки), но другие функции (например, log) не имеют таких строгое требование. IDK, что требования для округления для конверсий int- > float, поэтому IDK, если то, что делает MSVC, является строго законным (если вы не использовали /fp:fast или что-то еще).

См. также Брюс Доусон Сообщение о блоге дефиниции с плавающей запятой (часть его отличной серии о математике FP), хотя он не упоминает целочисленные ↔ FP преобразования.

Мы можем видеть в asm, связанном с OP тем, что MSVC сделал (разделился только на интересные инструкции и прокомментирован вручную):

; 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 = [email protected]  (a constant repeated in each of 4 float elements)


  ####### The vectorized unsigned->float conversion strategy:
  [email protected]:                                       ; 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 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 [email protected]    ; }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     [email protected]     ; if(i>=129): always false

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

... несколько более смешных, известных в момент компиляции, время спустя...

 ######## The scalar unsigned->float conversion strategy for the last element
[email protected]:
  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 [email protected][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 [email protected]
  # Yes, this is a loop, which always runs exactly once for the last element

Для сравнения, clang и gcc также не оптимизируют всю вещь во время компиляции, но они понимают, что они не нуждаются в цикле очистки, и просто делают одно скалярное хранилище или конвертируют после соответствующего петли. (clang фактически полностью разворачивает все, пока вы не скажете об этом не.)

См. код в проводник компилятора Godbolt.

gcc только преобразует верхнюю и нижнюю половинки 16b для поплавка отдельно и объединяет их с умножением на 65536 и добавляет.

Стратегия преобразования Clang unsigned → float интересна: она никогда не использует инструкцию cvt вообще. Я думаю, что это заставляет две 16-разрядные половинки беззнакового целого в мантиссу двух поплавков напрямую (с некоторыми трюками, чтобы установить экспоненты (побитовые булевы и ADDPS), а затем добавляет низкую и высокую половину вместе, как gcc.

Конечно, если вы скомпилируете 64-битный код, скалярное преобразование может только нулевое расширение uint32_t до 64-битного и преобразовать его как подписанный int64_t для float. Подписанный int64_t может представлять каждое значение uint32_t, а x86 может конвертировать 64-разрядный подписанный int, чтобы плавать эффективно. Но это не векторизация.