如何最快地夹紧实数（定点/浮点）值？

Question

有没有比使用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与我想要夹紧的值相同类型（例如上面的双倍）。

- Niklas

---

我认为你可以使用SSE3或类似的技术来实现这个，但不确定具体需要哪些指令或如何实现... 你可以参考一下：饱和算术 - rkj

抱歉，这个问题没有明确平台要求。我已经编辑了问题以使其更加清晰。 - Niklas

我知道你提出这个问题已经两年半了，但我希望你看看我的回答——三倍的改进是非常显著的。 - Mark Ransom

一个未指定的细节是你愿意为速度而交换哪种精度（相对或绝对）-如果有的话。如果代码要求在范围内返回一个确切的 a，那么许多答案都无法满足这个要求。如果精度不是问题，那么始终返回 (MY_MAX + MY_MIN)/2 将肯定是一个快速低精度的答案，也肯定是愚蠢的。建议容忍不超过1个 ULP 的误差。 - chux - Reinstate Monica

你会如何在SSE4变量（__m128）上执行它？ - Royi

有关 最高效优雅的截断数字的方法 - Trevor Boyd Smith

Answer 1

无论是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），没有分支。非常棒。

这是在Ubuntu 13.04上使用GCC 4.7和clang 3.2对Core i3 M 350进行测试的。顺便说一句，直接调用std::min和std::max的C++代码生成的汇编代码也是相同的。

这是针对double的。对于int，GCC和clang都会生成包含五个指令（不包括ret）且没有分支的汇编代码。同样很好。

我目前不使用定点数，因此我不会对定点数发表意见。

- Jorge

---

太好了。这个答案略优于好的回答，因为它对min和/或max中有一个或两个是非数字的情况进行了对称处理。此外，它还保留了d=-0.0的符号！ - chux - Reinstate Monica

使用 if (d < min) 和 if (d > max) 也会给我相同的汇编代码。然而，有趣的是，使用 if (d < min) 和 else if (d > max) 会生成不同的输出（有一个跳转指令）。 - elboulangero

这是关于问题的编译器分析：https://godbolt.org/z/ZW4W6F - Felipe Lavratti

已在MSVC 2019下测试，也可以编译为无分支代码（至少对于浮点数是这样）。 - Gabriele Giuseppini

Answer 2

虽然这是一个旧问题，但我今天正在解决它（使用双精度/浮点数）。

最佳方法是使用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++ 编译器不够聪明，无法消除多余的移动指令。 少了一条指令。 :)

- Spat

---

GCC有没有这些函数的等价物？ - Dan

是的，对于GCC x86，请使用__builtin_ia32_storess、__builtin_ia32_maxss和__builtin_ia32_minss等等相应的函数，并使用xmmintrin.h头文件来进行SSE1指令。在编译器中传递-mmmx -msse参数，您可能还需要-mfpmath=sse(,x87)参数。ARM Neon和AltiVec也提供了内部函数。有关更多详细信息，请参见X86内置函数。 - mctylr

编译器通常无法在一般情况下用内置函数替换std::min和std::max，因为内置函数提供了IEEE754规定的结果，例如min(2.0, NaN)和min(NaN, 2.0)（两种情况下都是2.0），而基于单个比较的天真实现将根据参数顺序返回不一致的结果。 C99和C++11提供了fmax和fmin，聪明的编译器会用高效的内联实现替换它们。 - strcat

使用SSE指令或将它们与标准浮点运算交错使用是否会有切换开销？ - Robinson

这似乎非常有帮助 --- 有人知道任何完整的实现吗？例如，对于gcc和clang等，请使用适当的#ifdef吗？ - Bob

Answer 3

如果您的处理器有一个快速的绝对值指令（如x86），则可以使用无分支的最小值和最大值，这将比if语句或三元操作更快。

如果您的处理器具有快速的绝对值指令（例如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 以节省一步。

在我的 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
}

- Mark Ransom

---

