acosf() 的精确矢量化实现
Accurate vectorizable implementation of acosf()
如果平台支持融合乘加 (FMA)，acosf() 的简单实现可以轻松实现 1.5 ulp 的误差界限（相对于无限精确的（数学）结果）。这意味着在舍入到最近或偶数模式下，结果与正确舍入的结果相差不会超过一个 ulp。
然而，这样的实现通常包括两个主要代码分支，它们将主要近似区间 [0,1] 大致分成两半，如下面的示例代码所示。当以 SIMD 架构为目标时，这种分支性会抑制编译器的自动矢量化。
是否有一种替代算法方法可以更轻松地进行自动矢量化，同时保持 1.5 ulps 的相同误差范围？可以假定平台支持 FMA。
/* approximate arcsin(a) on [-0.5625,+0.5625], max ulp err = 0.95080 */
float asinf_core(float a)
{
float r, s;
s = a * a;
r = 0x1.a7f260p-5f; // 5.17513156e-2
r = fmaf (r, s, 0x1.29a5cep-6f); // 1.81669723e-2
r = fmaf (r, s, 0x1.7f0842p-5f); // 4.67568673e-2
r = fmaf (r, s, 0x1.329256p-4f); // 7.48465881e-2
r = fmaf (r, s, 0x1.555728p-3f); // 1.66670144e-1
r = r * s;
r = fmaf (r, a, a);
return r;
}
/* maximum error = 1.45667 ulp */
float my_acosf (float a)
{
float r;
r = (a > 0.0f) ? (-a) : a; // avoid modifying the "sign" of NaNs
if (r > -0.5625f) {
/* arccos(x) = pi/2 - arcsin(x) */
r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, asinf_core (r));
} else {
/* arccos(x) = 2 * arcsin (sqrt ((1-x) / 2)) */
r = 2.0f * asinf_core (sqrtf (fmaf (0.5f, r, 0.5f)));
}
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);
}
return r;
}
我最接近满意的解决方案是基于 Robert Harley news posting 的想法，他在其中观察到对于 [0,1] 中的 x，acos(x) ≈ √ (2*(1-x))，多项式可以提供必要的比例因子，使其成为整个区间的精确近似值。从下面的代码中可以看出，这种方法导致 straight-line 代码仅两次使用三元运算符来处理负参数 half-plane.
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#define VECTORIZABLE 1
#define ARR_LEN (1 << 24)
#define MAX_ULP 1 /* deviation from correctly rounded result */
#if VECTORIZABLE
/*
Compute arccos(a) with a maximum error of 1.496766 ulp
This uses an idea from Robert Harley's posting in comp.arch.arithmetic on 1996/07/12
https://groups.google.com/forum/#!original/comp.arch.arithmetic/wqCPkCCXqWs/T9qCkHtGE2YJ
*/
float my_acosf (float a)
{
float r, s, t;
s = (a < 0.0f) ? 2.0f : (-2.0f);
t = fmaf (s, a, 2.0f);
s = sqrtf (t);
r = 0x1.c86000p-22f; // 4.25032340e-7
r = fmaf (r, t, -0x1.0258fap-19f); // -1.92483935e-6
r = fmaf (r, t, 0x1.90c5c4p-18f); // 5.97197595e-6
r = fmaf (r, t, -0x1.55668cp-19f); // -2.54363249e-6
r = fmaf (r, t, 0x1.c3f78ap-16f); // 2.69393295e-5
r = fmaf (r, t, 0x1.e8f446p-14f); // 1.16575764e-4
r = fmaf (r, t, 0x1.6df072p-11f); // 6.97973708e-4
r = fmaf (r, t, 0x1.3332a6p-08f); // 4.68746712e-3
r = fmaf (r, t, 0x1.555550p-05f); // 4.16666567e-2
r = r * t;
r = fmaf (r, s, s);
t = fmaf (0x1.ddcb02p+0f, 0x1.aee9d6p+0f, 0.0f - r); // PI-r
r = (a < 0.0f) ? t : r;
return r;
}
#else // VECTORIZABLE
/* approximate arcsin(a) on [-0.5625,+0.5625], max ulp err = 0.95080 */
float asinf_core(float a)
{
float r, s;
s = a * a;
r = 0x1.a7f260p-5f; // 5.17513156e-2
r = fmaf (r, s, 0x1.29a5cep-6f); // 1.81669723e-2
r = fmaf (r, s, 0x1.7f0842p-5f); // 4.67568673e-2
r = fmaf (r, s, 0x1.329256p-4f); // 7.48465881e-2
r = fmaf (r, s, 0x1.555728p-3f); // 1.66670144e-1
r = r * s;
r = fmaf (r, a, a);
return r;
}
/* maximum error = 1.45667 ulp */
float my_acosf (float a)
{
float r;
r = (a > 0.0f) ? (-a) : a; // avoid modifying the "sign" of NaNs
if (r > -0.5625f) {
/* arccos(x) = pi/2 - arcsin(x) */
r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, asinf_core (r));
} else {
/* arccos(x) = 2 * arcsin (sqrt ((1-x) / 2)) */
r = 2.0f * asinf_core (sqrtf (fmaf (0.5f, r, 0.5f)));
}
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);
}
return r;
}
#endif // VECTORIZABLE
int main (void)
{
double darg, dref;
float ref, *a, *b;
uint32_t argi, resi, refi;
printf ("%svectorizable implementation of acos\n",
VECTORIZABLE ? "" : "non-");
a = (float *)malloc (sizeof(a[0]) * ARR_LEN);
b = (float *)malloc (sizeof(b[0]) * ARR_LEN);
argi = 0x00000000;
do {
for (int i = 0; i < ARR_LEN; i++) {
memcpy (&a[i], &argi, sizeof(a[i]));
argi++;
}
for (int i = 0; i < ARR_LEN; i++) {
b[i] = my_acosf (a[i]);
}
for (int i = 0; i < ARR_LEN; i++) {
darg = (double)a[i];
dref = acos (darg);
ref = (float)dref;
memcpy (&refi, &ref, sizeof(refi));
memcpy (&resi, &b[i], sizeof(resi));
if (llabs ((long long int)resi - (long long int)refi) > MAX_ULP) {
printf ("error > 1 ulp a[i]=% 14.6a b[i]=% 14.6a ref=% 14.6a dref=% 21.13a\n",
a[i], b[i], ref, dref);
printf ("test FAILED\n");
return EXIT_FAILURE;
}
}
printf ("^^^^ argi = %08x\n", argi);
} while (argi);
printf ("test PASSED\n");
free (a);
free (b);
return EXIT_SUCCESS;
}
虽然这段代码的结构看起来有利于 auto-vectorization，但我在使用 Compiler Explorer 提供的编译器以 AVX2 为目标时运气并不好。唯一似乎能够在上面我的测试应用程序的内部循环上下文中对该代码进行矢量化的编译器是 Clang。但是 Clang 似乎只有在我指定 -ffast-math 时才能做到这一点，但是，它具有将 sqrtf() 调用转换为通过 [=15 计算的近似平方根的不希望的 side-effect =].我尝试了一些侵入性较小的开关，例如-fno-honor-nans、-fno-math-errno、-fno-trapping-math，但是 my_acosf() 即使我组合使用它们也没有矢量化。
目前，我已将上述代码手动翻译为 AVX2 + FMA 内在函数，如下所示：
#include "immintrin.h"
/* maximum error = 1.496766 ulp */
__m256 _mm256_acos_ps (__m256 x)
{
const __m256 zero= _mm256_set1_ps ( 0.0f);
const __m256 two = _mm256_set1_ps ( 2.0f);
const __m256 mtwo= _mm256_set1_ps (-2.0f);
const __m256 c0 = _mm256_set1_ps ( 0x1.c86000p-22f); // 4.25032340e-7
const __m256 c1 = _mm256_set1_ps (-0x1.0258fap-19f); // -1.92483935e-6
const __m256 c2 = _mm256_set1_ps ( 0x1.90c5c4p-18f); // 5.97197595e-6
const __m256 c3 = _mm256_set1_ps (-0x1.55668cp-19f); // -2.54363249e-6
const __m256 c4 = _mm256_set1_ps ( 0x1.c3f78ap-16f); // 2.69393295e-5
const __m256 c5 = _mm256_set1_ps ( 0x1.e8f446p-14f); // 1.16575764e-4
const __m256 c6 = _mm256_set1_ps ( 0x1.6df072p-11f); // 6.97973708e-4
const __m256 c7 = _mm256_set1_ps ( 0x1.3332a6p-8f); // 4.68746712e-3
const __m256 c8 = _mm256_set1_ps ( 0x1.555550p-5f); // 4.16666567e-2
const __m256 pi0 = _mm256_set1_ps ( 0x1.ddcb02p+0f); // 1.86637890e+0
const __m256 pi1 = _mm256_set1_ps ( 0x1.aee9d6p+0f); // 1.68325555e+0
__m256 s, r, t, m;
s = two;
t = mtwo;
m = _mm256_cmp_ps (x, zero, _CMP_LT_OQ);
t = _mm256_blendv_ps (t, s, m);
t = _mm256_fmadd_ps (x, t, s);
s = _mm256_sqrt_ps (t);
r = c0;
r = _mm256_fmadd_ps (r, t, c1);
r = _mm256_fmadd_ps (r, t, c2);
r = _mm256_fmadd_ps (r, t, c3);
r = _mm256_fmadd_ps (r, t, c4);
r = _mm256_fmadd_ps (r, t, c5);
r = _mm256_fmadd_ps (r, t, c6);
r = _mm256_fmadd_ps (r, t, c7);
r = _mm256_fmadd_ps (r, t, c8);
r = _mm256_mul_ps (r, t);
r = _mm256_fmadd_ps (r, s, s);
t = _mm256_sub_ps (zero, r);
t = _mm256_fmadd_ps (pi0, pi1, t);
r = _mm256_blendv_ps (r, t, m);
return r;
}
问题中代码的无分支版本是可能的（几乎没有任何多余的工作，只是一些 compare/blends 为 FMA 创建常量），但如果编译器会 auto-vectorize 的话，IDK。
主要的额外工作是无用的 sqrt / fma 如果所有元素都有 -|a| > -0.5625f，不幸的是在关键路径上。
asinf_core 的参数是 (r > -0.5625f) ? r : sqrtf (fmaf (0.5f, r, 0.5f))。
与此同时，您（或编译器）可以在输出中混合 FMA 的系数。
如果您通过将 pi/2 常量放入一个 float 而不是用 2 个常量乘数创建它以 fmaf 来牺牲 pi/2 常量的准确性，您可以
fmaf( condition?-1:2, asinf_core_result, condition ? pi/2 : 0)
所以你 select 在两个常量之间，或者 andps 一个常量与 SIMD 比较结果有条件地将它归零（例如 x86 SSE）。
最终修复基于原始输入的 range-check，因此 FP 混合与 asinf_core 的 FMA 工作之间再次存在 instruction-level 并行性。
事实上，我们可以通过在第二个条件下将常量输入混合到 asinf_core 的输出上，将其优化到之前的 FMA 中。我们要asinf_core作为它的被乘数之一，所以我们可以通过对常数取反来取反或不取反。 （SIMD 实现可能会 a_cmp = andnot( a>0.0f, a>=-1.0f)，然后 multiplier ^ (-0.0f & a_cmp)，其中 multiplier 之前是有条件地完成的。
该 FMA 在输出上的附加常数是 0、pi/2、pi 或 pi + pi/2。给定两个比较结果（在 a 和 r=-|a| 上 non-NaN 情况下），我们也许可以将其组合成一个 2 位整数并将其用作 [= 的混洗控制91=] 来自所有 4 个常量的向量的 FP 常量，例如使用 AVX vpermilps（使用变量控制快速 in-lane 随机播放）。即 不是混合 4 种不同的方式，而是使用 shuffle 作为 2 位 LUT!
如果我们这样做，我们也应该为乘法常数这样做，因为创建常数是主要成本。可变混合比 x86 上的随机播放更昂贵（通常是 2 微指令对 1）。在 Skylake 上，变量混合（如 vblendvps）可以使用任何端口（而随机播放仅在端口 5 上 运行）。有足够的 ILP，这可能是整个 uop 吞吐量或整个 ALU 端口的瓶颈，而不是端口 5。（Haswell 上的变量混合是端口 5 的 2 uops，所以它严格来说比 vpermilps ymm,ymm,ymm 差）。
我们将从 -1、1、-2 和 2 开始 select。
具有三元运算符的标量、auto-vectorizes（具有 8 个 vblendvps）和 gcc7.3 -O3 -march=skylake -ffast-math。 fast-math 需要自动矢量化：/ 不幸的是，gcc 仍然使用 rsqrtps + 牛顿迭代（没有 FMA？！？），即使 .
使用 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;
}
Auto-vectorization 使用一个循环进行测试，该循环 运行 将其遍历 1024 个对齐的 float 元素的数组；看到神箭 link.
TODO：内在函数版本。
这不是完全替代的算法方法，但是尽管如此
您可能会对这条扩展评论感兴趣。
看来，使用 gcc，函数 copysignf() 比
三元运算符。在下面的代码中，我重写了你的标量
与 copysignf()
而不是三元运算符。
即使使用相当老的 gcc 4.9 编译器，代码也可以矢量化，
选项 gcc -std=c99 -O3 -m64 -Wall -march=haswell -fno-math-errno。
sqrtf() 函数被矢量化为 vsqrtps 指令。 Godbolt link is here.
#include <stdio.h>
#include <immintrin.h>
#include <math.h>
float acosf_cpsgn (float a)
{
float r, s, t, pi2;
/* s = (a < 0.0f) ? 2.0f : (-2.0f); */
s = copysignf(2.0f, -a);
t = fmaf (s, a, 2.0f);
s = sqrtf (t);
r = 0x1.c86000p-22f; // 4.25032340e-7
r = fmaf (r, t, -0x1.0258fap-19f); // -1.92483935e-6
r = fmaf (r, t, 0x1.90c5c4p-18f); // 5.97197595e-6
r = fmaf (r, t, -0x1.55668cp-19f); // -2.54363249e-6
r = fmaf (r, t, 0x1.c3f78ap-16f); // 2.69393295e-5
r = fmaf (r, t, 0x1.e8f446p-14f); // 1.16575764e-4
r = fmaf (r, t, 0x1.6df072p-11f); // 6.97973708e-4
r = fmaf (r, t, 0x1.3332a6p-08f); // 4.68746712e-3
r = fmaf (r, t, 0x1.555550p-05f); // 4.16666567e-2
r = r * t;
r = fmaf (r, s, s);
/* t = fmaf (0x1.ddcb02p+0f, 0x1.aee9d6p+0f, 0.0f - r); // PI-r */
/* r = (a < 0.0f) ? t : r; */
r = copysignf(r, a);
pi2 = 0x1.ddcb02p+0f * 0.5f; /* no rounding here */
pi2 = pi2 - copysignf(pi2, a); /* no rounding here */
t = fmaf (pi2, 0x1.aee9d6p+0f, r); // PI-r
return t;
}
float my_acosf (float a)
{
float r, s, t;
s = (a < 0.0f) ? 2.0f : (-2.0f);
t = fmaf (s, a, 2.0f);
s = sqrtf (t);
r = 0x1.c86000p-22f; // 4.25032340e-7
r = fmaf (r, t, -0x1.0258fap-19f); // -1.92483935e-6
r = fmaf (r, t, 0x1.90c5c4p-18f); // 5.97197595e-6
r = fmaf (r, t, -0x1.55668cp-19f); // -2.54363249e-6
r = fmaf (r, t, 0x1.c3f78ap-16f); // 2.69393295e-5
r = fmaf (r, t, 0x1.e8f446p-14f); // 1.16575764e-4
r = fmaf (r, t, 0x1.6df072p-11f); // 6.97973708e-4
r = fmaf (r, t, 0x1.3332a6p-08f); // 4.68746712e-3
r = fmaf (r, t, 0x1.555550p-05f); // 4.16666567e-2
r = r * t;
r = fmaf (r, s, s);
t = fmaf (0x1.ddcb02p+0f, 0x1.aee9d6p+0f, 0.0f - r); // PI-r
r = (a < 0.0f) ? t : r;
return r;
}
/* The code from the next 2 functions is copied from the godbold link in Peter cordes' */
/* answer and modified */
int autovec_test_a (float *__restrict dst, float *__restrict src) {
dst = __builtin_assume_aligned(dst,32);
src = __builtin_assume_aligned(src,32);
for (int i=0 ; i<1024 ; i++ ) {
dst[i] = my_acosf(src[i]);
}
return 0;
}
int autovec_test_b (float *__restrict dst, float *__restrict src) {
dst = __builtin_assume_aligned(dst,32);
src = __builtin_assume_aligned(src,32);
for (int i=0 ; i<1024 ; i++ ) {
dst[i] = acosf_cpsgn(src[i]);
}
return 0;
}
如果平台支持融合乘加 (FMA)，acosf() 的简单实现可以轻松实现 1.5 ulp 的误差界限（相对于无限精确的（数学）结果）。这意味着在舍入到最近或偶数模式下，结果与正确舍入的结果相差不会超过一个 ulp。
然而，这样的实现通常包括两个主要代码分支，它们将主要近似区间 [0,1] 大致分成两半，如下面的示例代码所示。当以 SIMD 架构为目标时，这种分支性会抑制编译器的自动矢量化。
是否有一种替代算法方法可以更轻松地进行自动矢量化，同时保持 1.5 ulps 的相同误差范围？可以假定平台支持 FMA。
/* approximate arcsin(a) on [-0.5625,+0.5625], max ulp err = 0.95080 */
float asinf_core(float a)
{
float r, s;
s = a * a;
r = 0x1.a7f260p-5f; // 5.17513156e-2
r = fmaf (r, s, 0x1.29a5cep-6f); // 1.81669723e-2
r = fmaf (r, s, 0x1.7f0842p-5f); // 4.67568673e-2
r = fmaf (r, s, 0x1.329256p-4f); // 7.48465881e-2
r = fmaf (r, s, 0x1.555728p-3f); // 1.66670144e-1
r = r * s;
r = fmaf (r, a, a);
return r;
}
/* maximum error = 1.45667 ulp */
float my_acosf (float a)
{
float r;
r = (a > 0.0f) ? (-a) : a; // avoid modifying the "sign" of NaNs
if (r > -0.5625f) {
/* arccos(x) = pi/2 - arcsin(x) */
r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, asinf_core (r));
} else {
/* arccos(x) = 2 * arcsin (sqrt ((1-x) / 2)) */
r = 2.0f * asinf_core (sqrtf (fmaf (0.5f, r, 0.5f)));
}
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);
}
return r;
}
我最接近满意的解决方案是基于 Robert Harley news posting 的想法，他在其中观察到对于 [0,1] 中的 x，acos(x) ≈ √ (2*(1-x))，多项式可以提供必要的比例因子，使其成为整个区间的精确近似值。从下面的代码中可以看出，这种方法导致 straight-line 代码仅两次使用三元运算符来处理负参数 half-plane.
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#define VECTORIZABLE 1
#define ARR_LEN (1 << 24)
#define MAX_ULP 1 /* deviation from correctly rounded result */
#if VECTORIZABLE
/*
Compute arccos(a) with a maximum error of 1.496766 ulp
This uses an idea from Robert Harley's posting in comp.arch.arithmetic on 1996/07/12
https://groups.google.com/forum/#!original/comp.arch.arithmetic/wqCPkCCXqWs/T9qCkHtGE2YJ
*/
float my_acosf (float a)
{
float r, s, t;
s = (a < 0.0f) ? 2.0f : (-2.0f);
t = fmaf (s, a, 2.0f);
s = sqrtf (t);
r = 0x1.c86000p-22f; // 4.25032340e-7
r = fmaf (r, t, -0x1.0258fap-19f); // -1.92483935e-6
r = fmaf (r, t, 0x1.90c5c4p-18f); // 5.97197595e-6
r = fmaf (r, t, -0x1.55668cp-19f); // -2.54363249e-6
r = fmaf (r, t, 0x1.c3f78ap-16f); // 2.69393295e-5
r = fmaf (r, t, 0x1.e8f446p-14f); // 1.16575764e-4
r = fmaf (r, t, 0x1.6df072p-11f); // 6.97973708e-4
r = fmaf (r, t, 0x1.3332a6p-08f); // 4.68746712e-3
r = fmaf (r, t, 0x1.555550p-05f); // 4.16666567e-2
r = r * t;
r = fmaf (r, s, s);
t = fmaf (0x1.ddcb02p+0f, 0x1.aee9d6p+0f, 0.0f - r); // PI-r
r = (a < 0.0f) ? t : r;
return r;
}
#else // VECTORIZABLE
/* approximate arcsin(a) on [-0.5625,+0.5625], max ulp err = 0.95080 */
float asinf_core(float a)
{
float r, s;
s = a * a;
r = 0x1.a7f260p-5f; // 5.17513156e-2
r = fmaf (r, s, 0x1.29a5cep-6f); // 1.81669723e-2
r = fmaf (r, s, 0x1.7f0842p-5f); // 4.67568673e-2
r = fmaf (r, s, 0x1.329256p-4f); // 7.48465881e-2
r = fmaf (r, s, 0x1.555728p-3f); // 1.66670144e-1
r = r * s;
r = fmaf (r, a, a);
return r;
}
/* maximum error = 1.45667 ulp */
float my_acosf (float a)
{
float r;
r = (a > 0.0f) ? (-a) : a; // avoid modifying the "sign" of NaNs
if (r > -0.5625f) {
/* arccos(x) = pi/2 - arcsin(x) */
r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, asinf_core (r));
} else {
/* arccos(x) = 2 * arcsin (sqrt ((1-x) / 2)) */
r = 2.0f * asinf_core (sqrtf (fmaf (0.5f, r, 0.5f)));
}
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);
}
return r;
}
#endif // VECTORIZABLE
int main (void)
{
double darg, dref;
float ref, *a, *b;
uint32_t argi, resi, refi;
printf ("%svectorizable implementation of acos\n",
VECTORIZABLE ? "" : "non-");
a = (float *)malloc (sizeof(a[0]) * ARR_LEN);
b = (float *)malloc (sizeof(b[0]) * ARR_LEN);
argi = 0x00000000;
do {
for (int i = 0; i < ARR_LEN; i++) {
memcpy (&a[i], &argi, sizeof(a[i]));
argi++;
}
for (int i = 0; i < ARR_LEN; i++) {
b[i] = my_acosf (a[i]);
}
for (int i = 0; i < ARR_LEN; i++) {
darg = (double)a[i];
dref = acos (darg);
ref = (float)dref;
memcpy (&refi, &ref, sizeof(refi));
memcpy (&resi, &b[i], sizeof(resi));
if (llabs ((long long int)resi - (long long int)refi) > MAX_ULP) {
printf ("error > 1 ulp a[i]=% 14.6a b[i]=% 14.6a ref=% 14.6a dref=% 21.13a\n",
a[i], b[i], ref, dref);
printf ("test FAILED\n");
return EXIT_FAILURE;
}
}
printf ("^^^^ argi = %08x\n", argi);
} while (argi);
printf ("test PASSED\n");
free (a);
free (b);
return EXIT_SUCCESS;
}
虽然这段代码的结构看起来有利于 auto-vectorization，但我在使用 Compiler Explorer 提供的编译器以 AVX2 为目标时运气并不好。唯一似乎能够在上面我的测试应用程序的内部循环上下文中对该代码进行矢量化的编译器是 Clang。但是 Clang 似乎只有在我指定 -ffast-math 时才能做到这一点，但是，它具有将 sqrtf() 调用转换为通过 [=15 计算的近似平方根的不希望的 side-effect =].我尝试了一些侵入性较小的开关，例如-fno-honor-nans、-fno-math-errno、-fno-trapping-math，但是 my_acosf() 即使我组合使用它们也没有矢量化。
目前，我已将上述代码手动翻译为 AVX2 + FMA 内在函数，如下所示：
#include "immintrin.h"
/* maximum error = 1.496766 ulp */
__m256 _mm256_acos_ps (__m256 x)
{
const __m256 zero= _mm256_set1_ps ( 0.0f);
const __m256 two = _mm256_set1_ps ( 2.0f);
const __m256 mtwo= _mm256_set1_ps (-2.0f);
const __m256 c0 = _mm256_set1_ps ( 0x1.c86000p-22f); // 4.25032340e-7
const __m256 c1 = _mm256_set1_ps (-0x1.0258fap-19f); // -1.92483935e-6
const __m256 c2 = _mm256_set1_ps ( 0x1.90c5c4p-18f); // 5.97197595e-6
const __m256 c3 = _mm256_set1_ps (-0x1.55668cp-19f); // -2.54363249e-6
const __m256 c4 = _mm256_set1_ps ( 0x1.c3f78ap-16f); // 2.69393295e-5
const __m256 c5 = _mm256_set1_ps ( 0x1.e8f446p-14f); // 1.16575764e-4
const __m256 c6 = _mm256_set1_ps ( 0x1.6df072p-11f); // 6.97973708e-4
const __m256 c7 = _mm256_set1_ps ( 0x1.3332a6p-8f); // 4.68746712e-3
const __m256 c8 = _mm256_set1_ps ( 0x1.555550p-5f); // 4.16666567e-2
const __m256 pi0 = _mm256_set1_ps ( 0x1.ddcb02p+0f); // 1.86637890e+0
const __m256 pi1 = _mm256_set1_ps ( 0x1.aee9d6p+0f); // 1.68325555e+0
__m256 s, r, t, m;
s = two;
t = mtwo;
m = _mm256_cmp_ps (x, zero, _CMP_LT_OQ);
t = _mm256_blendv_ps (t, s, m);
t = _mm256_fmadd_ps (x, t, s);
s = _mm256_sqrt_ps (t);
r = c0;
r = _mm256_fmadd_ps (r, t, c1);
r = _mm256_fmadd_ps (r, t, c2);
r = _mm256_fmadd_ps (r, t, c3);
r = _mm256_fmadd_ps (r, t, c4);
r = _mm256_fmadd_ps (r, t, c5);
r = _mm256_fmadd_ps (r, t, c6);
r = _mm256_fmadd_ps (r, t, c7);
r = _mm256_fmadd_ps (r, t, c8);
r = _mm256_mul_ps (r, t);
r = _mm256_fmadd_ps (r, s, s);
t = _mm256_sub_ps (zero, r);
t = _mm256_fmadd_ps (pi0, pi1, t);
r = _mm256_blendv_ps (r, t, m);
return r;
}
问题中代码的无分支版本是可能的（几乎没有任何多余的工作，只是一些 compare/blends 为 FMA 创建常量），但如果编译器会 auto-vectorize 的话，IDK。
主要的额外工作是无用的 sqrt / fma 如果所有元素都有 -|a| > -0.5625f，不幸的是在关键路径上。
asinf_core 的参数是 (r > -0.5625f) ? r : sqrtf (fmaf (0.5f, r, 0.5f))。
与此同时，您（或编译器）可以在输出中混合 FMA 的系数。
如果您通过将 pi/2 常量放入一个 float 而不是用 2 个常量乘数创建它以 fmaf 来牺牲 pi/2 常量的准确性，您可以
fmaf( condition?-1:2, asinf_core_result, condition ? pi/2 : 0)
所以你 select 在两个常量之间，或者 andps 一个常量与 SIMD 比较结果有条件地将它归零（例如 x86 SSE）。
最终修复基于原始输入的 range-check，因此 FP 混合与 asinf_core 的 FMA 工作之间再次存在 instruction-level 并行性。
事实上，我们可以通过在第二个条件下将常量输入混合到 asinf_core 的输出上，将其优化到之前的 FMA 中。我们要asinf_core作为它的被乘数之一，所以我们可以通过对常数取反来取反或不取反。 （SIMD 实现可能会 a_cmp = andnot( a>0.0f, a>=-1.0f)，然后 multiplier ^ (-0.0f & a_cmp)，其中 multiplier 之前是有条件地完成的。
该 FMA 在输出上的附加常数是 0、pi/2、pi 或 pi + pi/2。给定两个比较结果（在 a 和 r=-|a| 上 non-NaN 情况下），我们也许可以将其组合成一个 2 位整数并将其用作 [= 的混洗控制91=] 来自所有 4 个常量的向量的 FP 常量，例如使用 AVX vpermilps（使用变量控制快速 in-lane 随机播放）。即 不是混合 4 种不同的方式，而是使用 shuffle 作为 2 位 LUT!
如果我们这样做，我们也应该为乘法常数这样做，因为创建常数是主要成本。可变混合比 x86 上的随机播放更昂贵（通常是 2 微指令对 1）。在 Skylake 上，变量混合（如 vblendvps）可以使用任何端口（而随机播放仅在端口 5 上 运行）。有足够的 ILP，这可能是整个 uop 吞吐量或整个 ALU 端口的瓶颈，而不是端口 5。（Haswell 上的变量混合是端口 5 的 2 uops，所以它严格来说比 vpermilps ymm,ymm,ymm 差）。
我们将从 -1、1、-2 和 2 开始 select。
具有三元运算符的标量、auto-vectorizes（具有 8 个 vblendvps）和 gcc7.3 -O3 -march=skylake -ffast-math。 fast-math 需要自动矢量化：/ 不幸的是，gcc 仍然使用 rsqrtps + 牛顿迭代（没有 FMA？！？），即使
使用 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;
}
Auto-vectorization 使用一个循环进行测试，该循环 运行 将其遍历 1024 个对齐的 float 元素的数组；看到神箭 link.
TODO：内在函数版本。
这不是完全替代的算法方法，但是尽管如此 您可能会对这条扩展评论感兴趣。
看来，使用 gcc，函数 copysignf() 比
三元运算符。在下面的代码中，我重写了你的标量
copysignf()
而不是三元运算符。
即使使用相当老的 gcc 4.9 编译器，代码也可以矢量化，
选项 gcc -std=c99 -O3 -m64 -Wall -march=haswell -fno-math-errno。
sqrtf() 函数被矢量化为 vsqrtps 指令。 Godbolt link is here.
#include <stdio.h>
#include <immintrin.h>
#include <math.h>
float acosf_cpsgn (float a)
{
float r, s, t, pi2;
/* s = (a < 0.0f) ? 2.0f : (-2.0f); */
s = copysignf(2.0f, -a);
t = fmaf (s, a, 2.0f);
s = sqrtf (t);
r = 0x1.c86000p-22f; // 4.25032340e-7
r = fmaf (r, t, -0x1.0258fap-19f); // -1.92483935e-6
r = fmaf (r, t, 0x1.90c5c4p-18f); // 5.97197595e-6
r = fmaf (r, t, -0x1.55668cp-19f); // -2.54363249e-6
r = fmaf (r, t, 0x1.c3f78ap-16f); // 2.69393295e-5
r = fmaf (r, t, 0x1.e8f446p-14f); // 1.16575764e-4
r = fmaf (r, t, 0x1.6df072p-11f); // 6.97973708e-4
r = fmaf (r, t, 0x1.3332a6p-08f); // 4.68746712e-3
r = fmaf (r, t, 0x1.555550p-05f); // 4.16666567e-2
r = r * t;
r = fmaf (r, s, s);
/* t = fmaf (0x1.ddcb02p+0f, 0x1.aee9d6p+0f, 0.0f - r); // PI-r */
/* r = (a < 0.0f) ? t : r; */
r = copysignf(r, a);
pi2 = 0x1.ddcb02p+0f * 0.5f; /* no rounding here */
pi2 = pi2 - copysignf(pi2, a); /* no rounding here */
t = fmaf (pi2, 0x1.aee9d6p+0f, r); // PI-r
return t;
}
float my_acosf (float a)
{
float r, s, t;
s = (a < 0.0f) ? 2.0f : (-2.0f);
t = fmaf (s, a, 2.0f);
s = sqrtf (t);
r = 0x1.c86000p-22f; // 4.25032340e-7
r = fmaf (r, t, -0x1.0258fap-19f); // -1.92483935e-6
r = fmaf (r, t, 0x1.90c5c4p-18f); // 5.97197595e-6
r = fmaf (r, t, -0x1.55668cp-19f); // -2.54363249e-6
r = fmaf (r, t, 0x1.c3f78ap-16f); // 2.69393295e-5
r = fmaf (r, t, 0x1.e8f446p-14f); // 1.16575764e-4
r = fmaf (r, t, 0x1.6df072p-11f); // 6.97973708e-4
r = fmaf (r, t, 0x1.3332a6p-08f); // 4.68746712e-3
r = fmaf (r, t, 0x1.555550p-05f); // 4.16666567e-2
r = r * t;
r = fmaf (r, s, s);
t = fmaf (0x1.ddcb02p+0f, 0x1.aee9d6p+0f, 0.0f - r); // PI-r
r = (a < 0.0f) ? t : r;
return r;
}
/* The code from the next 2 functions is copied from the godbold link in Peter cordes' */
/* answer and modified */
int autovec_test_a (float *__restrict dst, float *__restrict src) {
dst = __builtin_assume_aligned(dst,32);
src = __builtin_assume_aligned(src,32);
for (int i=0 ; i<1024 ; i++ ) {
dst[i] = my_acosf(src[i]);
}
return 0;
}
int autovec_test_b (float *__restrict dst, float *__restrict src) {
dst = __builtin_assume_aligned(dst,32);
src = __builtin_assume_aligned(src,32);
for (int i=0 ; i<1024 ; i++ ) {
dst[i] = acosf_cpsgn(src[i]);
}
return 0;
}