Эффективная целочисленная функция пола в С++

Ответ 1

Посмотрите на магические числа. Алгоритм, предложенный на веб-странице, должен быть гораздо более эффективным, чем простое приведение. Я никогда не использовал его сам, но вот сравнение производительности, которое они предлагают на сайте (xs_ToInt и xs_CRoundToInt - предлагаемые функции):

Performing 10000000 times:
simple cast           2819 ms i.e. i = (long)f;
xs_ToInt              1242 ms i.e. i = xs_ToInt(f); //numerically same as above
bit-twiddle(full)     1093 ms i.e. i = BitConvertToInt(f); //rounding from Fluid
fistp                  676 ms i.e. i = FISTToInt(f); //Herf, et al x86 Assembly rounding 
bit-twiddle(limited)   623 ms i.e. i = FloatTo23Bits(f); //Herf, rounding only in the range (0...1]  
xs_CRoundToInt         609 ms i.e. i = xs_CRoundToInt(f); //rounding with "magic" numbers

Кроме того, xs_ToInt, видимо, модифицируется, чтобы улучшить производительность:

Performing 10000000 times:
simple cast convert   3186 ms i.e. fi = (f*65536);
fistp convert         3031 ms i.e. fi = FISTToInt(f*65536);
xs_ToFix               622 ms i.e. fi = xs_Fix<16>::ToFix(f);

Краткое объяснение того, как работает метод "магических чисел":

"По сути, чтобы добавить два числа с плавающей запятой, ваш процессор" выстраивает "десятичные точки чисел так, чтобы он мог легко добавлять биты. Он делает это путем" нормализации "чисел таким образом, чтобы сохранялись старшие значащие биты. т.е. меньшее число "нормализуется", чтобы соответствовать большему. Таким образом, принцип преобразования "магического числа", которое использует xs_CRoundToInt(), заключается в следующем: мы добавляем достаточно большое число с плавающей запятой (число, которое настолько велико, что существует значащие цифры только до десятичной точки и ни одной после нее) до той, которую вы преобразуете так, что: (a) число нормализуется процессором до его целочисленного эквивалента и (b) добавление двух не стирает интеграл значащие биты в числе, которое вы пытались преобразовать (т.е. XX00 + 00YY = XXYY). "

Цитата взята с той же веб-страницы.

Ответ 2

Преобразование в int заведомо медленное.

Может быть, вы жили в безвыходном положении с x86-64 или иначе упустили из виду, что это не было верно в течение некоторого времени на x86. :)

SSE/SSE2 имеют инструкцию для преобразования с усечением (вместо режима округления по умолчанию). ISA эффективно поддерживает эту операцию именно потому, что преобразование с семантикой C не редкость в реальных кодовых базах. Код x86-64 использует регистры XMM SSE/SSE2 для скалярной математики FP, а не x87, из-за этого и других факторов, делающих его более эффективным. Даже современный 32-битный код использует XMM-регистры для скалярной математики.

При компиляции для x87 (без SSE3 fisttp) компиляторам приходилось менять режим округления x87 на усечение, сохранять FP в памяти, а затем снова изменять режим округления. (А затем перезагрузите целое число из памяти, как правило, из локального стека, если с ним что-то делать.) X87 был ужасен для этого.

Да, это было ужасно медленно, например, в 2006 году, когда была написана ссылка в ответе @Kirjain, если у вас все еще был 32-разрядный процессор или вы использовали процессор x86-64 для запуска 32-разрядного кода.

Конвертация с использованием режима округления, отличного от усечения или по умолчанию (ближайшего), напрямую не поддерживается, и до тех пор, пока в roundps/roundpd SSE4.1 не roundps roundpd всего использовать трюки с магическим числом, как в ссылке 2006 года из ответа @Kirjain.

Там есть несколько приятных трюков, но только для double → 32-разрядного целого числа. Вряд ли стоит расширяться в double если у вас есть float.

Или, как правило, просто добавьте число большой величины для запуска округления, а затем вычтите его снова, чтобы вернуться к исходному диапазону. Это может работать для float без расширения до double, но я не уверен, насколько легко сделать работу floor.

В любом случае, очевидным решением здесь являются _mm256_floor_ps() и _mm256_cvtps_epi32 (vroundps и vcvtps2dq). Не-AVX версия этого может работать с SSE4.1.

