用c ++中的内在函数检查nans

Question

我是新手使用内在函数但是我想编写一个函数,它接受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

Answer 1

我之前没有使用AVX内在函数进行编程,因此从文档中收集信息以快速汇总下面的代码.它似乎对我提供的一个测试用例有所帮助.

相关的观察是比较指令返回所有1的掩码(如果结果为TRUE)或全0(如果结果为FALSE).然后,可以使用此掩码通过使用结果来对掩码进行AND运算来有条理地将平方根的结果设置为零vsqrtpd.0.0IEEE-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.