在 C++ 中使用内在函数检查 nans
checking for nans with intrinsics in c++
我刚开始使用内在函数,但我想编写一个函数,该函数采用 4 个双精度计算的向量a > 1e-5 ? std::sqrt(a) : 0.0
我的第一直觉是按如下方式编写
#include <immintrin.h>
__m256d f(__m256d a)
{
__m256d is_valid = a > _mm256_set1_pd(1e-5);
__m256d sqrt_val = _mm256_sqrt_pd(a);
return is_valid * sqrt_val;
}
根据 gcc.godbolt.com 编译为以下内容
f(double __vector(4)):
vsqrtpd ymm1, ymm0
vcmpgtpd ymm0, ymm0, YMMWORD PTR .LC0[rip]
vmulpd ymm0, ymm1, ymm0
ret
.LC0:
.long 2296604913
.long 1055193269
.long 2296604913
.long 1055193269
.long 2296604913
.long 1055193269
.long 2296604913
.long 1055193269
但我担心如果 sqrt_val
包含 nan
会发生什么。我认为 0.0 * nan
行不通。此处的最佳做法是什么?
编辑
阅读@ChrisCooper(和@njuffa)的评论后,我被链接到另一个堆栈溢出答案,所以我将测试自我平等,然后 and
这与我的结果。
#include <immintrin.h>
__m256d f(__m256d a)
{
__m256d is_valid = a > _mm256_set1_pd(1e-5);
__m256d sqrt_val = _mm256_sqrt_pd(a);
__m256d result = is_valid * sqrt_val;
__m256d cmpeq = result == result;
return _mm256_and_pd(cmpeq, result);
}
编译为以下内容
f(double __vector(4)):
vsqrtpd ymm1, ymm0
vcmpgtpd ymm0, ymm0, YMMWORD PTR .LC0[rip]
vmulpd ymm0, ymm1, ymm0
vcmpeqpd ymm1, ymm0, ymm0
vandpd ymm0, ymm1, ymm0
ret
.LC0:
.long 2296604913
.long 1055193269
.long 2296604913
.long 1055193269
.long 2296604913
.long 1055193269
.long 2296604913
.long 1055193269
我之前没有使用 AVX 内在函数进行编程,因此从文档中收集信息以快速将下面的代码组合在一起。对于我提供的一个测试用例,它似乎可以正常工作。
相关的观察是比较指令return全1(如果结果为TRUE)或全0(如果结果为FALSE)的掩码。然后可以使用此掩码通过将掩码与 vsqrtpd
的结果进行与操作来有条件地将平方根的结果设置为零。 0.0
在IEEE-754双精度中的二进制表示全为0。
之前没有使用过这些内部函数,我发现比较谓词使用起来很棘手。据我了解,这里我们想使用有序比较来获得关于 NaN 的所需行为(也就是说,与 NaN 的比较应该导致 FALSE),因此 'O' 变体。我们也不希望 NaN 输入触发异常(也就是说,我们希望在这种情况下比较安静),因此 'Q' 变体。这意味着我们要使用 _CMP_GT_OQ.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <immintrin.h>
__m256d f (__m256d a)
{
double em5 = 1e-5;
__m256d v_em5 = _mm256_broadcast_sd (&em5);
__m256d v_sqrt = _mm256_sqrt_pd (a);
__m256d v_mask = _mm256_cmp_pd (a, v_em5, _CMP_GT_OQ);
__m256d v_res = _mm256_and_pd (v_sqrt, v_mask);
return v_res;
}
int main (void)
{
__m256d arg, res;
double args[4] = {2e-5, sqrt(-1.0), 1e-6, -1.0};
double ress [4] = {0};
memcpy (&arg, args, sizeof(arg));
res = f (arg);
memcpy (ress, &res, sizeof(res));
printf ("args = % 23.16e % 23.16e % 23.16e % 23.16e\n",
args[0], args[1], args[2], args[3]);
printf ("ress = % 23.16e % 23.16e % 23.16e % 23.16e\n",
ress[0], ress[1], ress[2], ress[3]);
return EXIT_SUCCESS;
}
我用 Intel C 编译器编译了上面的程序,输出结果是这样的:
args = 2.0000000000000002e-005 -1.#IND000000000000e+000 9.9999999999999995e-007 -1.0000000000000000e+000
ress = 4.4721359549995798e-003 0.0000000000000000e+000 0.0000000000000000e+000 0.0000000000000000e+000
这里,1.#IND000000000000e+000
是一个名为 INDEFINITE 的特定 QNaN。
我刚开始使用内在函数,但我想编写一个函数,该函数采用 4 个双精度计算的向量a > 1e-5 ? std::sqrt(a) : 0.0
我的第一直觉是按如下方式编写
#include <immintrin.h>
__m256d f(__m256d a)
{
__m256d is_valid = a > _mm256_set1_pd(1e-5);
__m256d sqrt_val = _mm256_sqrt_pd(a);
return is_valid * sqrt_val;
}
根据 gcc.godbolt.com 编译为以下内容
f(double __vector(4)):
vsqrtpd ymm1, ymm0
vcmpgtpd ymm0, ymm0, YMMWORD PTR .LC0[rip]
vmulpd ymm0, ymm1, ymm0
ret
.LC0:
.long 2296604913
.long 1055193269
.long 2296604913
.long 1055193269
.long 2296604913
.long 1055193269
.long 2296604913
.long 1055193269
但我担心如果 sqrt_val
包含 nan
会发生什么。我认为 0.0 * nan
行不通。此处的最佳做法是什么?
编辑
阅读@ChrisCooper(和@njuffa)的评论后,我被链接到另一个堆栈溢出答案,所以我将测试自我平等,然后 and
这与我的结果。
#include <immintrin.h>
__m256d f(__m256d a)
{
__m256d is_valid = a > _mm256_set1_pd(1e-5);
__m256d sqrt_val = _mm256_sqrt_pd(a);
__m256d result = is_valid * sqrt_val;
__m256d cmpeq = result == result;
return _mm256_and_pd(cmpeq, result);
}
编译为以下内容
f(double __vector(4)):
vsqrtpd ymm1, ymm0
vcmpgtpd ymm0, ymm0, YMMWORD PTR .LC0[rip]
vmulpd ymm0, ymm1, ymm0
vcmpeqpd ymm1, ymm0, ymm0
vandpd ymm0, ymm1, ymm0
ret
.LC0:
.long 2296604913
.long 1055193269
.long 2296604913
.long 1055193269
.long 2296604913
.long 1055193269
.long 2296604913
.long 1055193269
我之前没有使用 AVX 内在函数进行编程,因此从文档中收集信息以快速将下面的代码组合在一起。对于我提供的一个测试用例,它似乎可以正常工作。
相关的观察是比较指令return全1(如果结果为TRUE)或全0(如果结果为FALSE)的掩码。然后可以使用此掩码通过将掩码与 vsqrtpd
的结果进行与操作来有条件地将平方根的结果设置为零。 0.0
在IEEE-754双精度中的二进制表示全为0。
之前没有使用过这些内部函数,我发现比较谓词使用起来很棘手。据我了解,这里我们想使用有序比较来获得关于 NaN 的所需行为(也就是说,与 NaN 的比较应该导致 FALSE),因此 'O' 变体。我们也不希望 NaN 输入触发异常(也就是说,我们希望在这种情况下比较安静),因此 'Q' 变体。这意味着我们要使用 _CMP_GT_OQ.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <immintrin.h>
__m256d f (__m256d a)
{
double em5 = 1e-5;
__m256d v_em5 = _mm256_broadcast_sd (&em5);
__m256d v_sqrt = _mm256_sqrt_pd (a);
__m256d v_mask = _mm256_cmp_pd (a, v_em5, _CMP_GT_OQ);
__m256d v_res = _mm256_and_pd (v_sqrt, v_mask);
return v_res;
}
int main (void)
{
__m256d arg, res;
double args[4] = {2e-5, sqrt(-1.0), 1e-6, -1.0};
double ress [4] = {0};
memcpy (&arg, args, sizeof(arg));
res = f (arg);
memcpy (ress, &res, sizeof(res));
printf ("args = % 23.16e % 23.16e % 23.16e % 23.16e\n",
args[0], args[1], args[2], args[3]);
printf ("ress = % 23.16e % 23.16e % 23.16e % 23.16e\n",
ress[0], ress[1], ress[2], ress[3]);
return EXIT_SUCCESS;
}
我用 Intel C 编译器编译了上面的程序,输出结果是这样的:
args = 2.0000000000000002e-005 -1.#IND000000000000e+000 9.9999999999999995e-007 -1.0000000000000000e+000
ress = 4.4721359549995798e-003 0.0000000000000000e+000 0.0000000000000000e+000 0.0000000000000000e+000
这里,1.#IND000000000000e+000
是一个名为 INDEFINITE 的特定 QNaN。