42

有没有比使用 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与我想要钳制的值的类型相同(在上面的示例中为双精度)。

4

14 回答 14

53

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

三个指令(不包括 ret),没有分支。优秀的。

这是在带有 Core i3 M 350 的 Ubuntu 13.04 上使用 GCC 4.7 和 clang 3.2 测试的。顺便说一下,调用 std::min 和 std::max 的直接 C++ 代码生成了相同的程序集。

这是给双打的。对于 int,GCC 和 clang 都生成带有 5 条指令(不包括 ret)且没有分支的程序集。也很优秀。

我目前没有使用定点,所以我不会对定点发表意见。

于 2013-05-20T22:18:08.387 回答
40

老问题,但我今天正在解决这个问题(双打/浮点数)。

最好的方法是对浮点数使用 SSE MINSS/MAXSS,对双精度数使用 SSE2 MINSD/MAXSD。它们是无分支的,每个都需要一个时钟周期,并且由于编译器内在函数而易于使用。与使用 std::min/max 的钳位相比,它们赋予了超过一个数量级的性能提升。

你可能会感到惊讶。我当然做到了!不幸的是,即使启用了 /arch:SSE2 和 /FP:fast,VC++ 2010 也会对 std::min/max 使用简单的比较。我不能代表其他编译器。

这是在 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++ 编译器不够聪明,无法剔除多余的动作。少一个指令。:)

于 2011-07-03T22:52:14.173 回答
18

如果你的处理器有一个快速的绝对值指令(就像 x86 一样),你可以做一个无分支的 min 和 max,这将比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保存一个步骤。

在我的 Athlon II X2 上使用优化 FMIN=0 时,以下代码比三进制快 3 倍以上。

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
}
于 2011-09-15T01:26:50.710 回答
14

三元运算符确实是要走的路,因为大多数编译器都能够将它们编译成使用条件移动而不是分支的本机硬件操作(从而避免错误预测惩罚和管道泡沫等)。位操作可能会导致 load-hit-store

特别是,带有 SSE2 的 PPC 和 x86 有一个硬件操作,可以表示为这样的内在东西:

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 );
}

我强烈建议您避免使用整数运算对双精度进行位操作。在大多数现代 CPU 上,除了往返于 dcache 之外,没有直接的方法在 double 和 int 寄存器之间移动数据。这将导致称为 load-hit-store 的数据危害,它基本上清空 CPU 管道,直到内存写入完成(通常大约 40 个周期左右)。

例外情况是双精度值已经在内存中而不是在寄存器中:在这种情况下,不存在加载命中存储的危险。但是,您的示例表明您刚刚计算了双精度并从函数返回它,这意味着它可能仍在 XMM1 中。

于 2009-01-17T11:30:47.343 回答
9

对于 16.16 表示,简单的三进制不太可能在速度方面得到改进。

而对于双打,因为你需要它标准/便携式 C,任何类型的位摆弄都会以糟糕的方式结束。

即使有可能(我怀疑),你也会依赖双打的二进制表示。这(及其大小)取决于实施。

可能您可以使用 sizeof(double) 来“猜测”这一点,然后将各种 double 值的布局与其常见的二进制表示进行比较,但我认为您是在隐藏起来。

最好的规则是告诉编译器你想要什么(即三元),让它为你优化。

编辑:谦虚的馅饼时间。我刚刚测试了 quinmars 的想法(如下),它可以工作 - 如果你有 IEEE-754 浮点数。这给下面的代码带来了大约 20% 的加速。IObviously 不可移植,但我认为可能有一种标准化的方式来询问您的编译器是否使用带有#IF 的 IEEE754 浮点格式......?

  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);
于 2009-01-09T10:09:18.477 回答
7

IEEE 754 浮点位的排序方式是,如果您比较解释为整数的位,您将获得与直接将它们比较为浮点数相同的结果。因此,如果您找到或知道一种钳位整数的方法,您也可以将其用于 (IEEE 754) 浮点数。抱歉,我不知道更快的方法。

正如 rkj 所说,如果您将浮点数存储在数组中,则可以考虑使用一些 CPU 扩展,例如 SSE3。你可以看看 liboil 它为你做了所有的脏活。尽可能保持程序的可移植性并使用更快的 cpu 指令。(我不确定如何独立于操作系统/编译器的 liboil)。

