Самый быстрый способ зафиксировать реальное (фиксированное / с плавающей запятой) значение? - PullRequest
37 голосов
/ 09 января 2009

Есть ли более эффективный способ фиксирования действительных чисел, чем использование операторов if или троичных операторов? Я хочу сделать это как для двойников, так и для 32-битной реализации с фиксированной точкой (16.16). Я не прошу код, который может обрабатывать оба случая; они будут обрабатываться в отдельных функциях.

Очевидно, я могу сделать что-то вроде:

double clampedA;
double a = calculate();
clampedA = a > MY_MAX ? MY_MAX : a;
clampedA = a < MY_MIN ? MY_MIN : a;

или

double a = calculate();
double clampedA = a;
if(clampedA > MY_MAX)
    clampedA = MY_MAX;
else if(clampedA < MY_MIN)
    clampedA = MY_MIN;

Версия с фиксированной точкой будет использовать функции / макросы для сравнения.

Это сделано в критически важной для кода части, поэтому я ищу как можно более эффективный способ сделать это (что, я подозреваю, потребует битовых манипуляций)

РЕДАКТИРОВАТЬ: Это должен быть стандартный / переносимый C, функциональные возможности платформы здесь не представляют интереса. Кроме того, MY_MIN и MY_MAX имеют тот же тип, что и значение, которое я хочу зафиксировать (в приведенных выше примерах удваивается).

Ответы [ 14 ]

37 голосов
/ 04 июля 2011

Старый вопрос, но я работал над этой проблемой сегодня (с удвоениями / числами с плавающей запятой).

Наилучший подход - использовать SSE MINSS / MAXSS для поплавков и SSE2 MINSD / MAXSD для двойных. Они не требуют ответвлений и занимают один тактовый цикл, и их легко использовать благодаря встроенным функциям компилятора. Они обеспечивают увеличение производительности более чем на порядок по сравнению с зажимом с помощью std :: min / max.

Вы можете найти это удивительным. Я конечно сделал! К сожалению, VC ++ 2010 использует простые сравнения для std :: min / max, даже когда включены / arch: SSE2 и / FP: fast. Я не могу говорить за другие компиляторы.

Вот код, необходимый для этого в VC ++:

#include <mmintrin.h>

float minss ( float a, float b )
{
    // Branchless SSE min.
    _mm_store_ss( &a, _mm_min_ss(_mm_set_ss(a),_mm_set_ss(b)) );
    return a;
}

float maxss ( float a, float b )
{
    // Branchless SSE max.
    _mm_store_ss( &a, _mm_max_ss(_mm_set_ss(a),_mm_set_ss(b)) );
    return a;
}

float clamp ( float val, float minval, float maxval )
{
    // Branchless SSE clamp.
    // return minss( maxss(val,minval), maxval );

    _mm_store_ss( &val, _mm_min_ss( _mm_max_ss(_mm_set_ss(val),_mm_set_ss(minval)), _mm_set_ss(maxval) ) );
    return val;
}

Код двойной точности такой же, за исключением xxx_sd. ​​

Редактировать: Первоначально я написал функцию зажима в качестве комментария. Но, глядя на вывод ассемблера, я заметил, что компилятор VC ++ не был достаточно умен, чтобы отбросить избыточный ход. На одну инструкцию меньше. :)

36 голосов
/ 21 мая 2013

И GCC, и Clang генерируют красивую сборку для следующего простого, простого, переносимого кода:

double clamp(double d, double min, double max) {
  const double t = d < min ? min : d;
  return t > max ? max : t;
}

> gcc -O3 -march=native -Wall -Wextra -Wc++-compat -S -fverbose-asm clamp_ternary_operator.c

GCC-генерируемая сборка:

maxsd   %xmm0, %xmm1    # d, min
movapd  %xmm2, %xmm0    # max, max
minsd   %xmm1, %xmm0    # min, max
ret

> clang -O3 -march=native -Wall -Wextra -Wc++-compat -S -fverbose-asm clamp_ternary_operator.c

Сгенерированная Clang сборка:

