GPS数据矢量化JAVA_算法 – acosf()的精确矢量化实现

问题中代码的无分支版本是可能的(几乎没有任何冗余工作,只有一些比较/混合来为FMA创建常量),但IDK如果编译器将自动矢量化它.

如果所有元素都有 – | a |,则主要的额外工作是无用的sqrt / fma > -0.5625f,不幸的是在关键路径上.

asinf_core的arg是(r> -0.5625f)？ r：sqrtf(fmaf(0.5f,r,0.5f)).

与此同时,您(或编译器)可以在输出上混合FMA的系数.

如果你通过将pi / 2常量放入一个浮点而不是用2个常数被乘数来创建它来牺牲pi / 2常量的精度,你可以

fmaf( condition?-1:2, asinf_core_result, condition ? pi/2 : 0)

因此,您可以在两个常量之间进行选择,或者使用SIMD比较结果和一个常量来有条件地将其归零(例如x86 SSE).

最终修正基于原始输入的范围检查,因此FP混合和asinf_core的FMA工作之间再次存在指令级并行性.

事实上,我们可以通过在第二个条件下将常量输入混合到asinf_core的输出来优化它到之前的FMA.我们希望asinf_core作为它的被乘数之一,所以我们可以通过否定常数来否定. (SIMD实现可能执行a_cmp = andnot(a>> 0.0f,a> = – 1.0f),然后乘数^( – 0.0f& a_cmp),其中乘数在之前有条件地完成.

输出上该FMA的附加常数为0,pi / 2,pi或pi pi / 2.给定两个比较结果(对于非NaN情况,在a和on r = – | a |上),我们可以将其组合成2位整数并将其用作随机控制以从向量中选择FP常量所有4个常数,例如使用AVX vpermilps(带有可变控制的快速通道内洗牌).即不是混合4种不同的方式,而是使用shuffle作为2位LUT！

如果我们这样做,我们也应该为乘法常数做,因为创建常数是主要成本.可变混合比x86上的shuffle更昂贵(通常为2 uops对1).在Skylake上,可变混合(如vblendvps)可以使用任何端口(而shuffle只能在端口5上运行).有足够的ILP,这可能会影响整体uop吞吐量或整体ALU端口,而不是端口5.(Haswell上的可变混合对于端口5是2 uops,所以它严格比vpermilps ymm,ymm,ymm更差).

我们将从-1,1,-2和2中进行选择.

带有三元运算符的标量,使用gcc7.3 -O3 -march = skylake -ffast-math自动矢量化(使用8个vblendvps).自动向量化所需的快速数学：/不幸的是,gcc仍然使用rsqrtps进行牛顿迭代(没有FMA？！？),即使是-mrecip=none, which I thought was supposed to disable this也是如此.

使用clang5.0(具有相同选项)仅使用5个vblendvps进行自动向量化.请参见on the Godbolt compiler explorer.这会编译,看起来可能是正确数量的指令,但在其他方面未经测试.

// I think this is far more than enough digits for float precision, but wouldn't hurt to use a standard constant instead of what I typed from memory.

static const float pi_2 = 3.1415926535897932384626433 / 2;

static const float pi = 3.1415926535897932384626433;

//static const float pi_plus_pi_2 = 3.1415926535897932384626433 * 3.0 / 2;

/* maximum error UNKNOWN, completely UNTESTED */

float my_acosf_branchless (float a)

{

float r = (a > 0.0f) ? (-a) : a; // avoid modifying the "sign" of NaNs

bool a_in_range = !(a > 0.0f) && (a >= -1.0f);

bool rsmall = (r > -0.5625f);

float asinf_arg = rsmall ? r : sqrtf (fmaf (0.5f, r, 0.5f));

float asinf_res = asinf_core(asinf_arg);

#if 0

r = fmaf( rsmall?-1.0f:2.0f, asinf_res, rsmall ? pi_2 : 0);

if (!(a > 0.0f) && (a >= -1.0f)) { // avoid modifying the "sign" of NaNs

/* arccos (-x) = pi - arccos(x) */

r = fmaf (0x1.ddcb02p+0f, 0x1.aee9d6p+0f, -r);

}

#else

float fma_mul = rsmall? -1.0f:2.0f;

fma_mul = a_in_range ? -fma_mul : fma_mul;

float fma_add = rsmall ? pi_2 : 0;

fma_add = a_in_range ? fma_add + pi : fma_add;

// to vectorize, turn the 2 conditions into a 2-bit integer.

// Use vpermilps as a 2-bit LUT of float constants

// clang doesn't see the LUT trick, but otherwise appears non-terrible at this blending.

r = fmaf(asinf_res, fma_mul, fma_add);

#endif

return r;

}

使用循环测试自动矢量化,该循环在1024个对齐的浮点元素的数组上运行它;看看Godbolt的链接.

TODO：内在版本.