Почему floor () такой медленный?

Question

Недавно я написал некоторый код (ISO / ANSI C) и был удивлен низкой производительностью, которую он достиг. Короче говоря, оказалось, что виновником была функция floor(). Это было не только медленно, но и не векторизовано (с компилятором Intel, иначе ICL).

Вот некоторые тесты для определения пола для всех ячеек в 2D-матрице:

VC:  0.10
ICL: 0.20

Сравните это с простым составом:

VC:  0.04
ICL: 0.04

Как floor() может быть намного медленнее, чем простой актерский состав ?! По сути, это то же самое (кроме отрицательных чисел). 2-й вопрос: кто-нибудь знает о сверхбыстрой реализации floor()?

PS: Вот цикл, который я тестировал:

void Floor(float *matA, int *intA, const int height, const int width, const int width_aligned)
{
    float *rowA=NULL;
    int   *intRowA=NULL;
    int   row, col;

    for(row=0 ; row<height ; ++row){
        rowA = matA + row*width_aligned;
        intRowA = intA + row*width_aligned;
#pragma ivdep
        for(col=0 ; col<width; ++col){
            /*intRowA[col] = floor(rowA[col]);*/
            intRowA[col] = (int)(rowA[col]);
        }
    }
}

Answer 1

Несколько вещей делают пол медленнее, чем гипс, и предотвращают векторизацию.

Самый важный:

Этаж может изменить глобальное состояние. Если вы передаете значение, которое слишком велико, чтобы представлять его как целое число в формате с плавающей запятой, переменная errno будет установлена ​​на EDOM . Специальная обработка для NaN также выполняется. Все это относится к приложениям, которые хотят обнаружить случай переполнения и как-то справиться с ситуацией (не спрашивайте меня, как).

Обнаружение этих проблемных условий не простое и составляет более 90% времени исполнения пола. Фактическое округление дешево и может быть встроено / векторизовано. Также много кода, так что встраивание всей функции floor сделает вашу программу более медленной.

Некоторые компиляторы имеют специальные флаги компилятора, которые позволяют компилятору оптимизировать некоторые из редко используемых правил c-стандарта. Например, GCC можно сказать, что вы вообще не заинтересованы в errno. Для этого передайте -fno-math-errno или -ffast-math . ICC и VC могут иметь похожие флаги компилятора.

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

Answer 2

Если вы собираетесь преобразовать результат операции floor() в целое число, и если вас не беспокоит переполнение, то следующий код намного быстрее, чем (int)floor(x):

inline int int_floor(double x)
{
  int i = (int)x; /* truncate */
  return i - ( i > x ); /* convert trunc to floor */
}

Answer 3

Пол и потолок без ответвлений (лучше использовать пипилин) без проверки ошибок

int f(double x)
{
    return (int) x - (x < (int) x); // as dgobbi above, needs less than for floor
}

int c(double x)
{
    return (int) x + (x > (int) x);
}

или используя пол

int c(double x)
{
    return -(f(-x));
}

Answer 4

Реальная самая быстрая реализация для большого массива на современных процессорах x86 будет

• изменить режим округления MXCSR FP на округление в направлении -Infinity (он же floor) . В C это должно быть возможно с fenv stuff или _mm_getcsr / _mm_setcsr.