哇，好主意！我怀疑在某些CPU/编译器上，如果abs(a)没有被很好地内联/优化，这可能比三元运算符更慢... - Roddy

在C#中，使用Math.Abs这种方法会比较慢。 - Paul Chernoch

我会期望使用 fabs(value-FMAX) 而不是 int abs(int j)。 - chux - Reinstate Monica

@chux 我使用了一个 C++ 编译器进行测试，可能通过重载使用了适当的函数。 - Mark Ransom

弱点：这种方法可能会导致严重的精度损失。大于“value”的FMAX值在结果中可能会失去精度。如果FMAX是value的10倍，则有可能会失去1位小数。最糟糕的情况是，返回值始终为0.0。 - chux - Reinstate Monica

你需要说对于无符号值 min(a,b) = (a + b - abs(static_cast<int>(a-b))) / 2 吗？如果没有 static_cast<int>，那么如果 b > a，a-b 可能会是一个非常大的值，从而导致计算出错误的值。 - dgnuff

Answer 4

三目运算符确实是最好的选择，因为大多数编译器能够将它们编译为本机硬件操作，使用条件移动而不是分支（从而避免了错误预测惩罚、流水线泡沫等）。位操作可能会造成加载-存储。

特别是PPC和带有SSE2的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（数据高速缓存）之外，没有直接的方式在双精度浮点寄存器和整数寄存器之间移动数据。这将导致一种称为load-hit-store的数据危险，基本上会清空CPU流水线，直到内存写入完成（通常需要约40个周期左右）。

唯一的例外是如果双精度浮点值已经在内存中而不在寄存器中：这种情况下不存在load-hit-store的危险。然而，您的示例表明您刚刚计算出双精度浮点数并从函数返回它，这意味着它很可能仍然在XMM1中。

- Crashworks

---

有关三元运算符的一个注意事项：测试输入的类型和顺序如何影响优化后的输出。我曾在一个编译器上工作过，其中 A > B ? A : B 一直生成一个 MAX 指令，但 A < B ? B : A 却没有。 - AShelly

@AShelly ：你不得不想知道那个编译器的作者当时在想什么。 - Crashworks

所有的FP数字都可以很好地工作！它甚至保留了a == -0.0的符号！我只对一些不对称的值/限制有些担忧，这些值/限制涉及到Not-a-numbers。允许min是一个Not-a-number，并且很好地忽略了min。但是，如果max是NAN，则结果也是NAN。可以通过与return fsel( a - max, max, a );不同的代码使其对称。 - chux - Reinstate Monica

Answer 5

对于16.16表示法，简单的三元运算速度上不太可能有更好的表现。

对于双精度浮点数，因为你需要标准/可移植的C语言，任何位操作都会以失败告终。

即使位操作是可能的（我对此表示怀疑），你也要依赖于双精度浮点数的二进制表示。这个（以及它们的大小）是与实现相关的。

可能你可以使用sizeof(double)来“猜测”，然后将各种双精度浮点值的布局与它们的常见二进制表示进行比较，但我认为你会一无所获。

最好的方法是告诉编译器你想要什么（即三元运算符），让它为你进行优化。

编辑：请容许我承认错误。我刚测试了quinmars的想法（如下），它行得通——如果你有IEEE-754浮点数。这使下面的代码加速了约20%。显然不可移植，但我认为可能有一种标准化的方法来询问你的编译器是否使用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);

- Roddy

---

假设在我的情况下，使用IEEE-754浮点数是足够可移植的。感谢您抽出时间跟进。 - Niklas

当FMIN和*pfvalue都小于零时，带有int64_t版本将给出错误的结果，例如，FMIN=-1，FMAX=1，(*pfvalue)=-0.1; 请参见我的答案https://dev59.com/lnRC5IYBdhLWcg3wCMrX。 - jfs

@JFS 是的，IEE754使用符号/幅度编码，而不是2s补码。因此，与负数的比较是有缺陷的。如果FMIN和FMAX都>=0，则没问题（即使pfvalue为负）。如果FMIN或FMAX为零，则一切皆有可能... - Roddy

