C ++中的高效整数floor函数

Efficient integer floor function in C++

我想定义一个高效的整型 floor 函数,即从 float 或 double 执行向负无穷大截断的转换。

我们可以假设这些值不会发生整数溢出。到目前为止我有几个选择

转换为 int 是出了名的慢。 if 测试也是如此。 floor功能我没有计时过,但是看到帖子说它也很慢

你能想到在速度、准确性或允许范围方面更好的替代方案吗?它不需要便携。目标是最近的 x86/x64 架构。

看看magic numbers。网页上提出的算法应该远比简单的转换更高效。我自己从未使用过它,但这是他们在网站上提供的性能比较(xs_ToInt 和 xs_CRoundToInt 是建议的功能):

Performing 10000000 times:
simple cast           2819 ms i.e. i = (long)f;
xs_ToInt              1242 ms i.e. i = xs_ToInt(f); //numerically same as above
bit-twiddle(full)     1093 ms i.e. i = BitConvertToInt(f); //rounding from Fluid
fistp                  676 ms i.e. i = FISTToInt(f); //Herf, et al x86 Assembly rounding 
bit-twiddle(limited)   623 ms i.e. i = FloatTo23Bits(f); //Herf, rounding only in the range (0...1]  
xs_CRoundToInt         609 ms i.e. i = xs_CRoundToInt(f); //rounding with "magic" numbers

此外,xs_ToInt 显然已修改,从而提高了性能:

Performing 10000000 times:
simple cast convert   3186 ms i.e. fi = (f*65536);
fistp convert         3031 ms i.e. fi = FISTToInt(f*65536);
xs_ToFix               622 ms i.e. fi = xs_Fix<16>::ToFix(f);

简要说明 'magic numbers' 方法的工作原理:

"Basically, in order to add two floating point numbers, your processor "lines up" the decimal points of the numbers so that it can easily add the bits. It does this by "normalizing" the numbers such that the most significant bits are preserved, i.e. the smaller number "normalizes" to match the bigger one. So the principle of the "magic number" conversion that xs_CRoundToInt() uses is this: We add a big enough floating point number (a number that is so big that there are significant digits only UP TO the decimal point, and none after it) to the one you're converting such that: (a) the number gets normalized by the processor to its integer equivalent and (b) adding the two does not erase the integral significat bits in the number you were trying to convert (i.e. XX00 + 00YY = XXYY)."

引用自同一网页。

如果您批量执行此操作,如果您知道自己在做什么,编译器可能会对其进行自动矢量化。例如,这是一个小的实现,它在 GCC 上将浮点数自动向量化为整数:

#include <cmath>

// Compile with -O3 and -march=native to see autovectorization
__attribute__((optimize("-fno-trapping-math")))
void testFunction(float* input, int* output, int length) {
  // Assume the input and output are aligned on a 32-bit boundary.
  // Of course, you have  to ensure this when calling testFunction, or else
  // you will have problems.
  input = static_cast<float*>(__builtin_assume_aligned(input, 32));
  output = static_cast<int*>(__builtin_assume_aligned(output, 32));

  // Also assume the length is a multiple of 32.
  if (length & 31) __builtin_unreachable();

  // Do the conversion
  for (int i = 0; i < length; ++i) {
    output[i] = floor(input[i]);
  }
}

这是为 x86-64 生成的程序集(使用 AVX512 指令):

testFunction(float*, int*, int):
        test    edx, edx
        jle     .L5
        lea     ecx, [rdx-1]
        xor     eax, eax
.L3:
        # you can see here that the conversion was vectorized
        # to a vrndscaleps (that will round the float appropriately)
        # and a vcvttps2dq (thal will perform the conversion)
        vrndscaleps     ymm0, YMMWORD PTR [rdi+rax], 1
        vcvttps2dq      ymm0, ymm0
        vmovdqa64       YMMWORD PTR [rsi+rax], ymm0
        add     rax, 32
        cmp     rax, rdx
        jne     .L3
        vzeroupper
.L5:
        ret

如果您的目标不支持 AVX512,它仍会使用 SSE4.1 指令进行自动矢量化,前提是您有这些指令。这是 -O3 -msse4.1:

的输出
testFunction(float*, int*, int):
        test    edx, edx
        jle     .L1
        shr     edx, 2
        xor     eax, eax
        sal     rdx, 4