• зацикливает массив, делая _mm_cvtps_epi32 для векторов SIMD, преобразуя 4 float с в 32-разрядное целое число, используя текущий режим округления. (И сохраняя векторы результатов к месту назначения.)

  cvtps2dq xmm0, [rdi] - это одиночный микроплавкий переход на любом процессоре Intel или AMD начиная с K10 или Core 2. (https://agner.org/optimize/) То же самое для 256-битной версии AVX с векторами YMM .

• восстановить текущий режим округления до нормального режима IEEE по умолчанию, используя исходное значение MXCSR. (округление до ближайшего, даже с разрывом связи)

Это позволяет загружать + преобразовывать + сохранять 1 вектор SIMD результатов за такт, так же быстро, как и с усечением . (SSE2 имеет специальную инструкцию преобразования FP-> int для усечения, именно потому, что это очень часто требуется компиляторам C. В старые добрые времена с x87 даже (int)x требовало изменить режим округления x87 на усечение, а затем обратно. cvttps2dq для упакованного числа с плавающей запятой-> int с усечением (обратите внимание на дополнительные t в мнемоническом выражении). Или для скалярного перехода от XMM к целочисленным регистрам, cvttss2si или cvttsd2si для скаляра double до скалярного целого числа.

При некотором развертывании цикла и / или хорошей оптимизации это должно быть возможно без узких мест на внешнем интерфейсе, только пропускная способность хранилища 1 раз в час, при условии отсутствия узких мест в кэше. (И на Intel до Skylake, также узкое место с пропускной способностью пакетного преобразования 1 на такт.), То есть 16, 32 или 64 байта за цикл, используя SSE2, AVX или AVX512.

---

Без изменения текущего режима округления необходимо SSE4.1 roundps, чтобы округлить float до ближайшего целого числа float, используя выбранные вами режимы округления. Или вы можете использовать один из трюков, показанных в других ответах, которые работают для чисел с плавающей запятой с достаточно малой величиной, чтобы поместиться в 32-разрядное целое число со знаком, так как это ваш конечный формат назначения в любом случае.

---

(С правильными опциями компилятора, такими как -fno-math-errno и правильными опциями -march или -msse4, компиляторы могут встраивать floor, используя roundps, или скаляр и / или эквивалент двойной точности, например roundsd xmm1, xmm0, 1, но это стоит 2 мопа и имеет пропускную способность 1 на 2 такта на Haswell для скаляров или векторов. На самом деле, gcc8.2 встроит roundsd для floor даже без каких-либо опций быстрой математики, , как вы можно увидеть в проводнике компилятора Godbolt . Но это с -march=haswell. К сожалению, это не базовый уровень для x86-64, поэтому его нужно включить, если ваш компьютер его поддерживает.)

Answer 5

На самом деле версия без ответвлений, которая требует одного преобразования между областями с плавающей запятой и целочисленными доменами, сместит значение x на весь положительный или весь отрицательный диапазон, затем приведёт / усечёт и сместит его обратно.

long fast_floor(double x)
{
    const unsigned long offset = ~(ULONG_MAX >> 1);
    return (long)((unsigned long)(x + offset) - offset);
}

long fast_ceil(double x) {
    const unsigned long offset = ~(ULONG_MAX >> 1);
    return (long)((unsigned long)(x - offset) + offset );
}

Как указано в комментариях, эта реализация опирается на временное значение x +- offset, которое не переполняется.

На 64-битных платформах исходный код, использующий промежуточное значение int64_t, приведет к созданию ядра из трех команд, то же самое доступно для int32_t с уменьшенным диапазоном floor / ceil, где |x| < 0x40000000 -

inline int floor_x64(double x) {
   return (int)((int64_t)(x + 0x80000000UL) - 0x80000000LL);
}
inline int floor_x86_reduced_range(double x) {
   return (int)(x + 0x40000000) - 0x40000000;
}

Answer 6

Да, floor() чрезвычайно медленный на всех платформах, поскольку он должен реализовывать много поведения из спецификации IEEE fp. Вы не можете использовать его во внутренних циклах.

Я иногда использую макрос для аппроксимации этажа ():

#define PSEUDO_FLOOR( V ) ((V) >= 0 ? (int)(V) : (int)((V) - 1))

Он не ведет себя точно так же, как floor(): например, floor(-1) == -1, но PSEUDO_FLOOR(-1) == -2, но достаточно близко для большинства применений.

Answer 7

Быстрый двойной раунд

double round(double x)
{
    return double((x>=0.5)?(int(x)+1):int(x));
}

---

Терминал журнала

test custom_1 8.3837

test native_1 18.4989

test custom_2 8.36333

test native_2 18.5001

test custom_3 8.37316

test native_3 18.5012

---

Тест

void test(char* name, double (*f)(double))
{
    int it = std::numeric_limits<int>::max();

    clock_t begin = clock();

    for(int i=0; i<it; i++)
    {
        f(double(i)/1000.0);
    }
    clock_t end = clock();

    cout << "test " << name << " " << double(end - begin) / CLOCKS_PER_SEC << endl;

}

int main(int argc, char **argv)
{

    test("custom_1",round);
    test("native_1",std::round);
    test("custom_2",round);
    test("native_2",std::round);
    test("custom_3",round);
    test("native_3",std::round);
    return 0;
}

---

Результат

Приведение типов и использование вашего мозга примерно в 3 раза быстрее, чем использование нативных функций.

Answer 8

1. Они не делают то же самое. Этаж () является функцией. Поэтому его использование вызывает вызов функции, выделение кадра стека, копирование параметров и получение результата. Приведение не является вызовом функции, поэтому оно использует более быстрые механизмы (я полагаю, что оно может использовать регистры для обработки значений).
2. Возможно, floor () уже оптимизирован.
3. Можете ли вы выжать из своего алгоритма больше производительности? Может быть, переключение строк и столбцов может помочь? Можете ли вы кешировать общие значения? Все ли оптимизации вашего компилятора включены? Можете ли вы переключить операционную систему? компилятор? Jon Bentley's Programming Pearls имеет большой обзор возможных оптимизаций.