更快地近似数组的倒数平方根
Faster approximate reciprocal square root of an array
如何使用 popcnt 和 SSE4.2 在 cpu 上更快地计算数组的近似倒数平方根？
输入是存储在浮点数组中的正整数（范围从 0 到大约 200,000）。
输出是一个浮点数组。
两个数组都针对 sse 进行了正确的内存对齐。
以下代码仅使用 1 个 xmm 寄存器，在 linux 上运行，并且可以通过 gcc -O3 code.cpp -lrt -msse4.2
编译
谢谢。
#include <iostream>
#include <emmintrin.h>
#include <time.h>
using namespace std;
void print_xmm(__m128 xmm){
float out[4];
_mm_storeu_ps(out,xmm);
int i;
for (i = 0; i < 4; ++i) std::cout << out[i] << " ";
std::cout << std::endl;
}
void print_arr(float* ptr, size_t size){
size_t i;
for(i = 0; i < size; ++i){
cout << ptr[i] << " ";
}
cout << endl;
}
int main(void){
size_t size = 25000 * 4;
// this has to be multiple of 4
size_t repeat = 10000;
// test 10000 cycles of the code
float* ar_in = (float*)aligned_alloc(16, size*sizeof(float));
float* ar_out = (float*)aligned_alloc(16, size*sizeof(float));
//fill test data into the input array
//the data is an array of positive numbers.
size_t i;
for (i = 0; i < size; ++i){
ar_in[i] = (i+1) * (i+1);
}
//prepare for recipical square root.
__m128 xmm0;
size_t size_fix = size*sizeof(float)/sizeof(__m128);
float* ar_in_end = ar_in + size_fix;
float* ar_out_now;
float* ar_in_now;
//timing
struct timespec tp_start, tp_end;
i = repeat;
clock_gettime(CLOCK_MONOTONIC, &tp_start);
//start timing
while(--i){
ar_out_now = ar_out;
for(ar_in_now = ar_in;
ar_in_now != ar_in_end;
ar_in_now += 4, ar_out_now+=4){
//4 = sizeof(__m128)/sizeof(float);
xmm0 = _mm_load_ps(ar_in_now);
//cout << "load xmm: ";
//print_xmm(xmm0);
xmm0 = _mm_rsqrt_ps(xmm0);
//cout << "rsqrt xmm: ";
//print_xmm(xmm0);
_mm_store_ps(ar_out_now,xmm0);
}
}
//end timing
clock_gettime(CLOCK_MONOTONIC, &tp_end);
double timing;
const double nano = 0.000000001;
timing = ((double)(tp_end.tv_sec - tp_start.tv_sec )
+ (tp_end.tv_nsec - tp_start.tv_nsec) * nano)/repeat;
cout << " timing per cycle: " << timing << endl;
/*
cout << "input array: ";
print_arr(ar_in, size);
cout << "output array: ";
print_arr(ar_out,size);
*/
//free mem
free(ar_in);
free(ar_out);
return 0;
}
当然，你必须测量它，但是有已知代码可以计算（不是很精确）平方根倒数，检查 https://www.beyond3d.com/content/articles/8/
float InvSqrt (float x) {
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0x5f3759df - (i>>1);
x = *(float*)&i;
x = x*(1.5f - xhalf*x*x);
return x;
}
已测试（使用 VS2015 和 GCC 5.4.0）转换为 SSE2，link
__m128 InvSqrtSSE2(__m128 x) {
__m128 xhalf = _mm_mul_ps(x, _mm_set1_ps(0.5f));
x = _mm_castsi128_ps(_mm_sub_epi32(_mm_set1_epi32(0x5f3759df), _mm_srai_epi32(_mm_castps_si128(x), 1)));
return _mm_mul_ps(x, _mm_sub_ps(_mm_set1_ps(1.5f), _mm_mul_ps(xhalf, _mm_mul_ps(x, x))));
}
更新我
我认错了！感谢@EOF，他注意到不正确的转换，这些被替换为强制转换
因为这是一个算术强度很低的流式计算，所以你几乎可以肯定是内存受限的。如果您使用非时间加载和存储，您可能会发现可测量的加速。
// Hint to the CPU that we don't want to use the cache
// Be sure to update this if you use manual loop unrolling
_mm_prefetch(reinterpret_cast<char*>(ar_in+4), _MM_HINT_NTA);
// Hint to the CPU that we don't need to write back through the cache
_mm_stream_ps(ar_out_now,xmm0);
编辑:
我 运行 进行了一些测试，以了解不同硬件上的情况。不出所料，结果对所使用的硬件非常敏感。我认为这可能是由于现代 CPU 中 read/write 缓冲区数量的增加。
所有代码都是使用 gcc-6.1 编译的
gcc -std=c++14 -O3 -march=native -mavx -mfpmath=sse -ffast-math
英特尔酷睿 i3-3120M @ 2.5GHz； 3MB 缓存
OP's code: 350 +- 10 milliseconds
NTA Prefetch: 360 +- 5 milliseconds
NTA Prefetch+NTA store: 430 +- 10 milliseconds
英特尔酷睿 i7-6500U CPU @ 2.50GHz； 4MB 缓存
OP's code: 205 +- 5 milliseconds
NTA Prefetch: 220 +- 2 milliseconds
NTA Prefetch+NTA store: 200 +- 5 milliseconds
将问题大小增加到 2MB，与 OP 的解决方案相比，NTA Prefetch+NTA 存储的运行时间减少了约 30%。
结果：问题规模太小，无法从 NTA 中获益。在较旧的架构上，这是有害的。在较新的架构上，它与 OP 的解决方案不相上下。
结论：在这种情况下可能不值得付出额外的努力。
你的浮点数组有多大？如果它在 L1（或可能是 L2）中已经很热，则此代码的 gcc5.3 输出会成为现代 Intel CPU 上 uop 吞吐量的瓶颈，因为它使用 6 个融合域 uops 进行循环，每次迭代执行一个向量。 （因此它将 运行 每 2 个周期一个向量）。
要在现代 Intel CPU 上实现每个时钟 1 个向量的吞吐量，您需要展开循环（请参阅下文了解为什么非展开的 asm 无法工作）。让编译器为您完成它可能很好（而不是在 C++ 源代码中手动完成）。例如使用配置文件引导优化 (gcc -fprofile-use)，或者只是盲目地使用 -funroll-loops.
理论上每个时钟 16 个字节足以使一个内核的主内存带宽饱和。然而，IIRC Z Boson 观察到使用多核的带宽更好，这可能是因为多核保持更多未完成的请求，并且一个核上的停顿不会使内存空闲。但是，如果输入在内核的 L2 中很热，则最好使用该内核处理数据。
在 Haswell 或更高版本上，每个时钟加载和存储 16 个字节只是 L1 缓存带宽的一半，因此您需要一个 AVX 版本来最大化每个核心带宽。
如果你遇到内存瓶颈，you might want to do a Newton-Raphson iteration to get a nearly-full accuracy 1/sqrt(x)，特别是如果你对一个大数组使用多个线程。（因为如果单个线程不能每个时钟维持一个加载+存储。）
或者稍后加载此数据时可能只是即时使用 rsqrt。它非常便宜，具有高吞吐量，但延迟仍然类似于 FP 添加。同样，如果它是一个不适合缓存的大数组，通过减少对数据的单独传递来增加计算强度是一件大事。 （Cache Blocking aka Loop Tiling 也是一个好主意，如果您可以这样做：运行 您的算法在缓存大小的数据块上执行多个步骤。）
只有在找不到有效利用缓存的方法时，才使用绕过缓存的 NT 存储作为最后的手段。如果你能转换一些你将要使用的数据就更好了，这样下次使用时它就在缓存中。
主循环（不幸的是来自 .L31 to jne .L31 on the Godbolt compiler explorer) is 6 uops for Intel SnB-family CPUs, because indexed addressing modes don't micro-fuse. (This isn't yet documented in Agner Fog's microarch pdf。）
Nehalem 上有 4 个融合域微指令，只有三个 ALU 微指令，所以 Nehalem 应该 运行 每个时钟 1 个。
.L31: # the main loop: 6 uops on SnB-family, 4 uops on Nehalem
rsqrtps xmm0, XMMWORD PTR [rbx+rax] # tmp127, MEM[base: ar_in_now_10, index: ivtmp.51_61, offset: 0B]
movaps XMMWORD PTR [rbp+0+rax], xmm0 # MEM[base: ar_out_now_12, index: ivtmp.51_61, offset: 0B], tmp127
add rax, 16 # ivtmp.51,
cmp rax, 100000 # ivtmp.51,
jne .L31 #,
由于您想编写一个单独的目的地，因此无法将循环减少到 4 个融合域微指令，因此它可以 运行 每个时钟一个向量而不展开。 （加载和存储都需要是单寄存器寻址模式，因此使用由 current_dst 索引的 src - dst 而不是递增 src 的技巧不起作用）。
修改您的 C++，使 gcc 使用指针增量只会节省一个 uop，因为您必须增加 src 和 dst。即 float *endp = start + length; 和 for (p = start ; p < endp ; p+=4) {} 会形成一个像
这样的循环
.loop:
rsqrtps xmm0, [rsi]
add rsi, 16
movaps [rdi], xmm0
add rdi, 16
cmp rdi, rbx
jne .loop
希望 gcc 在展开时会做类似的事情，否则如果 rsqrtps + movaps 仍然使用索引寻址模式，它们将自己成为 4 个融合域微指令，再多的展开也不会让你的循环 运行 在每个时钟一个向量。
如何使用 popcnt 和 SSE4.2 在 cpu 上更快地计算数组的近似倒数平方根？
输入是存储在浮点数组中的正整数（范围从 0 到大约 200,000）。
输出是一个浮点数组。
两个数组都针对 sse 进行了正确的内存对齐。
以下代码仅使用 1 个 xmm 寄存器，在 linux 上运行，并且可以通过 gcc -O3 code.cpp -lrt -msse4.2
谢谢。
#include <iostream>
#include <emmintrin.h>
#include <time.h>
using namespace std;
void print_xmm(__m128 xmm){
float out[4];
_mm_storeu_ps(out,xmm);
int i;
for (i = 0; i < 4; ++i) std::cout << out[i] << " ";
std::cout << std::endl;
}
void print_arr(float* ptr, size_t size){
size_t i;
for(i = 0; i < size; ++i){
cout << ptr[i] << " ";
}
cout << endl;
}
int main(void){
size_t size = 25000 * 4;
// this has to be multiple of 4
size_t repeat = 10000;
// test 10000 cycles of the code
float* ar_in = (float*)aligned_alloc(16, size*sizeof(float));
float* ar_out = (float*)aligned_alloc(16, size*sizeof(float));
//fill test data into the input array
//the data is an array of positive numbers.
size_t i;
for (i = 0; i < size; ++i){
ar_in[i] = (i+1) * (i+1);
}
//prepare for recipical square root.
__m128 xmm0;
size_t size_fix = size*sizeof(float)/sizeof(__m128);
float* ar_in_end = ar_in + size_fix;
float* ar_out_now;
float* ar_in_now;
//timing
struct timespec tp_start, tp_end;
i = repeat;
clock_gettime(CLOCK_MONOTONIC, &tp_start);
//start timing
while(--i){
ar_out_now = ar_out;
for(ar_in_now = ar_in;
ar_in_now != ar_in_end;
ar_in_now += 4, ar_out_now+=4){
//4 = sizeof(__m128)/sizeof(float);
xmm0 = _mm_load_ps(ar_in_now);
//cout << "load xmm: ";
//print_xmm(xmm0);
xmm0 = _mm_rsqrt_ps(xmm0);
//cout << "rsqrt xmm: ";
//print_xmm(xmm0);
_mm_store_ps(ar_out_now,xmm0);
}
}
//end timing
clock_gettime(CLOCK_MONOTONIC, &tp_end);
double timing;
const double nano = 0.000000001;
timing = ((double)(tp_end.tv_sec - tp_start.tv_sec )
+ (tp_end.tv_nsec - tp_start.tv_nsec) * nano)/repeat;
cout << " timing per cycle: " << timing << endl;
/*
cout << "input array: ";
print_arr(ar_in, size);
cout << "output array: ";
print_arr(ar_out,size);
*/
//free mem
free(ar_in);
free(ar_out);
return 0;
}
当然，你必须测量它，但是有已知代码可以计算（不是很精确）平方根倒数，检查 https://www.beyond3d.com/content/articles/8/
float InvSqrt (float x) {
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0x5f3759df - (i>>1);
x = *(float*)&i;
x = x*(1.5f - xhalf*x*x);
return x;
}
已测试（使用 VS2015 和 GCC 5.4.0）转换为 SSE2，link
__m128 InvSqrtSSE2(__m128 x) {
__m128 xhalf = _mm_mul_ps(x, _mm_set1_ps(0.5f));
x = _mm_castsi128_ps(_mm_sub_epi32(_mm_set1_epi32(0x5f3759df), _mm_srai_epi32(_mm_castps_si128(x), 1)));
return _mm_mul_ps(x, _mm_sub_ps(_mm_set1_ps(1.5f), _mm_mul_ps(xhalf, _mm_mul_ps(x, x))));
}
更新我
我认错了！感谢@EOF，他注意到不正确的转换，这些被替换为强制转换
因为这是一个算术强度很低的流式计算，所以你几乎可以肯定是内存受限的。如果您使用非时间加载和存储，您可能会发现可测量的加速。
// Hint to the CPU that we don't want to use the cache
// Be sure to update this if you use manual loop unrolling
_mm_prefetch(reinterpret_cast<char*>(ar_in+4), _MM_HINT_NTA);
// Hint to the CPU that we don't need to write back through the cache
_mm_stream_ps(ar_out_now,xmm0);
编辑:
我 运行 进行了一些测试，以了解不同硬件上的情况。不出所料，结果对所使用的硬件非常敏感。我认为这可能是由于现代 CPU 中 read/write 缓冲区数量的增加。
所有代码都是使用 gcc-6.1 编译的
gcc -std=c++14 -O3 -march=native -mavx -mfpmath=sse -ffast-math
英特尔酷睿 i3-3120M @ 2.5GHz； 3MB 缓存
OP's code: 350 +- 10 milliseconds
NTA Prefetch: 360 +- 5 milliseconds
NTA Prefetch+NTA store: 430 +- 10 milliseconds
英特尔酷睿 i7-6500U CPU @ 2.50GHz； 4MB 缓存
OP's code: 205 +- 5 milliseconds
NTA Prefetch: 220 +- 2 milliseconds
NTA Prefetch+NTA store: 200 +- 5 milliseconds
将问题大小增加到 2MB，与 OP 的解决方案相比，NTA Prefetch+NTA 存储的运行时间减少了约 30%。
结果：问题规模太小，无法从 NTA 中获益。在较旧的架构上，这是有害的。在较新的架构上，它与 OP 的解决方案不相上下。
结论：在这种情况下可能不值得付出额外的努力。
你的浮点数组有多大？如果它在 L1（或可能是 L2）中已经很热，则此代码的 gcc5.3 输出会成为现代 Intel CPU 上 uop 吞吐量的瓶颈，因为它使用 6 个融合域 uops 进行循环，每次迭代执行一个向量。 （因此它将 运行 每 2 个周期一个向量）。
要在现代 Intel CPU 上实现每个时钟 1 个向量的吞吐量，您需要展开循环（请参阅下文了解为什么非展开的 asm 无法工作）。让编译器为您完成它可能很好（而不是在 C++ 源代码中手动完成）。例如使用配置文件引导优化 (gcc -fprofile-use)，或者只是盲目地使用 -funroll-loops.
理论上每个时钟 16 个字节足以使一个内核的主内存带宽饱和。然而，IIRC Z Boson 观察到使用多核的带宽更好，这可能是因为多核保持更多未完成的请求，并且一个核上的停顿不会使内存空闲。但是，如果输入在内核的 L2 中很热，则最好使用该内核处理数据。
在 Haswell 或更高版本上，每个时钟加载和存储 16 个字节只是 L1 缓存带宽的一半，因此您需要一个 AVX 版本来最大化每个核心带宽。
如果你遇到内存瓶颈，you might want to do a Newton-Raphson iteration to get a nearly-full accuracy 1/sqrt(x)，特别是如果你对一个大数组使用多个线程。（因为如果单个线程不能每个时钟维持一个加载+存储。）
或者稍后加载此数据时可能只是即时使用 rsqrt。它非常便宜，具有高吞吐量，但延迟仍然类似于 FP 添加。同样，如果它是一个不适合缓存的大数组，通过减少对数据的单独传递来增加计算强度是一件大事。 （Cache Blocking aka Loop Tiling 也是一个好主意，如果您可以这样做：运行 您的算法在缓存大小的数据块上执行多个步骤。）
只有在找不到有效利用缓存的方法时，才使用绕过缓存的 NT 存储作为最后的手段。如果你能转换一些你将要使用的数据就更好了，这样下次使用时它就在缓存中。
主循环（不幸的是来自 .L31 to jne .L31 on the Godbolt compiler explorer) is 6 uops for Intel SnB-family CPUs, because indexed addressing modes don't micro-fuse. (This isn't yet documented in Agner Fog's microarch pdf。）
Nehalem 上有 4 个融合域微指令，只有三个 ALU 微指令，所以 Nehalem 应该 运行 每个时钟 1 个。
.L31: # the main loop: 6 uops on SnB-family, 4 uops on Nehalem
rsqrtps xmm0, XMMWORD PTR [rbx+rax] # tmp127, MEM[base: ar_in_now_10, index: ivtmp.51_61, offset: 0B]
movaps XMMWORD PTR [rbp+0+rax], xmm0 # MEM[base: ar_out_now_12, index: ivtmp.51_61, offset: 0B], tmp127
add rax, 16 # ivtmp.51,
cmp rax, 100000 # ivtmp.51,
jne .L31 #,
由于您想编写一个单独的目的地，因此无法将循环减少到 4 个融合域微指令，因此它可以 运行 每个时钟一个向量而不展开。 （加载和存储都需要是单寄存器寻址模式，因此使用由 current_dst 索引的 src - dst 而不是递增 src 的技巧不起作用）。
修改您的 C++，使 gcc 使用指针增量只会节省一个 uop，因为您必须增加 src 和 dst。即 float *endp = start + length; 和 for (p = start ; p < endp ; p+=4) {} 会形成一个像
.loop:
rsqrtps xmm0, [rsi]
add rsi, 16
movaps [rdi], xmm0
add rdi, 16
cmp rdi, rbx
jne .loop
希望 gcc 在展开时会做类似的事情，否则如果 rsqrtps + movaps 仍然使用索引寻址模式，它们将自己成为 4 个融合域微指令，再多的展开也不会让你的循环 运行 在每个时钟一个向量。