在 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。