为什么floor()函数执行速度如此缓慢？

Question

最近我写了一些代码（ISO/ANSI C），然而它的性能表现令人惊讶。长话短说，罪魁祸首是 floor() 函数。不仅速度慢，而且在使用英特尔编译器（即 ICL）进行向量化时也无法实现。

以下是对在二维矩阵中的所有单元格执行 floor 的一些基准测试：

VC:  0.10
ICL: 0.20

与简单的强制类型转换进行比较：

VC:  0.04
ICL: 0.04

floor() 比一个简单的转换要慢这么多，为什么？它本质上做的是相同的事情(除了负数)。 第二个问题：有人知道一个超快的 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]);
        }
    }
}

- Roger

Answer 1

几个因素使得floor函数比一个转换显式的整数类型更慢，同时也阻止了向量化。

其中最重要的因素是：

floor函数可以修改全局状态。如果你传递一个太大无法表示为整数的浮点值，那么errno变量会被设置为EDOM。还会进行NaNs的特殊处理。所有这些行为都是为了应用程序能检测到溢出情况并以某种方式处理该情况（别问我怎么处理）。

检测这些问题不简单，并且占据floor函数执行时间的90%以上。实际的舍入操作是廉价的，可以被内联/向量化。由于代码量很大，所以完全内联floor函数会使你的程序运行更慢。

一些编译器有特殊的编译器标志，允许编译器优化掉一些很少使用的C标准规则。例如，GCC 可以告诉编译器你根本不关心errno。为了做到这一点，请传递-fno-math-errno或-ffast-math。ICC和VC可能有类似的编译器标志。

顺便说一下 - 您可以使用简单的转换方式编写自己的floor函数。只是您需要分别处理负数和正数的情况。如果您不需要处理溢出和NaNs的特殊情况，那么这可能会更快。

- Nils Pipenbrinck

Answer 2

如果你要将 floor() 操作的结果转换为 int 类型，并且不担心溢出问题，那么以下代码比 (int)floor(x) 更快：

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

- dgobbi

---

如果你想把这个函数放到头文件中，那么最好使用 static inline 而不是 inline - 参见 https://dev59.com/HWkv5IYBdhLWcg3w8FSY#10245969 - Christoph

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

或者使用 floor

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

- PolarBear2015

---

嗯，floor 对于负整数会给出错误的答案，而 ceil 对于正整数会给出错误的答案。 - geometrian

谢谢imallett。代码现在应该没问题了。 - PolarBear2015

Answer 4

现代x86 CPU上实际最快的实现方式，对于一个大数组来说，是：

• 将MXCSR FP舍入模式更改为向下舍入（也称为floor）。在C语言中，可以通过fenv或_mm_getcsr/_mm_setcsr来实现。

• 循环遍历数组，在SIMD向量上执行_mm_cvtps_epi32，使用当前的舍入模式将4个float转换为32位整数。(并将结果向量存储到目标位置。)

  cvtps2dq xmm0，[rdi]是自K10或Core 2以来任何Intel或AMD CPU上的单条微融合uop。(https://agner.org/optimize/)256位AVX版本也是如此，具有YMM向量。

• 将当前舍入模式恢复为正常的IEEE默认模式，使用MXCSR的原始值。(四舍五入，平均数作为打结点)

这使得每个时钟周期可以加载、转换和存储1个SIMD向量的结果，与截断一样快。 （SSE2具有用于截断的特殊FP-> int转换指令，正是因为它非常常见于C编译器。在旧的x87时代，即使(int)x也需要更改x87舍入模式为截断，然后再改回来。cvttps2dq用于带截断的打包浮点数->整数（请注意助记符中额外的t）。或者对于标量，从XMM到整数寄存器，cvttss2si或cvttsd2si用于标量double到标量整数。

通过一些循环展开和/或良好的优化，这应该是可能的，而不会在前端造成瓶颈，只要假设没有缓存未命中的瓶颈，每个时钟周期的存储吞吐量为1。 （在Skylake之前的英特尔上，还受到每时钟周期打包转换吞吐量为1的瓶颈的限制。）即每个周期16、32或64字节，使用SSE2、AVX或AVX512。

---

如果不改变当前的舍入模式，您需要使用SSE4.1 roundps 将float四舍五入为最近的整数float，并选择所需的舍入模式。或者，您可以使用其他答案中显示的技巧来处理大小足够小以适合有符号32位整数的浮点数，因为那是您最终的目标格式。

---

（使用正确的编译器选项，例如-fno-math-errno和正确的-march或-msse4选项，编译器可以使用roundps或标量和/或双精度等效的roundsd xmm1，xmm0，1内联floor，但这将花费2个uops，并且Haswell的标量或矢量每2个时钟周期具有1个吞吐量。实际上，即使没有任何快速数学选项，gcc8.2也会在-march=haswell的情况下为floor内联roundsd，正如您在Godbolt编译器浏览器上看到的那样。但是这不是x86-64的基线，因此如果您的机器支持它，则需要启用它。）

- Peter Cordes

---

附注：某种程度上，icc似乎不知道vcvtps2dq指令依赖于MXCSR控制和状态寄存器的值。在这个例子中，icc交换了x=_mm_cvtps_epi32(y);和_MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST);的顺序。 - wim

