钳制真实（固定/浮点）值的最快方法？

Posted 2023-05-08

【中文标题】钳制真实（固定/浮点）值的最快方法？

【英文标题】：Fastest way to clamp a real (fixed/floating point) value?

【发布时间】：2010-09-30 11:02:21

【问题描述】：

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

【问题讨论】：

我觉得你可以使用SSE3或者一些类似的技术来做这个，但是不知道具体是哪些命令/如何...你可以看看：Saturation arithmetic

抱歉，关于平台要求的问题不清楚。我已将问题编辑为有点清楚。

我知道你问这个问题已经两年半了，但我希望你检查我的回答 - 3 倍的改进是显着的。

一个未指定的细节是您愿意以何种精度（相对或绝对）来换取速度 - 如果有的话。如果代码要求返回范围内的a 与a 完全一样，那么许多答案都不会遇到这个障碍。如果精度与 no 无关，那么总是返回 (MY_MAX + MY_MIN)/2 肯定是一个快速的低精度答案，而且肯定是愚蠢的。建议容忍不超过 1 个ULP 错误。

你将如何处理 SSE4 变量 (__m128)？

【参考方案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），没有分支。太棒了。

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

这是双打。对于 int，GCC 和 clang 都生成带有五个指令（不包括 ret）且没有分支的程序集。也很棒。

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

【讨论】：

太棒了。比good answer 稍微好一点，因为它对称地处理min 和/或max 当一个/两个都是非数字时。它还保留带有d = -0.0的符号！

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

准确。这应该是正确的答案。这是对问题的编译器分析：godbolt.org/z/ZW4W6F

使用 MSVC 2019 测试，也编译为无分支（至少对于浮点数）。

【参考方案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++ 编译器不够聪明，无法剔除多余的动作。少一个指令。 :)

【讨论】：

GCC 的这些功能是否有等价物？