maxsd   %xmm0, %xmm1
minsd   %xmm1, %xmm2
movaps  %xmm2, %xmm0
ret

Три инструкции (не считая рет), веток нет. Отлично.

Это было протестировано с GCC 4.7 и clang 3.2 на Ubuntu 13.04 с Core i3 M 350. Кстати, простой код C ++, вызывающий std :: min и std :: max, сгенерировал одну и ту же сборку.

Это для двойников. А для int и GCC, и clang генерируют сборку с пятью инструкциями (не считая ret) и без ветвей. Также отлично.

В настоящее время я не использую фиксированную точку, поэтому я не буду давать заключение о фиксированной точке.

15 голосов
/ 15 сентября 2011

Если ваш процессор имеет быструю инструкцию для абсолютного значения (как у x86), вы можете сделать минимальные и максимальные значения без ветвей, которые будут быстрее, чем оператор if или троичная операция.

min(a,b) = (a + b - abs(a-b)) / 2
max(a,b) = (a + b + abs(a-b)) / 2

Если один из терминов равен нулю (как это часто бывает при фиксации), код немного упрощается:

max(a,0) = (a + abs(a)) / 2

Когда вы объединяете обе операции, вы можете заменить две /2 на одну /4 или *0.25, чтобы сохранить шаг.

Следующий код более чем в 3 раза быстрее троичного на моем Athlon II X2 при использовании оптимизации для FMIN = 0.

double clamp(double value)
{
    double temp = value + FMAX - abs(value-FMAX);
#if FMIN == 0
    return (temp + abs(temp)) * 0.25;
#else
    return (temp + (2.0*FMIN) + abs(temp-(2.0*FMIN))) * 0.25;
#endif
}
14 голосов
/ 17 января 2009

Тернарный оператор - действительно путь, потому что большинство компиляторов могут скомпилировать их в собственную аппаратную операцию, которая использует условное перемещение вместо ветви (и, таким образом, избегает штрафов за неправильное предсказание и пузырьков конвейера и т. Д.). Манипулирование битами может вызвать загрузку хранилища .

В частности, PPC и x86 с SSE2 имеют аппаратную операционную систему, которая может быть выражена как нечто вроде этого:

double fsel( double a, double b, double c ) {
  return a >= 0 ? b : c; 
}

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

inline double clamp ( double a, double min, double max ) 
{
   a = fsel( a - min , a, min );
   return fsel( a - max, max, a );
}

Я настоятельно рекомендую вам избегать битовых манипуляций с двойными числами с помощью целочисленных операций . На большинстве современных процессоров нет прямого способа перемещения данных между двойными и целыми регистрами, кроме как в обратном направлении к dcache. Это приведет к опасности данных, называемой хранилищем попадания нагрузки, которая в основном очищает конвейер ЦП до тех пор, пока не завершится запись в память (обычно около 40 циклов или около того).

Исключением из этого является случай, когда двойные значения уже находятся в памяти, а не в регистре: в этом случае нет опасности сохранения загрузки. Однако ваш пример показывает, что вы только что вычислили double и вернули его из функции, что означает, что он, вероятно, все еще будет в XMM1.

8 голосов
/ 09 января 2009

Для представления 16.16 простая троичная фигура вряд ли будет улучшена по скорости.

А для парных разрядов, потому что вам нужен стандартный / переносной язык Си, любая битовая битва плохо закончится.

Даже если бы была возможна бит-скрипка (в чем я сомневаюсь), вы бы полагались на двоичное представление парных чисел. ЭТО (и их размер) зависит от реализации.

Возможно, вы могли бы "угадать" это, используя sizeof (double), а затем сравнивать расположение различных значений double с их общими двоичными представлениями, но я думаю, что вы скрываетесь от нуля.

Лучшее правило - СКАЖИТЕ КОМПИЛЕРУ, ЧТО ВЫ ХОТИТЕ (т. Е. Троичному), и позвольте ему оптимизировать для вас.