@wim：是的，我也想知道这是否会成为一个问题。如果#pragma STDC FENV_ACCESS ON对任何实际编译器有效，我应该添加一些内容。 （C++11及更高版本中是否存在FENV_ACCESS pragma？）。或者尝试ICC编译选项，例如-fp-model strict，告诉它您修改了FP舍入模式。（ICC默认值为-fp-model fast=1。） - Peter Cordes

Answer 5

是的，floor()在所有平台上都非常慢，因为它必须实现IEEE fp规范中的许多行为。您无法在内部循环中使用它。

我有时会使用一个宏来近似floor():

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

PSEUDO_FLOOR()和floor()的行为并不完全相同：例如，floor(-1) == -1，但是PSEUDO_FLOOR(-1) == -2，但对于大多数用途来说已经足够接近了。

- jrgc

---

天真的实现。PSEUDO_FLOOR(x++) 会破坏这个实现。 - user47322

是的，Charlie。最好将其制作为内联函数。 - Triang3l

Answer 6

一个实际上不需要分支的版本，只需要在浮点数和整数域之间进行单个转换，将值 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;
}

- Aki Suihkonen

---

这是否依赖于“long”比“int”更宽以正确处理完整范围的“int”结果？这在许多32位平台上都不是这种情况，例如在x86-64 Windows（其中int和long均为32位的LLP64 ABI）。因此，也许你应该使用“long long”。但仍然是一个好主意。 - Peter Cordes

是的（也就是长整型比整型宽），但我认为这可以通过强制转换为无符号整型来缓解。 - Aki Suihkonen

在x86上，将double转换为unsigned long有些慢。https://godbolt.org/z/1UqaQw。直到AVX512，x86-64才有这样的指令，只能将`double`转换为有符号整数。在32位x86上，其中`unsigned long是32位类型，x87 fistp可以执行FP-> 64位有符号整数，您可以使用其低半部分作为unsigned int。但截断需要SSE3 fisttp`或更改舍入模式。SSE2也无法将截断为32位无符号整数或64位有符号整数。其他答案可能更有效率。 - Peter Cordes

Answer 7

1. 它们并不做同样的事情。floor()是一个函数。因此，使用它会产生一个函数调用，分配一个堆栈帧，复制参数并检索结果。 强制转换不是一个函数调用，所以它使用更快的机制（我相信它可能使用寄存器来处理值）。
2. 可能floor()已经被优化了。
3. 你能从你的算法中挤出更多的性能吗？也许交换行和列可以帮助？你能缓存常见的值吗？你的编译器的所有优化都开启了吗？你能切换操作系统吗？编译器？ Jon Bentley's Programming Pearls有一个关于可能的优化的很好的评论。

- Yuval F

---

永远不要假设标准库已经被优化了。它们几乎总是非常慢的。有时候，通过使用自己的定制代码，您可以获得重大的速度提升。 - Joe

floor()是一个函数，但它被广泛使用，以至于编译器将其视为内置函数，例如memcpy或sqrt，并在需要时进行内联。例如，GCC -O2用于x86-64时，即使需要多个指令，也会将其内联，而没有SSE4.1的情况下，无法进行roundss / roundps（https://godbolt.org/z/5jdTvcx7x）。但是，没有SSE4.1，与具有截断的fp->int相比，速度要慢得多，后者具有更快的硬件支持。 - Peter Cordes

Answer 8

快速双精度取整

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

---

终端日志

测试自定义_1 8.3837

测试本地_1 18.4989

测试自定义_2 8.36333

测试本地_2 18.5001

测试自定义_3 8.37316

测试本地_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倍。

- Jan Cajthaml

---

你的 round() 函数不起作用。你需要使用浮点数取模来检查小数部分是否大于 0.5，或者你可以使用旧的 (int) (double_value + 0.5) 技巧来进行四舍五入。 - Jason R

关于四舍五入的FP->int，请参阅https://dev59.com/rXRB5IYBdhLWcg3w3K4J#47347224。 - Peter Cordes