Фактическое деление во время выполнения (FP и целое число) очень медленное, поэтому вы определенно хотите этого избежать. То, как вы написали это выражение, возможно, компилируется для включения деления, потому что математика FP не ассоциативна (без -ffast-math);компилятор не может превратить x / foo * bar в x * (bar/foo) для вас, хотя это очень хорошо с инвариантом цикла bar/foo. Вам нужны либо числа с плавающей точкой, либо 64-разрядные целые числа, чтобы избежать переполнения при умножении, но только FP позволяет повторно использовать нецелочисленный результат с циклическим инвариантным делением.
_mm256_fmadd_ps выглядит как очевидный путьс предварительно вычисленным инвариантным значением цикла для множителя (d - c) / (b - a). Если float округление не является проблемой для того, чтобы сделать это строго по порядку (умножить, а затем разделить), то, вероятно, хорошо сначала выполнить это неточное деление вне цикла. Как
_mm256_set1_ps((d - c) / (double)(b - a)). Использование double для этого вычисления позволяет избежать ошибки округления при преобразовании в FP операндов деления.
Вы используете те же самые a, b, c, d для многих x, предположительно поступающих из смежной памяти. Вы используете результат как часть адреса памяти, поэтому, к сожалению, вам в конечном итоге понадобятся результаты из SIMD в целочисленные регистры. (Возможно, с разбросанными хранилищами AVX512 вы могли бы избежать этого.)
Современные процессоры x86 имеют пропускную способность 2 / тактовой загрузки, поэтому, вероятно, лучшим вариантом для возврата 8x uint32_t в целочисленные регистры является векторная загрузка / целочисленная перезагрузка вместопотратив 2 мопа за элемент для АЛУ перемешать вещи. Это имеет некоторую задержку, поэтому я бы предложил преобразовать в буфер tmp, возможно, 16 или 32 дюйма (64 или 128 байт), то есть 2x или 4x __m256i перед циклом через этот скаляр.
Или, возможно, альтернативное преобразованиеи сохранение одного вектора с последующим циклом по 8 элементам другого , который вы конвертировали ранее. т.е. программный конвейер. Внеочередное выполнение может скрыть задержку, но вы уже собираетесь расширить возможности скрытия задержек для пропусков кэша для всего, что вы делаете с памятью.
В зависимости от вашего процессора (например, Haswell или какой-то другой)Skylake), использование 256-битных векторных инструкций может ограничить ваш максимальный турбо немного ниже, чем в противном случае. Вы можете подумать о том, чтобы делать только векторы по 4, но тогда вы тратите больше мопов на элемент.
Если не SIMD, то даже скалярный C ++ fma() все еще хорош, для vfmadd213sd, но с использованием встроенных функцийочень удобный способ получить округление (вместо усечения) из float -> int (vcvtps2dq, а не vcvttps2dq).
Обратите внимание, что uint32_t <-> float преобразованиене доступен напрямую до AVX512. Для скаляра вы можете просто преобразовать в / из int64_t с усечением / расширением нуля для младшей половины без знака.
Очень удобно, что (как обсуждалось в комментариях) ваши входные данные ограничены по диапазону, так что если вы интерпретируете их какцелые числа со знаком имеют одинаковое значение (неотрицательное число со знаком). Как x, так и x-a (и b-a), как известно, являются положительными и <= INT32_MAX, т.е. <code>0x7FFFFFFF. (Или, по крайней мере, не отрицательно. Ноль в порядке.)
Округление с плавающей запятой
Для SIMD одинарная точность float равна очень хорошо для пропускной способности SIMD. Эффективное упакованное преобразование в / из подписанного int32_t. Но не каждый int32_t может быть точно представлен как float. Большие значения округляются до ближайшего четного, ближайшего кратного 2 ^ 2, 2 ^ 3 или более, чем дальше 2 ^ 24, то значение равно.
Использование SIMD double возможно, но требует некоторого тасования.
Я не думаю, что float обычно является проблемой для формулы, написанной с (float)(x-a). Если входной диапазон b-a большой, это означает, что оба диапазона велики, и ошибка округления не приведет к отображению всех возможных значений x в один и тот же выход. В зависимости от множителя ошибка округления на входе может быть хуже, чем ошибка округления на выходе, возможно, оставляя некоторые представимые выходные значения с плавающей запятой неиспользованными для более высоких значений x-a.
Но если мы хотим выделить часть -a * (d - c) / (b - a)и объедините его с +c в конце, затем
- У нас потенциально может быть потеря точности из-за катастрофической отмены в этом значении.
- Нам нужно сделать
(float)x для необработанного входного значения. Если a огромен, а b-a мал, т. Е. Небольшой диапазон в верхней части возможного входного диапазона, ошибка округления может отображать все возможные значения x в один и тот же тип с плавающей запятой. - Для лучшегоПри использовании FMA мы хотим сделать
+c перед преобразованием обратно в целое число, что снова может привести к ошибке округления выходных данных, если d-c - небольшой выходной диапазон, но c огромен. В твоем случае не проблема;с d <= 2 ^ <sup>20 - 1 мы знаем, что float может точно представлять каждое выходное целое значение в этом диапазоне c..d.
Если у вас не было ограничения входного диапазона, вы могли бы сместить диапазон до / от знака до масштабирования, используя целое число (x-a)+0x80000000U на входе и ...+c+0x80000000U на выходе (после округления до ближайшего int32_t). Но это привело бы к огромной float ошибке округления для небольших uint32_t входов (близких к 0), которые смещены по диапазону, чтобы приблизиться к INT_MIN.
Нам не нужно сдвигать диапазон дляb-a или d-c, потому что + или - или XOR с 0x80000000U отменят в вычитаниях.
Пример:
Векторы const должны быть поднятыпосле этого он может выйти из цикла, или вы можете сделать это вручную.
Для этого требуется AVX1 + FMA (например, AMD Piledriver или Intel Haswell или более поздняя версия). Не проверено, извините, я даже не бросил это на Godbolt, чтобы посмотреть, скомпилируется ли он.
// fastest but not safe if b-a is small and a > 2^24
static inline
__m256i range_scale_fast_fma(__m256i data, uint32_t a, uint32_t b, uint32_t c, uint32_t d)
{
// avoid rounding errors when computing the scale factor, but convert double->float on the final result
double scale_scalar = (d - c) / (double)(b - a);
const __m256 scale = _mm256_set1_ps(scale_scalar);
const __m256 add = _m256_set1_ps(-a*scale_scalar + c);
// (x-a) * scale + c
// = x * scale + (-a*scale + c) but with different rounding error from doing -a*scale + c
__m256 in = _mm256_cvtepi32_ps(data);
__m256 out = _mm256_fmadd_ps(in, scale, add);
return _mm256_cvtps_epi32(out); // convert back with round to nearest-even
// _mm256_cvttps_epi32 truncates, matching C rounding; maybe good for scalar testing
}
Или более безопасная версия, с помощью которой вводится диапазонное смещение с целым числом: Вы можете легкопри необходимости избегайте FMA для переносимости (просто AVX1) и используйте целочисленное сложение для вывода. Но мы знаем, что выходной диапазон достаточно мал, чтобы он всегда мог точно представлять любое целое число
static inline
__m256i range_scale_safe_fma(__m256i data, uint32_t a, uint32_t b, uint32_t c, uint32_t d)
{
// avoid rounding errors when computing the scale factor, but convert double->float on the final result
const __m256 scale = _mm256_set1_ps((d - c) / (double)(b - a));
const __m256 cvec = _m256_set1_ps(c);
__m256i in_offset = _mm256_add_epi32(data, _mm256_set1_epi32(-a)); // add can more easily fold a load of a memory operand than sub because it's commutative. Only some compilers will do this for you.
__m256 in_fp = _mm256_cvtepi32_ps(in_offset);
__m256 out = _mm256_fmadd_ps(in_fp, scale, _mm256_set1_ps(c)); // in*scale + c
return _mm256_cvtps_epi32(out);
}
Без FMA вы все равно могли бы использовать vmulps. Вы можете также преобразовать обратно в целое число перед добавлением c, если вы делаете это, хотя vaddps будет безопасным.
Вы можете использовать это в цикле, например
void foo(uint32_t *arr, ptrdiff_t len)
{
if (len < 24) special case;
alignas(32) uint32_t tmpbuf[16];
// peel half of first iteration for software pipelining / loop rotation
__m256i arrdata = _mm256_loadu_si256((const __m256i*)&arr[0]);
__m256i outrange = range_scale_safe_fma(arrdata);
_mm256_store_si256((__m256i*)tmpbuf, outrange);
// could have used an unsigned loop counter
// since we probably just need an if() special case handler anyway for small len which could give len-23 < 0
for (ptrdiff_t i = 0 ; i < len-(15+8) ; i+=16 ) {
// prep next 8 elements
arrdata = _mm256_loadu_si256((const __m256i*)&arr[i+8]);
outrange = range_scale_safe_fma(arrdata);
_mm256_store_si256((__m256i*)&tmpbuf[8], outrange);
// use first 8 elements
for (int j=0 ; j<8 ; j++) {
use tmpbuf[j] which corresponds to arr[i+j]
}
// prep 8 more for next iteration
arrdata = _mm256_loadu_si256((const __m256i*)&arr[i+16]);
outrange = range_scale_safe_fma(arrdata);
_mm256_store_si256((__m256i*)&tmpbuf[0], outrange);
// use 2nd 8 elements
for (int j=8 ; j<16 ; j++) {
use tmpbuf[j] which corresponds to arr[i+j]
}
}
// use tmpbuf[0..7]
// then cleanup: one vector at a time until < 8 or < 4 with 128-bit vectors, then scalar
}
Эти имена переменных звучат глупо, но я не мог придумать ничего лучшего.
Этот программный конвейер является оптимизацией;Вы можете просто заставить его работать / опробовать его с одним вектором за один раз, используемый сразу. (Оптимизируйте перезагрузку первого элемента с перезагрузки до vmovd, используя _mm_cvtsi128_si32(_mm256_castsi256_si128(outrange)), если хотите.)
Особые случаи
Если есть случаи, когда вы знаете (b - a)это степень 2, вы можете отсканировать биты с помощью tzcnt или bsf, а затем умножить. (Для них есть, например, GNU C __builtin_ctz() для подсчета конечных нулей.)
Или вы можете гарантировать , что (b - a) всегда является степенью 2?
Или лучше, , если (b - a) / (d - c) - это точная степень 2, все это можно просто сдвинуть / сдвинуть вправо / добавить.
Если вы не можете всегда убедитесь, что вам все еще иногда требуется общий случай, но, возможно, возможно сделать это эффективно.