Да, деление двух многочленов часто может дать вам лучший компромисс между скоростью и точностью, чем один огромный многочлен. Пока есть достаточно работы, чтобы скрыть пропускную способность 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) вместо 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, как это делает скалярная версия. (Потому что вас не волнует, что произойдет, если двойные числа слишком велики, чтобы поместиться в целое число.)

Инструкции 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.

Связанный:

• Самая быстрая реализация экспоненциальной функции с использованием AVX и Самая быстрая реализация экспоненциальной функции с использованием SSE имеет версии с другими компромиссами между скоростью и точностью, для _ps (упакованный одинарная точность float).

• Быстрая экспоненциальная SSE низкой точности с использованием операций двойной точности на другом конце спектра, но все же для double.

• Сколько тактов стоит ли экспоненциация AVX/SSE на современном процессоре x86_64? обсуждает некоторые существующие библиотеки, такие как SVML и VCL Agner Fog (под лицензией GPL). И libmvec glibc.

Тогда, конечно, мне нужно проверить все, так как я не совсем уверен, что преобразовал его правильно.

перебор всех 2^64 битовых шаблонов double нецелесообразен, в отличие от float, где их всего 4 миллиарда, но, возможно, перебор всех double, у которых младшие 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);

Как правило, вы можете векторизовать простую настройку условия значения, используя сравнение без ветвления + смешивание. Вместо if(x) y=0; у вас есть SIMD-эквивалент y = (condition) ? 0 : y; для каждого элемента. Сравнение SIMD создает маску из нулевых/всех элементов, чтобы вы могли использовать ее для смешивания.

например в этом случае cmppd ввод и blendvpd вывод, если у вас SSE4.1. Или просто с SSE2, и/и не/или блендить. См. встроенные параметры SSE для сравнения (_mm_cmpeq_ps) и операции присваивания для _ps версии обоих, _pd идентично.

На ассемблере это будет выглядеть так:

; 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 представлен битовым шаблоном со всеми нулями, поэтому сделайте сравнение, которое даст true, если оно находится в диапазоне, и И результат с этим. (Чтобы оставить его без изменений или заставить его +0.0). Это лучше, чем _mm_blendv_pd, который стоит 2 мкп на большинстве процессоров Intel (а 128-битная версия AVX всегда стоит 2 мкп на Intel). И не хуже на AMD или Skylake.

+-Inf представлен битовым шаблоном мантиссы=0, показатель степени=все единицы. (Любое другое значение в мантиссе представляет +-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 с результатом сравнения и + Infinity.

Использование целочисленного арифметического сдвига для трансляции бита знака в каждую позицию в 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 и неотрицательный в противном случае.

Агнер Фог сообщает, что vfixupimmpd находится всего в 1 моп на Skylake-X.