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, чтобы плавать эффективно. Но это не векторизация.
