Да, деление двух полиномов часто может дать вам лучший компромисс между скоростью и точностью, чем один огромный полином. Пока есть достаточно работы, чтобы скрыть пропускную способность divpd. (Последние процессоры x86 имеют довольно приличную пропускную способность деления FP. Все еще плохо против умножения, но это всего 1 моп, поэтому он не останавливает конвейер, если вы используете его достаточно редко, то есть смешиваете с большим количеством умножений. Включая в окружающий код использует exp)
Однако _mm_cvtepi64_pd(_mm_cvttpd_epi64(px)); не будет работать с SSE2. Для встроенных преобразований в / из 64-разрядных целых чисел требуется AVX512DQ .
Чтобы выполнить округленное округление до ближайшего целого числа, в идеале вы должны использовать SSE4.1 _mm_round_pd(x, _MM_FROUND_TO_NEAREST_INT |_MM_FROUND_NO_EXC), (или усечение в направлении нуля, или пол или потолок в направлении - + Inf).
Но нам на самом деле это не нужно.
Скалярный код заканчивается на int n и double px, представляющих одно и то же числовое значение. Он использует bad / buggy floor(val+0.5) idiom вместо rint(val) или nearbyint(val) для округления до ближайшего, а затем преобразует это уже целое число double в int (с семантикой усечения в C ++ , но это не имеет значения, потому что значение double уже является точным целым числом.)
При использовании SIMD, проще всего просто конвертировать в 32-битное целое и обратно.
__m128i n = _mm_cvtpd_epi32( _mm_mul_pd(log2e, x) ); // round to nearest
__m128d px = _mm_cvtepi32_pd( n );
Округление до int с требуемым режимом, затем преобразование обратно в double, эквивалентно double-> double округлению, а затем получение версии int, как это делает скалярная версия. (Потому что вас не волнует, что происходит для двойников, слишком больших, чтобы поместиться в int.)
Команды cvtsd2si и si2sd по 2 мопа каждая и перемешивают 32-битные целые числа, упакованные в младшие 64 бита вектора. Поэтому, чтобы установить 64-битные целочисленные сдвиги, чтобы снова вставить биты в double, вам нужно будет перемешать. Старшие 64 бита n будут нулями, поэтому мы можем использовать это для создания 64-битного целого числа n, выровненного с двойными числами:
n = _mm_shuffle_epi32(n, _MM_SHUFFLE(3,1,2,0)); // 64-bit integers
Но только с SSE2 есть обходные пути. Преобразование в 32-разрядное целое число и обратно является одним из вариантов: вам не нужны входные данные, слишком маленькие или слишком большие. Но пакетное преобразование между double и int стоит как минимум 2 моп на процессорах Intel в каждом направлении, то есть всего 4. Но только 2 из этих мопов нуждаются в блоках FMA, и ваш код, вероятно, не является узким местом на порту 5 со всеми этими умножениями и сложениями.
Или добавьте очень большое число и вычтите его снова: достаточно большое, чтобы каждый double был на расстоянии 1 целое, поэтому обычное округление FP делает то, что вы хотите. (Это работает для входных данных, которые не умещаются в 32 бита, но не double> 2 ^ 52. Так что в любом случае это будет работать.) Также см. Как эффективно выполнять преобразования double / int64 с SSE / AVX? который использует этот трюк. Хотя я не смог найти пример на SO.
Похожие:
Тогда, конечно, мне нужно проверить все, так как я не совсем уверен, что преобразовал это правильно.
итерация по всем 2 ^ 64 double битовым шаблонам нецелесообразна, в отличие от float, где существует только 4 миллиарда, но, возможно, итерация по всем double s, у которых младшие 32 бита их мантиссы равны нулю было бы хорошим началом. то есть проверка в цикле с
bitpatterns = _mm_add_epi64(bitpatterns, _mm_set1_epi64x( 1ULL << 32 ));
doubles = _mm_castsi128_pd(bitpatterns);
https://randomascii.wordpress.com/2014/01/27/theres-only-four-billion-floatsso-test-them-all/
Для тех нескольких последних строк, корректирующих ввод для входов вне диапазона:
Версия float, которую вы цитируете, полностью исключает проверку диапазона.Это, очевидно, самый быстрый способ , если ваши входы всегда будут в диапазоне или если вас не волнует, что происходит с входами вне диапазона.
Альтернативная более дешевая проверка диапазона (возможно, толькодля отладки) будет превращать значения вне диапазона в NaN, ИЛИ результат упакованного сравнения в результат.(Битовый массив «все единицы» представляет NaN.)
__m128d out_of_bounds = _mm_cmplt_pd( limit, abs(initial_x) ); // abs = mask off the sign bit
result = _mm_or_pd(result, out_of_bounds);
В общем случае вы можете векторизовать простую установку условия для значения, используя сравнение без ветвлений + blend .Вместо if(x) y=0; у вас есть SIMD-эквивалент y = (condition) ? 0 : y; для каждого элемента. SIMD сравнивает, создает маску из элементов «все ноль / все-один», так что вы можете использовать ее для смешивания.
Например, в этом случае cmppd для ввода и blendvpd для вывода, если у вас SSE4.1.Или только с SSE2, и / или без / или слияния.См. Встроенные функции SSE для сравнения (_mm_cmpeq_ps) и операции назначения для _ps версии обоих, _pd идентичен.
В asm это будет выглядеть так:
; result in xmm0 (in need of fixups for out of range inputs)
; initial_x in xmm2
; constants:
; xmm5 = limit
; xmm6 = +Inf
cmpltpd xmm2, xmm5 ; xmm2 = input_x < limit ? 0xffff... : 0
andpd xmm0, xmm2 ; result = result or 0
andnpd xmm2, xmm6 ; xmm2 = 0 or +Inf (In that order because we used ANDN)
orpd xmm0, xmm2 ; result |= 0 or +Inf
; xmm0 = (input < limit) ? result : +Inf
(В более ранней версии ответа я думал, что, возможно, я сохраняю movaps для копирования регистра, но это просто стандартная смесь для болота. Она уничтожает initial_x, поэтому компилятор должен скопировать этозарегистрируйтесь в какой-то момент при расчете result, хотя.)
Оптимизации для этого специального условия
Или в этом случае 0.0представлен битовым шаблоном со всеми нулями , так что сделайте сравнение, которое даст истинное значение, если оно находится в пределах диапазона, и И с ним результат.(Чтобы оставить его без изменений или принудительно +0.0).Это лучше, чем _mm_blendv_pd, который стоит 2 моп на большинстве процессоров Intel (а 128-разрядная версия AVX всегда стоит 2 моп на Intel).И это не хуже на AMD или Skylake.
+-Inf представлен битовой комбинацией значенийи = 0, экспонента = все-единицы.(Любое другое значение в значенииand представляет + -NaN.) Поскольку слишком большие входные данные, по-видимому, будут по-прежнему оставлять ненулевые значения, мы не можем просто сравнивать результат И и ИЛИ с конечным результатом.Я думаю, нам нужно сделать обычное смешивание или что-то более дорогое (3 мопа и векторная константа).
Это добавляет 2 цикла задержки к конечному результату;И ANDNPD, и ORPD находятся на критическом пути.CMPPD и ANDPD не являются;они могут работать параллельно с тем, что вы делаете для вычисления результата.
Надеюсь, ваш компилятор будет использовать ANDPS и т. д., а не PD, для всего, кроме CMP, потому что он на 1 байт короче, но идентичен, потому что они 'оба поразрядных операций.Я написал ANDPD просто для того, чтобы мне не приходилось объяснять это в комментариях.
Возможно, вы сможете сократить задержку критического пути, комбинируя обе исправления перед применением к результату Таким образом, у вас есть только одна смесь.Но тогда я думаю, что вам также нужно объединить результаты сравнения.
Или, поскольку ваши верхняя и нижняя границы имеют одинаковую величину, может быть, вы можете сравнить абсолютное значение?(замаскируйте бит знака initial_x и сделайте _mm_cmplt_pd(abs_initial_x, _mm_set1_pd(details::EXP_LIMIT))).Но тогда вам нужно разобраться, обнулить ли или установить + Inf.
Если бы у вас был SSE4.1 для _mm_blendv_pd, вы могли бы использовать initial_x в качестве элемента управления смешиванием для исправления, которое может потребоватьсяприменяется, потому что blendv заботится только о знаковом бите управления смешиванием (в отличие от версии AND / ANDN / OR, где все биты должны совпадать.)
__m128d fixup = _mm_blendv_pd( _mm_setzero_pd(), _mm_set1_pd(INFINITY), initial_x ); // fixup = (initial_x signbit) ? 0 : +Inf
// see below for generating fixup with an SSE2 integer arithmetic-shift
const signbit_mask = _mm_castsi128_pd(_mm_set1_epi64x(0x7fffffffffffffff)); // ~ set1(-0.0)
__m128d abs_init_x = _mm_and_pd( initial_x, signbit_mask );
__m128d out_of_range = _mm_cmpgt_pd(abs_init_x, details::EXP_LIMIT);
// Conditionally apply the fixup to result
result = _mm_blendv_pd(result, fixup, out_of_range);
Возможно использовать cmplt вместо cmpgt и переставьте, если вам небезразлично, что произойдет, если initial_x будет NaN .Выбор сравнения так, чтобы false применял исправление вместо true, будет означать, что неупорядоченное сравнение приводит к 0 или + Inf для ввода -NaN или + NaN.Это все еще не делает распространение NaN.Вы можете _mm_cmpunord_pd(initial_x, initial_x) и ИЛИ это в fixup, если хотите, чтобы это произошло.
Особенно на Skylake и AMD Bulldozer / Ryzen, где SSE2 blendvpd всего 1 моп, это должно быть довольно неплохо.(Код VEX, vblendvpd равен 2 моп, с 3 входами и отдельным выходом.)
Вы все еще можете использовать некоторые из этих идей только с SSE2, возможно, создав fixup, выполнивсравнить с нулем, а затем _mm_and_pd или _mm_andnot_pd с результатом сравнения и + Бесконечность.
Использование целочисленного арифметического сдвига для широковещательной передачи бита знака на каждую позицию в double неэффективный: psraq не существует, только psraw/d.Только логические сдвиги имеют размер 64-битного элемента.
Но вы могли бы создать fixup с помощью всего одного целочисленного сдвига и маски и побитового инвертирования
__m128i ix = _mm_castsi128_pd(initial_x);
__m128i ifixup = _mm_srai_epi32(ix, 11); // all 11 bits of exponent field = sign bit
ifixup = _mm_and_si128(ifixup, _mm_set1_epi64x(0x7FF0000000000000ULL) ); // clear other bits
// ix = the bit pattern for 0 (non-negative x) or +Inf (negative x)
__m128d fixup = _mm_xor_si128(ifixup, _mm_set1_epi32(-1)); // bitwise invert
Затем смешайте fixup с result для входов вне диапазона как обычно.
Недорого проверяйте abs(initial_x) > details::EXP_LIMIT
Если алгоритм expбыл уже в квадрате initial_x, вы можете сравнить с EXP_LIMIT в квадрате.Но это не так, xx = x*x происходит только после некоторого вычисления, чтобы создать x.
Если у вас AVX512F / VL, VFIXUPIMMPD может пригодиться здесь.Он предназначен для функций, в которых выходные данные для особых случаев поступают от «специальных» входов, таких как NaN и + -Inf, отрицательных, положительных или нулевых значений, сохраняя сравнение для этих случаев.(например, после обратной Ньютона-Рафсона (x) для x = 0.)
Но оба ваших особых случая нуждаются в сравнении.Или они?
Если вы возводите в квадрат свои входные данные и вычитаете, это стоит только одного FMA, чтобы сделать initial_x * initial_x - details::EXP_LIMIT * details::EXP_LIMIT, чтобы создать результат, отрицательный для abs(initial_x) < details::EXP_LIMIT, и неотрицательный в противном случае.
Agner Fog сообщает, что vfixupimmpd - это всего 1 моп на Skylake-X.