.L3:
        roundps xmm0, XMMWORD PTR [rdi+rax], 1
        cvttps2dq       xmm0, xmm0
        movaps  XMMWORD PTR [rsi+rax], xmm0
        add     rax, 16
        cmp     rax, rdx
        jne     .L3
.L1:
        ret

See it live on godbolt

Casting to int is notoriously slow.

也许您从 x86-64 开始就一直生活在岩石下,或者以其他方式错过了 x86 上一段时间以来并非如此。 :)

SSE/SSE2 有一条指令将 转换为 t运行cation(而不是默认的舍入模式)。 ISA 有效地支持此操作正是因为使用 C 语义的转换在实际代码库中并不少见。 x86-64 代码使用 SSE/SSE2 XMM 寄存器用于标量 FP 数学,而不是 x87,因为这个和其他使其更高效的东西。即使是现代 32 位代码也使用 XMM 寄存器进行标量数学计算。

为 x87(没有 SSE3 fisttp)编译时,编译器过去必须将 x87 舍入模式更改为 t运行cation,FP 存储到内存,然后再次将舍入模式更改回来. (然后从内存中重新加载整数,通常是从堆栈上的本地,如果用它做进一步的事情。)x87 在这方面 很糟糕

是的, 非常慢,例如在 2006 年写下@Kirjain 的回答中的 link 时,如果您仍然有 32 位 CPU 或者正在使用 x86-64 CPU 到 运行 32 位代码。


不直接支持使用 t运行cation 或默认(最近)以外的舍入模式进行转换,在 SSE4.1 roundps/roundpd 之前,您最好的选择是magic-number @Kirjain 的回答中 the 2006 link 中的技巧。

那里有一些不错的技巧,但仅适用于 double -> 32 位整数。如果您有 float.

,不太可能值得扩展到 double

或者更常见的是,只需添加一个 large-magnitude 数字即可触发舍入,然后再次减去它以返回到原始范围。这可以用于 float 而无需扩展到 double,但我不确定让 floor 工作有多容易。


无论如何,这里明显的解决方案是_mm256_floor_ps()_mm256_cvtps_epi32vroundpsvcvtps2dq)。 non-AVX 版本可以与 SSE4.1 一起使用。

我不确定我们是否可以做得更好;如果你有一个巨大的数组要处理(并且无法设法将这项工作与其他工作交错),你可以将 MXCSR 舍入模式设置为 "towards -Inf"(floor)并简单地使用 vcvtps2dq(使用当前舍入模式)。然后把它放回去。但最好是 cache-block 您的转换或在生成数据时即时进行转换,大概来自其他需要将 FP 舍入模式设置为默认最近的 FP 计算。

roundps/pd/ss/sd 在 Intel CPUs 上是 2 微指令,但在 AMD Ryzen 上只有 1 微指令(每 128 位通道)。 cvtps2dq 也是 1 微指令。 packed double->int 转换还包括 shuffle。标量 FP->int 转换(复制到整数寄存器)通常也需要额外的 uop。

所以在某些情况下 magic-number 技巧有可能获胜;如果 _mm256_floor_ps() + cvt 是关键瓶颈的一部分(或者如果你有 double 并且想要 int32 更有可能)。


@Cássio Renan 的 int foo = floorf(f) 实际上会 auto-vectorize 如果使用 gcc -O3 -fno-trapping-math(或 -ffast-math)编译,-march= 具有 SSE4.1 或AVX。 https://godbolt.org/z/ae_KPv

如果您将其与其他未手动向量化的标量代码一起使用,这可能会有用。特别是如果您希望编译器 auto-vectorize 整个事情。

为什么不直接使用这个:

#include <cmath>

auto floor_(float const x) noexcept
{
  int const t(x);

  return t - (t > x);
}

这里是对 Cássio Renan 的精彩回答的修改。它用标准 C++ 替换了所有 compiler-specific 扩展,并且在理论上可以移植到任何符合标准的编译器。此外,它检查参数是否正确对齐,而不是假设如此。它优化了相同的代码。

#include <assert.h>
#include <cmath>
#include <stddef.h>
#include <stdint.h>

#define ALIGNMENT alignof(max_align_t)
using std::floor;

// Compiled with: -std=c++17 -Wall -Wextra -Wpedantic -Wconversion -fno-trapping-math -O -march=cannonlake -mprefer-vector-width=512

