Странное преобразование массива uint32_t в float
У меня есть следующий фрагмент кода:
#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. Есть идеи?
3 ответа
Я провел расследование по поводу реализации PowerPC (Freescale MCP7450), так как они ИМХО гораздо лучше задокументированы, чем любая вуду, которую придумывает Intel.
Как оказалось, модуль с плавающей запятой, FPU и векторный модуль могут иметь различное округление для операций с плавающей запятой. FPU может быть настроен на использование одного из четырех режимов округления; округление до ближайшего (по умолчанию) усечения до положительной бесконечности и отрицательной бесконечности. Тем не менее, векторный блок может округляться только до ближайшего, с несколькими инструкциями выбора, имеющими определенные правила округления. Внутренняя точность FPU составляет 106 бит. Векторное устройство соответствует стандарту IEEE-754, но в документации не указано намного больше.
Если посмотреть на ваш результат, преобразование 2570980608 ближе к исходному целому числу, что говорит о том, что FPU имеет лучшую внутреннюю точность, чем векторный блок ИЛИ другие режимы округления.
Преобразование между unsigned а также float не просто на x86; для него нет единой инструкции (до AVX512). Обычная техника - преобразовать как подписанный, а затем исправить результат. Есть несколько способов сделать это. (См. Этот раздел вопросов и ответов для некоторых векторизованных вручную методов с внутренними Си-кодами, не все из которых имеют идеально округленные результаты.)
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 = __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
... более смешные прыжки, известные во время компиляции позже...
######## 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
Для сравнения, clang и gcc также не оптимизируют все это во время компиляции, но они понимают, что им не нужен цикл очистки, и просто делают одно скалярное хранилище или конвертируют после соответствующих циклов. (на самом деле лязг полностью разворачивает все, если только вы не скажете этого)
Смотрите код в проводнике компилятора Godbolt.
gcc просто конвертирует верхнюю и нижнюю половину 16b в плавающее отдельно и объединяет их с умножением на 65536 и сложением.
Clang-х unsigned -> float Стратегия конверсии интересна: она никогда не использует cvt инструкция на всех. Я думаю, что он вставляет две 16-разрядные половины целого числа без знака в мантиссу двух чисел напрямую (с некоторыми приемами, чтобы установить показатели (побитовое булево и ADDPS)), а затем складывает половину младшего и старшего разрядов, как это делает gcc.
Конечно, если вы компилируете в 64-битный код, скалярное преобразование может просто увеличить нулю uint32_t в 64-битный и преобразовать это как подписанное int64_t, чтобы плавать. Подписанное int64_t может представлять каждое значение uint32_t, а x86 может эффективно преобразовать 64-разрядное подписанное int в float. Но это не векторизация.
32-разрядный двоичный код IEEE-754, такой как MSVC++, обеспечивает точность только 6-7 десятичных цифр. Ваше начальное значение находится в пределах диапазона этого типа, но кажется, что оно не может быть точно представлено этим типом, как на самом деле имеет место для большинства значений типа uint32_t,
В то же время модуль с плавающей запятой процессора x86 или x86_64 использует более широкое представление, чем 64-битный MSVC++ double, Вероятно, что после выхода из цикла последний вычисленный элемент массива остается в регистре FPU в его расширенной форме точности. Затем программа может использовать это значение непосредственно из регистра, а не считывать его обратно из памяти, что она обязана делать с предыдущими элементами.
Если программа выполняет == Сравнение путем продвижения более узкого представления к более широкому, а не наоборот, тогда эти два значения могут действительно сравниваться неравно, как обратное путешествие от расширенной точности к float и обратно теряет точность. В любом случае оба значения преобразуются в тип double когда перешел к printf(); если они действительно сравниваются неравно, то, вероятно, результаты этих преобразований также различаются.
Я не разбираюсь в опциях компиляции MSVC++, но, скорее всего, есть такой, который бы отменил это поведение. Такие параметры иногда называются "строгая математика" или "строгий фп". Имейте в виду, однако, что включение такой опции (или отключение ее противоположности) может быть очень дорогостоящим в программе с большим количеством FP.