Я не уверен, сможем ли мы сделать еще лучше; Если у вас был огромный массив для обработки (и вам не удавалось чередовать эту работу с другой работой), вы могли бы установить режим округления MXCSR на "в сторону -Inf" (этаж) и просто использовать vcvtps2dq (который использует текущее округление Режим). Затем установите его обратно. Но, вероятно, лучше кешировать блокировку вашего преобразования или делать это на лету, когда вы генерируете данные, предположительно из других вычислений FP, для которых режим округления FP должен быть установлен в значение по умолчанию Ближайший.

roundps/pd/ss/sd - 2 моп на процессорах Intel, но только 1 моп (на 128-битную дорожку) на AMD Ryzen. cvtps2dq также 1 моп. Пакетное double-> int преобразование также включает случайное перемешивание. Скалярное преобразование FP-> int (которое копирует в регистр целых чисел) обычно также стоит для этого дополнительного мопа.

Так что в некоторых случаях есть вероятность того, что трюки с магическими числами будут выигрышем; возможно, стоит разобраться, являются ли _mm256_floor_ps() + cvt частью критического узкого места (или, более вероятно, если у вас есть double и вы хотите int32).

@Cássio Renan int foo = floorf(f) будет фактически векторизоваться, если скомпилирован с помощью gcc -O3 -fno-trapping-math (или -ffast-math), с -march= чем-то, что имеет SSE4.1 или AVX. https://godbolt.org/z/ae_KPv

Это может быть полезно, если вы используете это с другим скалярным кодом, который не векторизован вручную. Особенно, если вы надеетесь, что компилятор будет автоматически векторизовать все это.

Ответ 3

Если вы делаете это в пакетном режиме, компилятор может автоматически векторизовать его, если вы знаете, что делаете. Например, вот небольшая реализация, которая автоматически векторизует преобразование чисел с плавающей точкой в GCC:

#include <cmath>

// Compile with -O3 and -march=native to see autovectorization
__attribute__((optimize("-fno-trapping-math")))
void testFunction(float* input, int* output, int length) {
  // Assume the input and output are aligned on a 32-bit boundary.
  // Of course, you have  to ensure this when calling testFunction, or else
  // you will have problems.
  input = static_cast<float*>(__builtin_assume_aligned(input, 32));
  output = static_cast<int*>(__builtin_assume_aligned(output, 32));

  // Also assume the length is a multiple of 32.
  if (length & 31) __builtin_unreachable();

  // Do the conversion
  for (int i = 0; i < length; ++i) {
    output[i] = floor(input[i]);
  }
}

Это сгенерированная сборка для x86-64 (с инструкциями AVX512):

testFunction(float*, int*, int):
        test    edx, edx
        jle     .L5
        lea     ecx, [rdx-1]
        xor     eax, eax
.L3:
        # you can see here that the conversion was vectorized
        # to a vrndscaleps (that will round the float appropriately)
        # and a vcvttps2dq (thal will perform the conversion)
        vrndscaleps     ymm0, YMMWORD PTR [rdi+rax], 1
        vcvttps2dq      ymm0, ymm0
        vmovdqa64       YMMWORD PTR [rsi+rax], ymm0
        add     rax, 32
        cmp     rax, rdx
        jne     .L3
        vzeroupper
.L5:
        ret

Если ваша цель не поддерживает AVX512, она все равно будет автоматически выполнять векторизацию с использованием инструкций SSE4.1, если они у вас есть. Это вывод с -O3 -msse4.1:

testFunction(float*, int*, int):
        test    edx, edx
        jle     .L1
        shr     edx, 2
        xor     eax, eax
        sal     rdx, 4
.L3:
        roundps xmm0, XMMWORD PTR [rdi+rax], 1
        cvttps2dq       xmm0, xmm0
        movaps  XMMWORD PTR [rsi+rax], xmm0
        add     rax, 16
        cmp     rax, rdx
        jne     .L3
.L1:
        ret

Посмотрите это в прямом эфире на Godbolt

Ответ 4

Почему бы просто не использовать это:

#include <cmath>

auto floor_(float const x) noexcept
{
  int const t(x);

  return t - (t > x);
}

Ответ 5

Вот модификация Кассио Ренанс отличный ответ. Он заменяет все специфичные для компилятора расширения стандартным C++ и теоретически переносим на любой соответствующий компилятор. Кроме того, он проверяет, что аргументы правильно выровнены, а не предполагают это. Это оптимизирует к тому же коду.

#include <assert.h>
#include <cmath>
#include <stddef.h>
#include <stdint.h>

#define ALIGNMENT alignof(max_align_t)
using std::floor;

// Compiled with: -std=c++17 -Wall -Wextra -Wpedantic -Wconversion -fno-trapping-math -O -march=cannonlake -mprefer-vector-width=512