РЕДАКТИРОВАТЬ: Скромное время пирога. Я только что проверил идею Quinmars (ниже), и она работает - если у вас есть поплавки IEEE-754. Это дало ускорение примерно на 20% по коду ниже. Очевидно, что он не переносимый, но я думаю, что может существовать стандартизированный способ запроса вашего компилятора, использует ли он форматы float IEEE754 с #IF ...?

  double FMIN = 3.13;
  double FMAX = 300.44;

  double FVAL[10] = {-100, 0.23, 1.24, 3.00, 3.5, 30.5, 50 ,100.22 ,200.22, 30000};
  uint64  Lfmin = *(uint64 *)&FMIN;
  uint64  Lfmax = *(uint64 *)&FMAX;

    DWORD start = GetTickCount();

    for (int j=0; j<10000000; ++j)
    {
        uint64 * pfvalue = (uint64 *)&FVAL[0];
        for (int i=0; i<10; ++i)
            *pfvalue++ = (*pfvalue < Lfmin) ? Lfmin : (*pfvalue > Lfmax) ? Lfmax : *pfvalue;
    }

    volatile DWORD hacktime = GetTickCount() - start;

    for (int j=0; j<10000000; ++j)
    {
        double * pfvalue = &FVAL[0];
        for (int i=0; i<10; ++i)
            *pfvalue++ = (*pfvalue < FMIN) ? FMIN : (*pfvalue > FMAX) ? FMAX : *pfvalue;
    }

    volatile DWORD normaltime = GetTickCount() - (start + hacktime);
7 голосов
/ 21 марта 2010

Вместо тестирования и ветвления я обычно использую этот формат для зажима:

clampedA = fmin(fmax(a,MY_MIN),MY_MAX);

Хотя я никогда не проводил анализ производительности скомпилированного кода.

7 голосов
/ 09 января 2009

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

Если у вас есть числа с плавающей запятой, хранящиеся в массивах, вы можете использовать некоторые расширения ЦП, такие как SSE3, как сказал rkj. Вы можете взглянуть на liboil, он делает всю грязную работу за вас. Держит вашу программу переносимой и использует более быстрые инструкции процессора, если это возможно. (Я не уверен, насколько независим от OS / компилятора liboil).

4 голосов
/ 09 января 2009

Реально, ни один достойный компилятор не будет иметь разницы между оператором if () и выражением? :. Код достаточно прост, чтобы они могли определить возможные пути. Тем не менее, ваши два примера не идентичны. Эквивалентный код, использующий?: Будет

a = (a > MAX) ? MAX : ((a < MIN) ? MIN : a);

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

Если зажим является редким, вы можете проверить необходимость зажима с помощью одного теста:

if (abs(a - (MAX+MIN)/2) > ((MAX-MIN)/2)) ...

например. при MIN = 6 и MAX = 10 это сначала сместит вниз на 8, а затем проверит, находится ли оно между -2 и +2. Сохраняет ли это что-нибудь, во многом зависит от относительной стоимости ветвления.

2 голосов
/ 09 января 2009

Вот, возможно, более быстрая реализация, похожая на @ ответ Родди :

typedef int64_t i_t;
typedef double  f_t;

static inline
i_t i_tmin(i_t x, i_t y) {
  return (y + ((x - y) & -(x < y))); // min(x, y)
}

static inline
i_t i_tmax(i_t x, i_t y) {
  return (x - ((x - y) & -(x < y))); // max(x, y)
}

f_t clip_f_t(f_t f, f_t fmin, f_t fmax)
{
#ifndef TERNARY
  assert(sizeof(i_t) == sizeof(f_t));
  //assert(not (fmin < 0 and (f < 0 or is_negative_zero(f))));
  //XXX assume IEEE-754 compliant system (lexicographically ordered floats)
  //XXX break strict-aliasing rules
  const i_t imin = *(i_t*)&fmin;
  const i_t imax = *(i_t*)&fmax;
  const i_t i    = *(i_t*)&f;
  const i_t iclipped = i_tmin(imax, i_tmax(i, imin));

#ifndef INT_TERNARY
  return *(f_t *)&iclipped;
#else /* INT_TERNARY */
  return i < imin ? fmin : (i > imax ? fmax : f); 
#endif /* INT_TERNARY */

#else /* TERNARY */
  return fmin > f ? fmin : (fmax < f ? fmax : f);
#endif /* TERNARY */
}