void testFunction(const float in[], int32_t out[], const ptrdiff_t length)
{
  static_assert(sizeof(float) == sizeof(int32_t), "");
  assert((uintptr_t)(void*)in % ALIGNMENT == 0);
  assert((uintptr_t)(void*)out % ALIGNMENT == 0);
  assert((size_t)length % (ALIGNMENT/sizeof(int32_t)) == 0);

  alignas(ALIGNMENT) const float* const input = in;
  alignas(ALIGNMENT) int32_t* const output = out;

  // Do the conversion
  for (int i = 0; i < length; ++i) {
    output[i] = static_cast<int32_t>(floor(input[i]));
  }
}

这在 GCC 上的优化效果不如使用 non-portable 扩展的原始版本。 C++ 标准确实支持 alignas 说明符、对对齐数组的引用以及 std::align 函数 returns 缓冲区内的对齐范围。 None 然而,这些使得我测试的任何编译器生成对齐而不是未对齐的向量加载和存储。

虽然 alignof(max_align_t) 在 x86_64 上只有 16,并且可以将 ALIGNMENT 定义为常量 64,但这无助于任何编译器生成更好的代码,所以我追求便携性。强制编译器假定指针对齐的最接近可移植方式的方法是使用 <immintrin.h> 中的类型,大多数 x86 编译器都支持该类型,或者用 [=13] 定义 struct =] 说明符。通过检查预定义的宏,您还可以将宏扩展为 Linux 编译器上的 __attribute__ ((aligned (ALIGNMENT))),或 Windows 编译器上的 __declspec (align (ALIGNMENT)),以及我们不知道的编译器上的安全内容关于,但是 GCC 需要 type 上的属性来实际生成对齐的加载和存储。

此外,原始示例调用了 bulit-in 来告诉 GCC length 不可能不是 32 的倍数。如果您 assert() 这个或调用标准abort()这样的函数,GCC、Clang、ICC都不会做同样的推导。因此,他们生成的大部分代码将处理 length 不是向量宽度的整数倍的情况。

一个可能的原因是,这两种优化都无法让您获得那么快的速度:具有对齐地址的未对齐内存指令在 Intel CPU 上速度很快,并且处理 length 不是很好的回合的情况的代码number 有几个字节长,运行s 在常数时间内。

作为脚注,GCC 能够优化 <cmath> 中的内联函数,优于 <math.c> 中实现的宏。

GCC 9.1 需要一组特定的选项来生成 AVX512 代码。默认情况下,即使使用 -march=cannonlake,它也会更喜欢 256 位向量。它需要 -mprefer-vector-width=512 来生成 512 位代码。 (感谢 Peter Cordes 指出这一点。)它用展开的代码跟进矢量化循环以转换数组的任何剩余元素。

这是矢量化的主循环,减去一些 constant-time 初始化,error-checking 和 clean-up 代码只会 运行 一次:

.L7:
        vrndscaleps     zmm0, ZMMWORD PTR [rdi+rax], 1
        vcvttps2dq      zmm0, zmm0
        vmovdqu32       ZMMWORD PTR [rsi+rax], zmm0
        add     rax, 64
        cmp     rax, rcx
        jne     .L7

eagle-eyed 会注意到与 Cássio Renan 程序生成的代码有两个不同之处:它使用 %zmm 而不是 %ymm 寄存器,并且它使用未对齐的 vmovdqu32 而不是对齐 vmovdqa64.

具有相同标志的 Clang 8.0.0 对展开循环做出了不同的选择。每次迭代都对 8 个 512 位向量(即 128 single-precision 浮点数)进行操作,但拾取剩菜的代码并未展开。如果在那之后至少剩下 64 个浮点数,它会为这些使用另外四个 AVX512 指令,然后使用非矢量化循环清除任何额外的部分。

如果你用 Clang++ 编译原始程序,它会毫无怨言地接受它,但不会做同样的优化:它仍然不会假设 length 是矢量宽度的倍数,也不是指针对齐。

它更喜欢 AVX512 代码而不是 AVX256,即使没有 -mprefer-vector-width=512

        test    rdx, rdx
        jle     .LBB0_14
        cmp     rdx, 63
        ja      .LBB0_6
        xor     eax, eax
        jmp     .LBB0_13
.LBB0_6:
        mov     rax, rdx
        and     rax, -64
        lea     r9, [rax - 64]
        mov     r10, r9
        shr     r10, 6
        add     r10, 1
        mov     r8d, r10d
        and     r8d, 1
        test    r9, r9
        je      .LBB0_7
        mov     ecx, 1
        sub     rcx, r10
        lea     r9, [r8 + rcx]
        add     r9, -1
        xor     ecx, ecx