void testFunction(const float in[], int32_t out[], const ptrdiff_t length)
{
  static_assert(sizeof(float) == sizeof(int32_t), "");
  assert((uintptr_t)(void*)in % ALIGNMENT == 0);
  assert((uintptr_t)(void*)out % ALIGNMENT == 0);
  assert((size_t)length % (ALIGNMENT/sizeof(int32_t)) == 0);

  alignas(ALIGNMENT) const float* const input = in;
  alignas(ALIGNMENT) int32_t* const output = out;

  // Do the conversion
  for (int i = 0; i < length; ++i) {
    output[i] = static_cast<int32_t>(floor(input[i]));
  }
}

Это не так хорошо оптимизирует GCC, как оригинал, который использовал непереносимые расширения. Стандарт C++ поддерживает спецификатор alignas, ссылки на выровненные массивы и функцию std::align которая возвращает выровненный диапазон в буфере. Однако ни один из них не позволяет любому протестированному компилятору генерировать выровненные, а не выровненные векторные нагрузки и хранилища.

Хотя alignof(max_align_t) составляет только 16 для x86_64, и можно определить ALIGNMENT как константу 64, это не поможет любому компилятору генерировать лучший код, поэтому я остановился на переносимости. Наиболее близким к переносимому способу заставить компилятор предположить, что poitner выровнен, будет использование типов из <immintrin.h>, которые большинство компиляторов для x86 поддерживают, или определение struct с помощью alignas. Проверяя предопределенные макросы, вы также можете расширить макрос до __attribute__ ((aligned (ALIGNMENT))) в компиляторах Linux или __declspec (align (ALIGNMENT)) в компиляторах Windows, и что-то более безопасное для компилятора, о котором мы не знаем, но GCC нуждается в атрибуте типа для фактического создания выровненных загрузок и хранилищ.

Кроме того, в оригинальном примере вызывался скрытый вход, чтобы сообщить GCC, что length не может быть кратна 32. Если вы assert() это или вызываете стандартную функцию, такую как abort(), ни GCC, ни Clang, ни ICC сделает такой же вывод. Таким образом, большая часть кода, который они генерируют, будет обрабатывать случай, когда length не является кратным кратным векторной ширине.

Вероятная причина этого заключается в том, что ни одна оптимизация не дает вам такой большой скорости: инструкции невыровненной памяти с выровненными адресами быстры на процессорах Intel, а код для обработки случая, когда length не является хорошим круглым числом, имеет длину несколько байтов и выполняется в постоянное время

В качестве сноски GCC может оптимизировать встроенные функции из <cmath> лучше, чем макросы, реализованные в <math.c>.

GCC 9.1 требуется определенный набор опций для генерации кода AVX512. По умолчанию, даже с -march=cannonlake, он предпочтет 256-битные векторы. Для генерации 512-битного кода требуется параметр -mprefer-vector-width=512. (Спасибо Peter Cordes за указание на это.) Далее следует векторизованный цикл с развернутым кодом для преобразования любых оставшихся элементов массива.

Вот векторизованный основной цикл, за исключением некоторой инициализации с постоянным временем, проверки ошибок и кода очистки, который будет выполняться только один раз:

.L7:
        vrndscaleps     zmm0, ZMMWORD PTR [rdi+rax], 1
        vcvttps2dq      zmm0, zmm0
        vmovdqu32       ZMMWORD PTR [rsi+rax], zmm0
        add     rax, 64
        cmp     rax, rcx
        jne     .L7

Eagle-eyed заметит два отличия от кода, сгенерированного программой Cássio Renans: он использует% zmm вместо регистров% ymm и сохраняет результаты с выровненным vmovdqu32 а не с выровненным vmovdqa64.

Clang 8.0.0 с одинаковыми флагами делает разные варианты развертывания циклов. Каждая итерация работает с восемью 512-битными векторами (то есть 128 плавающими одинарной точности), но код для сбора остатков не разворачивается. Если после этого осталось не менее 64 операций с плавающей запятой, он использует для них еще четыре инструкции AVX512, а затем очищает все дополнения с помощью невекторизованного цикла.

Если вы скомпилируете исходную программу в Clang++, она примет ее без жалоб, но не будет выполнять те же оптимизации: она по-прежнему не будет предполагать, что length кратна ширине вектора или что указатели выровнены.

Он предпочитает код AVX512 AVX256, даже без -mprefer-vector-width=512.

        test    rdx, rdx
        jle     .LBB0_14
        cmp     rdx, 63
        ja      .LBB0_6
        xor     eax, eax
        jmp     .LBB0_13
