为什么 floor() 这么慢？

Question

我最近编写了一些代码 (ISO/ANSI C)，但对它实现的糟糕性能感到惊讶。长话短说，原来罪魁祸首是floor() 函数。它不仅速度慢，而且没有矢量化（使用英特尔编译器，也就是 ICL）。

以下是对 2D 矩阵中的所有单元格执行地板的一些基准：

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

Answer 1

有几件事会使 floor 比 cast 慢并阻止矢量化。

最重要的一个：

floor 可以修改全局状态。如果您传递的值太大而无法以浮点格式表示为整数，则 errno 变量将设置为 EDOM。还对 NaN 进行了特殊处理。所有这些行为都适用于想要检测溢出情况并以某种方式处理这种情况的应用程序（不要问我如何）。

检测这些有问题的条件并不简单，占地板执行时间的 90% 以上。实际的舍入很便宜，可以内联/矢量化。还有很多代码，所以内联整个 floor-function 会让你的程序运行得更慢。

一些编译器具有特殊的编译器标志，允许编译器优化掉一些很少使用的 c 标准规则。例如 GCC 可以被告知您对 errno 根本不感兴趣。为此，请通过 -fno-math-errno 或 -ffast-math。 ICC 和 VC 可能有类似的编译器标志。

顺便说一句-您可以使用简单的演员来滚动自己的地板功能。你只需要以不同的方式处理消极和积极的情况。如果您不需要对溢出和 NaN 进行特殊处理，这可能会快很多。

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

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

或使用地板

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

Answer 4

大型数组在现代 x86 CPU 上的实际最快实现是

• 将 MXCSR FP 舍入模式更改为向 -Infinity（又名 floor）舍入。在 C 中，这应该可以通过 fenv 或 _mm_getcsr / _mm_setcsr 实现。

• 循环遍历数组，对 SIMD 向量执行 _mm_cvtps_epi32，使用当前舍入模式将 4 个 floats 转换为 32 位整数。 （并将结果向量存储到目的地。）

  cvtps2dq xmm0, [rdi] 是自 K10 或 Core 2 以来任何 Intel 或 AMD CPU 上的单个微融合微指令。(https://agner.org/optimize/) 与 256 位 AVX 版本相同，带有 YMM 向量。

• 使用 MXCSR 的原始值将当前舍入模式恢复为正常的 IEEE 默认模式。 （四舍五入，甚至作为抢七）

这允许在每个时钟周期加载 + 转换 + 存储 1 个 SIMD 结果向量，与截断一样快。 （SSE2 有一个特殊的 FP->int 转换指令用于截断，正是因为 C 编译器非常需要它。在 x87 的糟糕过去，即使(int)x 也需要将 x87 舍入模式更改为截断然后返回。@987654323 @（注意助记符中额外的 t）。或者对于标量，从 XMM 到整数寄存器，cvttss2si 或 cvttsd2si 对于标量 double 到标量整数。

通过一些循环展开和/或良好的优化，这应该可以在前端没有瓶颈的情况下实现，假设没有缓存未命中瓶颈，每时钟只有 1 个存储吞吐量。 （在 Skylake 之前的 Intel 上，每时钟 1 个打包转换吞吐量也存在瓶颈。）即 每个周期 16、32 或 64 字节，使用 SSE2、AVX 或 AVX512。

---

在不更改当前舍入模式的情况下，您需要 SSE4.1 roundps 使用您选择的舍入模式将 float 舍入到最接近的整数 float。或者您可以使用其他答案中显示的技巧之一，这些技巧适用于幅度足够小的浮点数以适合有符号的 32 位整数，因为无论如何这都是您的最终目标格式。）

---

（使用正确的编译器选项，例如-fno-math-errno，以及正确的-march 或-msse4 选项，编译器可以使用roundps 或标量和/或双精度等价物内联floor，例如roundsd xmm1, xmm0, 1，但这需要 2 微秒，并且在 Haswell 上每 2 个时钟吞吐量为 1 个用于标量或向量。实际上，gcc8.2 将内联 roundsd 用于 floor，即使没有任何快速数学选项 as you can see on the Godbolt compiler explorer。但是那是-march=haswell。不幸的是它不是x86-64的基线，所以如果你的机器支持它，你需要启用它。）

Answer 5

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

我有时会使用宏来近似 floor()：

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

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

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

正如 cmets 中所指出的，此实现依赖于临时值 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;
}

Answer 7

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

Answer 8

快速双轮

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

---

终端日志

测试 custom_1 8.3837

测试 native_1 18.4989

测试 custom_2 8.36333

测试 native_2 18.5001

测试 custom_3 8.37316

测试 native_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 倍。