.LBB0_9:                                # =>This Inner Loop Header: Depth=1
        vrndscaleps     zmm0, zmmword ptr [rdi + 4*rcx], 9
        vrndscaleps     zmm1, zmmword ptr [rdi + 4*rcx + 64], 9
        vrndscaleps     zmm2, zmmword ptr [rdi + 4*rcx + 128], 9
        vrndscaleps     zmm3, zmmword ptr [rdi + 4*rcx + 192], 9
        vcvttps2dq      zmm0, zmm0
        vcvttps2dq      zmm1, zmm1
        vcvttps2dq      zmm2, zmm2
        vmovups zmmword ptr [rsi + 4*rcx], zmm0
        vmovups zmmword ptr [rsi + 4*rcx + 64], zmm1
        vmovups zmmword ptr [rsi + 4*rcx + 128], zmm2
        vcvttps2dq      zmm0, zmm3
        vmovups zmmword ptr [rsi + 4*rcx + 192], zmm0
        vrndscaleps     zmm0, zmmword ptr [rdi + 4*rcx + 256], 9
        vrndscaleps     zmm1, zmmword ptr [rdi + 4*rcx + 320], 9
        vrndscaleps     zmm2, zmmword ptr [rdi + 4*rcx + 384], 9
        vrndscaleps     zmm3, zmmword ptr [rdi + 4*rcx + 448], 9
        vcvttps2dq      zmm0, zmm0
        vcvttps2dq      zmm1, zmm1
        vcvttps2dq      zmm2, zmm2
        vcvttps2dq      zmm3, zmm3
        vmovups zmmword ptr [rsi + 4*rcx + 256], zmm0
        vmovups zmmword ptr [rsi + 4*rcx + 320], zmm1
        vmovups zmmword ptr [rsi + 4*rcx + 384], zmm2
        vmovups zmmword ptr [rsi + 4*rcx + 448], zmm3
        sub     rcx, -128
        add     r9, 2
        jne     .LBB0_9
        test    r8, r8
        je      .LBB0_12
.LBB0_11:
        vrndscaleps     zmm0, zmmword ptr [rdi + 4*rcx], 9
        vrndscaleps     zmm1, zmmword ptr [rdi + 4*rcx + 64], 9
        vrndscaleps     zmm2, zmmword ptr [rdi + 4*rcx + 128], 9
        vrndscaleps     zmm3, zmmword ptr [rdi + 4*rcx + 192], 9
        vcvttps2dq      zmm0, zmm0
        vcvttps2dq      zmm1, zmm1
        vcvttps2dq      zmm2, zmm2
        vcvttps2dq      zmm3, zmm3
        vmovups zmmword ptr [rsi + 4*rcx], zmm0
        vmovups zmmword ptr [rsi + 4*rcx + 64], zmm1
        vmovups zmmword ptr [rsi + 4*rcx + 128], zmm2
        vmovups zmmword ptr [rsi + 4*rcx + 192], zmm3
.LBB0_12:
        cmp     rax, rdx
        je      .LBB0_14
.LBB0_13:                               # =>This Inner Loop Header: Depth=1
        vmovss  xmm0, dword ptr [rdi + 4*rax] # xmm0 = mem[0],zero,zero,zero
        vroundss        xmm0, xmm0, xmm0, 9
        vcvttss2si      ecx, xmm0
        mov     dword ptr [rsi + 4*rax], ecx
        add     rax, 1
        cmp     rdx, rax
        jne     .LBB0_13
.LBB0_14:
        pop     rax
        vzeroupper
        ret
.LBB0_7:
        xor     ecx, ecx
        test    r8, r8
        jne     .LBB0_11
        jmp     .LBB0_12

ICC 19 也生成 AVX512 指令,但与 clang 有很大不同。它使用魔术常数做了更多 set-up,但不展开任何循环,而是在 512 位向量上运行。

此代码也适用于其他编译器和体系结构。 (尽管 MSVC 仅支持 AVX2 及以下的 ISA,不能 auto-vectorize 循环。)例如,在具有 -march=armv8-a+simd 的 ARM 上,它会生成具有 frintm v0.4s, v0.4sfcvtzs v0.4s, v0.4s 的矢量化循环.

Try it for yourself.