.LBB0_6:
        mov     rax, rdx
        and     rax, -64
        lea     r9, [rax - 64]
        mov     r10, r9
        shr     r10, 6
        add     r10, 1
        mov     r8d, r10d
        and     r8d, 1
        test    r9, r9
        je      .LBB0_7
        mov     ecx, 1
        sub     rcx, r10
        lea     r9, [r8 + rcx]
        add     r9, -1
        xor     ecx, ecx
.LBB0_9:                                # =>This Inner Loop Header: Depth=1
        vrndscaleps     zmm0, zmmword ptr [rdi + 4*rcx], 9
        vrndscaleps     zmm1, zmmword ptr [rdi + 4*rcx + 64], 9
        vrndscaleps     zmm2, zmmword ptr [rdi + 4*rcx + 128], 9
        vrndscaleps     zmm3, zmmword ptr [rdi + 4*rcx + 192], 9
        vcvttps2dq      zmm0, zmm0
        vcvttps2dq      zmm1, zmm1
        vcvttps2dq      zmm2, zmm2
        vmovups zmmword ptr [rsi + 4*rcx], zmm0
        vmovups zmmword ptr [rsi + 4*rcx + 64], zmm1
        vmovups zmmword ptr [rsi + 4*rcx + 128], zmm2
        vcvttps2dq      zmm0, zmm3
        vmovups zmmword ptr [rsi + 4*rcx + 192], zmm0
        vrndscaleps     zmm0, zmmword ptr [rdi + 4*rcx + 256], 9
        vrndscaleps     zmm1, zmmword ptr [rdi + 4*rcx + 320], 9
        vrndscaleps     zmm2, zmmword ptr [rdi + 4*rcx + 384], 9
        vrndscaleps     zmm3, zmmword ptr [rdi + 4*rcx + 448], 9
        vcvttps2dq      zmm0, zmm0
        vcvttps2dq      zmm1, zmm1
        vcvttps2dq      zmm2, zmm2
        vcvttps2dq      zmm3, zmm3
        vmovups zmmword ptr [rsi + 4*rcx + 256], zmm0
        vmovups zmmword ptr [rsi + 4*rcx + 320], zmm1
        vmovups zmmword ptr [rsi + 4*rcx + 384], zmm2
        vmovups zmmword ptr [rsi + 4*rcx + 448], zmm3
        sub     rcx, -128
        add     r9, 2
        jne     .LBB0_9
        test    r8, r8
        je      .LBB0_12
.LBB0_11:
        vrndscaleps     zmm0, zmmword ptr [rdi + 4*rcx], 9
        vrndscaleps     zmm1, zmmword ptr [rdi + 4*rcx + 64], 9
        vrndscaleps     zmm2, zmmword ptr [rdi + 4*rcx + 128], 9
        vrndscaleps     zmm3, zmmword ptr [rdi + 4*rcx + 192], 9
        vcvttps2dq      zmm0, zmm0
        vcvttps2dq      zmm1, zmm1
        vcvttps2dq      zmm2, zmm2
        vcvttps2dq      zmm3, zmm3
        vmovups zmmword ptr [rsi + 4*rcx], zmm0
        vmovups zmmword ptr [rsi + 4*rcx + 64], zmm1
        vmovups zmmword ptr [rsi + 4*rcx + 128], zmm2
        vmovups zmmword ptr [rsi + 4*rcx + 192], zmm3
.LBB0_12:
        cmp     rax, rdx
        je      .LBB0_14
.LBB0_13:                               # =>This Inner Loop Header: Depth=1
        vmovss  xmm0, dword ptr [rdi + 4*rax] # xmm0 = mem[0],zero,zero,zero
        vroundss        xmm0, xmm0, xmm0, 9
        vcvttss2si      ecx, xmm0
        mov     dword ptr [rsi + 4*rax], ecx
        add     rax, 1
        cmp     rdx, rax
        jne     .LBB0_13
.LBB0_14:
        pop     rax
        vzeroupper
        ret
.LBB0_7:
        xor     ecx, ecx
        test    r8, r8
        jne     .LBB0_11
        jmp     .LBB0_12

ICC 19 также генерирует инструкции AVX512, но сильно отличается от clang. Он больше настраивается с помощью магических констант, но не разворачивает никаких циклов, работая вместо этого на 512-битных векторах.

Этот код также работает на других компиляторах и архитектурах. (Хотя MSVC поддерживает только ISA вплоть до AVX2 и не может автоматически векторизовать цикл.) Например, в ARM с -march=armv8-a+simd он генерирует векторизованный цикл с frintm v0.4s, v0.4s и fcvtzs v0.4s, v0.4s.

Попробуйте сами.