是的，对于 GCC x86，使用 __builtin_ia32_storess、__builtin_ia32_maxss、__builtin_ia32_minss` 是等效的函数和 SSE1 指令的 xmmintrin.h 标头。将-mmmx -msse 传递给编译器，您可能还需要-mfpmath=sse(,x87)。内在函数也可用于 ARM Neon 和 AltiVec。有关详细信息，请参阅X86 Built-in functions。

在一般情况下，编译器不可能用内部函数替换 std::min 和 std::max，因为内部函数为 min(2.0, NaN) 和 min(NaN, 2.0) 提供了 IEEE754 指定的结果（即2.0 在这两种情况下），而基于单个比较的幼稚实现将根据参数顺序返回不一致的结果。 C99 和 C++11 提供了 fmax 和 fmin，一个聪明的编译器将用高效的内联实现替换它们。

使用 SSE 指令或将它们与标准浮点操作交错是否有切换惩罚？

这似乎真的很有帮助 --- 有没有人知道任何地方的完整实现，例如，为 gcc 和 clang 等使用适当的#ifdef？

【参考方案3】：

如果您的处理器具有绝对值的快速指令（如 x86 那样），您可以执行无分支的 min 和 max，这将比 if 语句或三元运算更快。

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

如果其中一项为 0（当您进行钳位时通常是这种情况），则代码会进一步简化：

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

【讨论】：

哇——好主意！我怀疑在某些 CPU/编译器上，如果 abs(a) 没有很好地内联/优化，这实际上可能比三进制慢...

在 C# 中，使用 Math.Abs​​，这种方法比较慢。

我希望fabs(value-FMAX) 而不是int abs(int j)。

@chux 我用 C++ 编译器测试过，它可能通过重载使用了正确的函数。

弱点：这种方法会导致精度损失。大于value 的FMAX 值可能会丢失结果的精度。如果 FMAX 是 value 的 10 倍，则可能会丢失 1 个十进制数字。最坏的情况，钳制的返回值总是 0.0。

【参考方案4】：

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

特别是，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 之外，没有直接的方法可以在 double 和 int 寄存器之间移动数据。这将导致称为 load-hit-store 的数据危害，它基本上清空 CPU 管道，直到内存写入完成（通常大约 40 个周期左右）。

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

【讨论】：

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

@AShelly ：您一定想知道编写该编译器的人在想什么。

适用于所有 FP 号码！它甚至保留了a == -0.0 的标志！只有我担心的值/限制涉及非数字的一些不对称性 - 允许 min 成为非数字并很好地忽略 min。然而，如果max 是 NAN，则结果是 NAN。可以使用不同于return fsel( a - max, max, a );

的代码实现对称

【参考方案5】：

对于 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);

【讨论】：

假设浮点情况下的 IEEE-754 对于我的情况来说足够便携。感谢您抽出宝贵时间跟进。

当FMIN 和*pfvalue 都小于零时，带有int64_t 的版本会给出错误的结果，例如，FMIN=-1, FMAX=1, (*pfvalue)=-0.1;看我的回答***.com/questions/427477/…

@JFS 啊，是的。 IEE754 使用符号/幅度编码，而不是 2s 补码。所以与负数的比较是有缺陷的。如果 FMIN 和 FMAX 都 >= 零，那么你很好（即使 pfvalue 是负数）。如果 FMIN 或 FMAX 为零，则所有赌注都关闭...

我想知道您是否有时间将我的无分支最小/最大解决方案与您的进行比较？我想要一些独立的验证，特别是因为我无法使用 quinmars 版本复制您的结果。

@Mark - 我会看看我能做什么。不同的结果可能是因为你的编译器刚刚优化了一个负载比我的好！

【参考方案6】：

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

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

【讨论】：

仅适用于正浮点数。如果符号可能混合，则需要注意它们，如果不同则尽早返回，如果为负则取绝对值并反向排序。简而言之，优化只适用于正浮点数。

【参考方案7】：

我通常使用这种格式进行钳制，而不是测试和分支：

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

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

【讨论】：

不错。任何替代代码都应以此为标准进行测试，以在性能上超越标准，但在功能上相匹配。

【参考方案8】：

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

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

因为当 a > MAX 时避免 A

如果钳位很少，可以单次测试是否需要钳位：

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

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

【讨论】：

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

我喜欢用一个测试来夹紧的想法；）

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

1) 期望fabs() 比int abs(int) 2) 边缘条件问题与fabs(a - (MAX+MIN)/2) > ((MAX-MIN)/2) 中的精度损失有关。第一种方法没有这些问题。

@Ganindu：编译器在过去十年中确实得到了改进。我不会担心的。更糟糕的情况是它会出现在分析中。

【参考方案9】：

这是一个可能更快的实现，类似于@Roddy's answer：

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

见Compute the minimum (min) or maximum (max) of two integers without branching和Comparing floating point numbers

IEEE 浮点和双精度格式是 设计使数字是 “按字典顺序排列”，其中—— 用 IEEE 架构师威廉的话来说 Kahan 的意思是“如果两个浮点 相同格式的数字是有序的 （比如 x

一个测试程序：

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

【讨论】：

Re: is_negative_zero，请问你有没有使用C99的math.h signbit函数（结合x == 0），而是使用1.0 / x < 0检查符号的原因零？

@mctylr：不记得了。 signbit 似乎也是work。

【参考方案10】：

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

【参考方案11】：

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

【参考方案12】：

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

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

【讨论】：

只为遇到此问题的人提供：C++17 的 <algorithm> 标头引入了 std::clamp(n, low, high [, compare])。

【参考方案13】：

如果我理解正确，您希望将值“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;

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

【讨论】：

即使将MY_MIN,MY_MAX 用作int，如果a 不在int 范围附近，这种方法也会失败，因为(int)a 是一个问题。

【参考方案14】：

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

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

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

fabs(x)

而另一个反映1.0左右为

1.0-fabs(x-1.0)

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

【讨论】：

很多精度损失，从大约x < 0.25开始。使用values < DBL__EPSILON，结果会丢失所有精度。