我想知道你是否有时间比较一下我的无分支min/max解决方案和你的解决方案？我很希望得到一些独立的验证，特别是因为我无法使用quinmars版本复制你的结果。 - Mark Ransom

@Mark - 我会尽力而为。不同的结果可能是因为你的编译器比我的优化得更好！ - Roddy

Answer 6

与其进行测试和分支，我通常使用以下格式进行夹紧：

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

虽然我从未对编译后的代码进行过任何性能分析。

- Linasses

---

不错。任何替代代码都应该以此为标准进行性能测试，同时在功能上匹配。 - chux - Reinstate Monica

+1 不清楚速度如何，但是当性能不是一个问题时，这绝对比涉及三元运算符等解决方案更简洁。 - charlescochran

Answer 7

IEEE 754浮点数的位以一种方式排序，使得如果将其解释为整数进行比较，则可获得与直接将其解释为浮点数进行比较相同的结果。因此，如果您找到或知道限制整数的方法，也可以将其用于（IEEE 754）浮点数。很抱歉，我不知道更快的方法。

如果您在数组中存储了浮点数，则可以考虑使用一些CPU扩展，例如SSE3，正如rkj所说。您可以查看liboil，它会为您完成所有繁琐的工作。保持程序的可移植性，并在可能时使用更快的CPU指令。（我不确定liboil在操作系统/编译器方面是否独立。）

- quinmars

---

仅适用于正浮点数。如果符号可能混合，您需要注意它们，如果不同则提前返回，如果为负数，则取绝对值并反转顺序。简而言之，优化仅适用于正浮点数。 - Potatoswatter

Answer 8

实际上，任何一个好的编译器都不会区分 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 之间。这是否会节省任何东西取决于分支的相对成本。

- MSalters

---

你会惊讶——上次我查看它的反汇编时，我的编译器将三元表达式正确地转换为相应的条件移动操作码，但它却将等效的if/else代码块转换为两个比较和分支。 - Crashworks

我喜欢使用单个测试夹紧的想法 ;) - e.tadeu

我正在寻找一种快速的方法来测试一个点是否在边界框内。这意味着测试X值是否在最大值和最小值之间，Y值也是如此。您的建议看起来很有前途。 - Paul Chernoch

1. 使用 fabs() 要比 int abs(int) 更好。
2. 在 fabs(a - (MAX+MIN)/2) > ((MAX-MIN)/2) 中存在精度损失的边缘情况问题。第一种方法没有这些问题。

- chux - Reinstate Monica

@MSalters，我在想这个#define CLAMP(VAL, LO, HI) VAL = ((VAL < LO) ? LO : ((VAL > HI) ? HI : VAL)，你认为编译器会充分优化吗？我正在使用64位ARM GCC。 - Ganindu

@Ganindu：编译器在过去的十年中肯定有所改进。我不会太担心。最坏的情况是它会出现在性能分析中。 - MSalters

Answer 9

这里有一个可能更快的实现方法，类似于@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架构师威廉卡汉的话来说：“如果在相同格式中排序了两个浮点数（例如x<y），则当它们的位被重新解释为符号-大小整数时，它们以相同的方式排序。”

一个测试程序：

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

- J.F. Sebastian

---

关于 is_negative_zero，你为什么没有使用 C99 标准库中的 math.h 中的 signbit 函数（与 x == 0 结合使用），而是使用了 1.0 / x < 0 来检查零的符号？ - mctylr

@mctylr：不记得了。 signbit似乎也可以正常工作（http://ideone.com/zqMl0）。 - jfs

Answer 10

我在C++中的两分钱。可能和使用三元运算符没有什么区别，希望不会生成任何分支代码。

template <typename T>
inline T clamp(T val, T lo, T hi) {
    return std::max(lo, std::min(hi, val));
}

- wcochran

---

只是给偶然看到这篇文章的人提供一些信息：C++17 的 <algorithm> 头文件引入了 std::clamp(n, low, high [, compare])。 - Tony Delroy