于 2009-01-09T10:07:29.650 回答
7

我通常不使用测试和分支,而是使用这种格式进行钳位:

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

虽然我从来没有对编译后的代码做过任何性能分析。

于 2010-03-21T15:23:19.313 回答
4

实际上,没有像样的编译器会在 if() 语句和 ?: 表达式之间产生差异。代码很简单,他们将能够发现可能的路径。也就是说,您的两个示例并不相同。使用 ?: 的等效代码是

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

因为当 a > MAX 时避免 A < MIN 测试。现在这可能会有所作为,否则编译器将不得不发现两个测试之间的关系。

如果夹紧很少见,您可以通过一次测试来测试是否需要夹紧:

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

例如,当 MIN=6 和 MAX=10 时,这将首先将 a 向下移动 8,然后检查它是否位于 -2 和 +2 之间。这是否节省任何东西在很大程度上取决于分支的相对成本。

于 2009-01-09T10:06:13.507 回答
2

这是一个可能更快的实现,类似于@Roddy 的答案

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 */
}

请参阅计算两个整数的最小值 (min) 或最大值 (max) 而无需分支比较浮点数

IEEE 浮点和双精度格式的设计使数字是“按字典顺序排列的”,用 IEEE 架构师 William Kahan 的话来说,这意味着“如果对相同格式的两个浮点数进行排序(比如 x < y),那么当它们的位被重新解释为符号幅度整数时,它们的排序方式相同。”</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)
于 2009-01-09T19:36:09.200 回答
1

我自己尝试了 SSE 方法,并且汇编输出看起来更干净,所以一开始我很受鼓舞,但在计时了数千次之后,它实际上慢了很多。看起来 VC++ 编译器确实不够聪明,无法知道你真正想要什么,而且它似乎在不应该的时候在 XMM 寄存器和内存之间来回移动。也就是说,当编译器似乎对所有浮点计算都使用 SSE 指令时,我不知道为什么编译器不够聪明,无法在三元运算符上使用 SSE min/max 指令。另一方面,如果您正在为 PowerPC 进行编译,则可以在 FP 寄存器上使用 fsel 内在函数,而且速度更快。

于 2011-08-19T17:39:20.610 回答
1

如上所述,fmin/fmax 函数运行良好(在 gcc 中,使用 -ffast-math)。尽管 gfortran 具有使用对应于 max/min 的 IA 指令的模式,但 g++ 没有。在 icc 中必须使用 std::min/max,因为 icc 不允许缩短 fmin/fmax 如何与非有限操作数一起工作的规范。

于 2015-11-14T17:22:09.137 回答
1

我在 C++ 中的 2 美分。可能与使用三元运算符没有什么不同,希望不会生成分支代码

template <typename T>
inline T clamp(T val, T lo, T hi) {
    return std::max(lo, std::min(hi, val));
}
于 2017-10-26T20:32:43.807 回答
0

如果我理解正确,您希望将值“a”限制在 MY_MIN 和 MY_MAX 之间的范围内。“a”的类型是双精度。您没有指定 MY_MIN 或 MY_MAX 的类型。

简单的表达:

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

应该做的伎俩。

如果 MY_MAX 和 MY_MIN 恰好是整数,我认为可能需要进行一些小的优化:

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

通过更改为整数比较,您可能会获得轻微的速度优势。

于 2009-01-09T09:21:45.517 回答
0

如果您想使用快速绝对值指令,请查看我在 minicomputer 中找到的这段代码片段,它将浮点数限制在 [0,1] 范围内

clamped = 0.5*(fabs(x)-fabs(x-1.0f) + 1.0f);

(我稍微简化了代码)。我们可以认为它取两个值,一个反映为 >0

fabs(x)

另一个反映大约 1.0 为 <1.0

1.0-fabs(x-1.0)

我们取它们的平均值。如果它在范围内,那么这两个值将与 x 相同,因此它们的平均值将再次为 x。如果超出范围,则其中一个值为 x,另一个值为 x 翻转“边界”点,因此它们的平均值将恰好是边界点。

于 2012-05-20T07:46:22.450 回答