См. Вычисление минимального (минимального) или максимального (максимального) двух целых чисел без ветвления и Сравнение чисел с плавающей запятой

IEEE float и double форматы были разработан так, чтобы числа «Лексикографически упорядоченный», который - по словам IEEE архитектора Уильяма Кахан означает «если два с плавающей точкой номера в том же формате упорядочены (скажем, х <у), то они упорядочены так же, когда их биты переосмысливается как Величина Знака целые числа.»</p>

Тестовая программа:

/** gcc -std=c99 -fno-strict-aliasing -O2 -lm -Wall *.c -o clip_double && clip_double */
#include <assert.h> 
#include <iso646.h>  // not, and
#include <math.h>    // isnan()
#include <stdbool.h> // bool
#include <stdint.h>  // int64_t
#include <stdio.h>

static 
bool is_negative_zero(f_t x) 
{
  return x == 0 and 1/x < 0;
}

static inline 
f_t range(f_t low, f_t f, f_t hi) 
{
  return fmax(low, fmin(f, hi));
}

static const f_t END = 0./0.;

#define TOSTR(f, fmin, fmax, ff) ((f) == (fmin) ? "min" :       \
                  ((f) == (fmax) ? "max" :      \
                   (is_negative_zero(ff) ? "-0.":   \
                    ((f) == (ff) ? "f" : #f))))

static int test(f_t p[], f_t fmin, f_t fmax, f_t (*fun)(f_t, f_t, f_t)) 
{
  assert(isnan(END));
  int failed_count = 0;
  for ( ; ; ++p) {
    const f_t clipped  = fun(*p, fmin, fmax), expected = range(fmin, *p, fmax);
    if(clipped != expected and not (isnan(clipped) and isnan(expected))) {
      failed_count++;
      fprintf(stderr, "error: got: %s, expected: %s\t(min=%g, max=%g, f=%g)\n", 
          TOSTR(clipped,  fmin, fmax, *p), 
          TOSTR(expected, fmin, fmax, *p), fmin, fmax, *p);
    }
    if (isnan(*p))
      break;
  }
  return failed_count;
}  

int main(void)
{
  int failed_count = 0;
  f_t arr[] = { -0., -1./0., 0., 1./0., 1., -1., 2, 
        2.1, -2.1, -0.1, END};
  f_t minmax[][2] = { -1, 1,  // min, max
               0, 2, };

  for (int i = 0; i < (sizeof(minmax) / sizeof(*minmax)); ++i) 
    failed_count += test(arr, minmax[i][0], minmax[i][1], clip_f_t);      

  return failed_count & 0xFF;
}

В консоли:

$ gcc -std=c99 -fno-strict-aliasing -O2 -lm *.c -o clip_double && ./clip_double 

Он печатает:

error: got: min, expected: -0.  (min=-1, max=1, f=0)
error: got: f, expected: min    (min=-1, max=1, f=-1.#INF)
error: got: f, expected: min    (min=-1, max=1, f=-2.1)
error: got: min, expected: f    (min=-1, max=1, f=-0.1)
1 голос
/ 19 августа 2011

Я сам попробовал подход SSE к этому, и вывод сборки выглядел немного чище, поэтому сначала я был воодушевлен, но после синхронизации с тысячами раз он был на самом деле немного медленнее. Он действительно выглядит так, как будто компилятор VC ++ не достаточно умен, чтобы знать, что вы на самом деле намереваетесь, и, кажется, он перемещает вещи назад и назад между регистрами XMM и памятью, когда этого не следует делать. Тем не менее, я не знаю, почему компилятор не настолько умен, чтобы использовать инструкции SSE min / max для троичного оператора, когда он все равно использует инструкции SSE для всех вычислений с плавающей запятой. С другой стороны, если вы компилируете для PowerPC, вы можете использовать встроенную функцию fsel в регистрах FP, и это намного быстрее.

...