我想定义一个有效的整数底函数,即从浮点数或双精度执行截断到负无穷大的转换。
我们可以假设这些值不会发生整数溢出。到目前为止,我有几个选择
-
转换为 int; 这需要对负值进行特殊处理,因为强制转换会截断到零;
I= int(F); if (I < 0 && I != F) I--;
-
将地板的结果转换为 int;
int(floor(F));
-
以较大的偏移转换为 int 以获得正数(对于较大的值,这可能会返回错误的结果);
int(F + double(0x7fffffff)) - 0x7fffffff;
众所周知,转换为 int 的速度非常慢。如果测试也是如此。我没有计时地板功能,但看到帖子声称它也很慢。
你能想到在速度、精度或允许范围方面更好的替代方案吗?它不需要是便携式的。目标是最新的 x86/x64 体系结构。
转换为 int 是出了名的慢。
也许你从 x86-64 开始就一直生活在岩石下,或者错过了在 x86 上已经有一段时间不是这样了。 :)
SSE/SSE2 有一个使用截断进行转换的指令(而不是默认的舍入模式)。 ISA 之所以有效地支持此操作,正是因为使用 C 语义进行转换在实际代码库中并不罕见。 x86-64 代码使用 SSE/SSE2 XMM 寄存器进行标量 FP 数学运算,而不是 x87,因为这个和其他使其更有效率的东西。 即使是现代 32 位代码也使用 XMM 寄存器进行标量数学运算。
当编译为 x87(没有 SSE3fisttp
)时,编译器过去必须将 x87 舍入模式更改为截断,FP 存储到内存,然后再次更改舍入模式。 (然后从内存中重新加载整数,通常是从堆栈上的本地重新加载,如果用它做进一步的事情。 x87 为此很糟糕。
是的,这非常慢,例如,在 2006 年,当您仍然拥有 32 位 CPU 或使用 x86-64 CPU 运行 32 位代码时,@Kirjain答案中的链接被写入时。
不直接支持使用截断或默认值(最近)以外的舍入模式进行转换,直到 SSE4.1roundps
/roundpd
您最好的选择是魔术数字技巧,例如 @Kirjain 答案中的 2006 年链接。
那里有一些不错的技巧,但仅适用于double
-> 32 位整数。 如果你有float
,不太可能值得扩展到double
.
或者更常见的是,只需添加一个大数量级数字来触发舍入,然后再次减去它以返回到原始范围。 这可以适用于float
而无需扩展到double
,但我不确定使floor
工作有多容易。
无论如何,这里明显的解决方案是_mm256_floor_ps()
和_mm256_cvtps_epi32
(vroundps
和vcvtps2dq
)。 非 AVX 版本可以与 SSE4.1 一起使用。
我不确定我们是否可以做得更好;如果你有一个巨大的数组要处理(并且无法设法将这项工作与其他工作交错),您可以将MXCSR舍入模式设置为"朝向-Inf"(地板),只需使用vcvtps2dq
(使用当前的舍入模式)。 然后将其放回原处。 但是,最好缓存阻止转换或在生成数据时动态进行转换,大概来自需要将 FP 舍入模式设置为默认"最近"的其他 FP 计算。
roundps
/pd/ss/sd 在英特尔 CPU 上为 2 uops,但在 AMD 锐龙上仅为 1 uop(每 128 位通道)。cvtps2dq
也是 1 uop。 打包的双>整数转换还包括随机播放。 标量FP->int转换(复制到整数寄存器)通常还需要额外的uop。
因此,在某些情况下,魔术数字技巧有可能获胜;也许值得研究_mm256_floor_ps()
+ cvt是否是关键瓶颈的一部分(或者如果你有双精度并且想要int32,则更有可能)。
@Cássio 如果用gcc -O3 -fno-trapping-math
(或-ffast-math
)编译,Renan 的int foo = floorf(f)
实际上会自动矢量化,-march=
具有 SSE4.1 或 AVX 的东西。 https://godbolt.org/z/ae_KPv
如果您将其与其他未手动矢量化的标量代码一起使用,这可能很有用。 特别是如果你希望编译器会自动矢量化整个事情。
看看幻数。网页上提出的算法应该比简单的投射效率高得多。我自己从未使用过它,但这是他们在网站上提供的性能比较(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);
简要说明"幻数"方法的工作原理:
">基本上,为了添加两个浮点数,您的处理器将数字的小数点"排列"起来,以便它可以轻松地添加位。它通过"规范化"数字来实现这一点,以便保留最重要的位,即较小的数字"规范化"以匹配较大的数字。因此,xs_CRoundToInt() 使用的"幻数"转换的原理是这样的:我们将一个足够大的浮点数(一个非常大的数字,以至于只有小数点以下有有效数字,而在它之后没有)添加到您要转换的数字中,以便:(a) 处理器将数字规范化为其整数等价物,并且 (b) 将两者相加不会擦除您尝试转换的数字(即 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
在Godbolt上观看直播
为什么不直接使用它:
#include <cmath>
auto floor_(float const x) noexcept
{
int const t(x);
return t - (t > x);
}
这是对卡西奥·雷南(Cássio Renan)出色答案的修改。 它用标准C++替换所有特定于编译器的扩展,理论上可以移植到任何符合标准的编译器。 此外,它还检查参数是否正确对齐,而不是假设正确对齐。 它针对相同的代码进行优化。
#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上的优化效果不如使用非可移植扩展的原始版本。 C++ 标准确实支持alignas
说明符、对对齐数组的引用以及返回缓冲区内对齐范围的std::align
函数。 但是,这些都不会使我测试的任何编译器生成对齐而不是未对齐的向量加载和存储。
尽管alignof(max_align_t)
在x86_64上只有 16,并且可以将ALIGNMENT
定义为常量 64,但这无助于任何编译器生成更好的代码,所以我选择了可移植性。 强制编译器假定 poitner 对齐的最接近可移植方法的方法是使用<immintrin.h>
中的类型,大多数 x86 编译器都支持这些类型,或者使用alignas
说明符定义struct
。 通过检查预定义的宏,您还可以扩展宏以在 Linux 编译器上__attribute__ ((aligned (ALIGNMENT)))
,或在 Windows 编译器上__declspec (align (ALIGNMENT))
宏,以及在我们不知道的编译器上安全的东西,但 GCC 需要类型上的属性来实际生成对齐的加载和存储。
此外,原始示例称为 bulit-in,以告诉 GCC,length
不可能不是 32 的倍数。 如果您assert()
这个或调用标准函数,例如abort()
,GCC,Clang和ICC都不会进行相同的扣除。 因此,他们生成的大多数代码将处理length
不是矢量宽度的好舍入倍数的情况。
造成这种情况的一个可能原因是,这两种优化都没有为您提供那么快的速度:在英特尔 CPU 上,具有对齐地址的未对齐内存指令速度很快,并且处理length
不是一个好整数的情况的代码只有几个字节长,并且在恒定时间内运行。
作为脚注,GCC 能够从<cmath>
优化内联函数,比<math.c>
中实现的宏更好。
GCC 9.1 需要一组特定的选项来生成 AVX512 代码。 默认情况下,即使使用-march=cannonlake
,它也更喜欢 256 位向量。 它需要-mprefer-vector-width=512
来生成 512 位代码。 (感谢Peter Cordes指出这一点。 它使用展开的代码跟踪矢量化循环,以转换数组的任何剩余元素。
下面是矢量化的主循环,减去一些只运行一次的常量时间初始化、错误检查和清理代码:
.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
眼尖的人会注意到与 Cássio Renan 程序生成的代码有两个不同之处:它使用 %zmm 而不是 %ymm 寄存器,并且它使用未对齐的vmovdqu32
而不是对齐的vmovdqa64
存储结果。
具有相同标志的 Clang 8.0.0 对展开循环做出了不同的选择。 每次迭代都对 8 个 512 位向量(即 128 个单精度浮点数)进行操作,但拾取剩余向量的代码不会展开。 如果之后至少剩下 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
有很大不同。 它使用魔术常量进行更多设置,但不展开任何循环,而是在 512 位向量上运行。
此代码也适用于其他编译器和体系结构。 (尽管 MSVC 仅支持 AVX2 以下的 ISA,并且无法自动矢量化环路。 例如,在带有-march=armv8-a+simd
的ARM上,它会生成一个带有frintm v0.4s, v0.4s
和fcvtzs v0.4s, v0.4s
的矢量化循环。
亲自尝试一下。