如果平台支持融合乘加(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)的新闻帖子中的一个想法，他在新闻帖子中观察到，对于[0,1]中的x，acos(x)≈ √(2*(1-x))，并且多项式可以提供必要的比例因子，使其在整个区间内成为准确的近似值。从下面的代码中可以看出，这种方法产生直线代码，只需两次使用三元运算符来处理负半平面中的参数。
#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 acosn",
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.13an",
a[i], b[i], ref, dref);
printf ("test FAILEDn");
return EXIT_FAILURE;
}
}
printf ("^^^^ argi = %08xn", argi);
} while (argi);
printf ("test PASSEDn");
free (a);
free (b);
return EXIT_SUCCESS;
}
尽管此代码的结构似乎有利于自动矢量化，但在使用 Compiler Explorer 提供的编译器针对AVX2时，我运气不佳。唯一似乎能够在上面测试应用程序的内部循环上下文中对此代码进行矢量化的编译器是 Clang。但是 Clang 似乎只有在我指定-ffast-math的情况下才能做到这一点，但是，这具有将sqrtf()调用转换为通过rsqrt计算的近似平方根的不良副作用。我尝试了一些侵入性较小的开关，例如-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;
}
问题中代码的无分支版本是可能的(几乎没有任何冗余工作，只是一些比较/混合为 FMA 创建常量)，但如果编译器会自动矢量化它，则 IDK 会自动矢量化。
如果所有元素都-|a| > -0.5625f，主要的额外工作是无用的sqrt/fma，不幸的是在关键路径上。
asinf_core的参数是(r > -0.5625f) ? r : sqrtf (fmaf (0.5f, r, 0.5f)).
与此同时，您(或编译器)可以在输出上混合 FMA 的系数。
如果您通过将pi/2常数放入一个float而不是用 2 个常量乘法来创建它来牺牲fmaf的准确性，则可以
fmaf( condition?-1:2, asinf_core_result, condition ? pi/2 : 0)
因此，您可以在两个常量之间进行选择，或者andps一个带有 SIMD 比较结果的常量，以有条件地将其归零(例如 x86 SSE)。
最终的修正基于原始输入的范围检查，因此FP混合和asinf_core的FMA工作之间再次存在指令级并行性。
事实上，我们可以将其优化到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。 鉴于两个比较结果(对于非 NaN 情况的a和r=-|a|)，我们也许可以将其组合成一个 2 位整数，并将其用作随机控制，以从所有 4 个常量的向量中选择一个 FP 常量，例如使用 AVXvpermilps(带有变量控制的快速车道内洗牌)。 即，与其混合 4 种不同的方式，不如使用随机播放作为 2 位 LUT！
如果我们这样做，我们也应该为乘法常数这样做，因为创建常量是主要成本。 可变混合比 x86 上的随机播放更昂贵(通常为 2 uops 对 1)。 在 Skylake 上，可变混合(如vblendvps)可以使用任何端口(而随机播放仅在端口 5 上运行)。 有足够的 ILP，这可能会成为整体 uop 吞吐量或整体 ALU 端口的瓶颈，而不是端口 5。 (对于端口 5，Haswell 上的可变混合为 2 uops，因此严格来说比vpermilps ymm,ymm,ymm差)。
我们将从 -1、1、-2 和 2 中进行选择。
带有三元运算符的标量，使用gcc7.3 -O3 -march=skylake -ffast-math自动矢量化(使用 8vblendvps)。 自动矢量化所需的快速数学：/不幸的是，gcc 仍然使用rsqrtps+ 牛顿迭代(没有 FMA?!?)，即使使用-mrecip=none，我认为应该禁用它。
使用 clang5.0 只需 5vblendvps即可自动矢量化(使用相同的选项)。 在 Godbolt 编译器资源管理器上查看两者。 这编译了，看起来可能是适量的指令，但未经测试。
// 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 个对齐的float元素数组上运行;请参阅 Godbolt 链接。
待办事项：内部版本。
这不是一种替代算法方法，但尽管如此 您可能会对这句扩展评论感兴趣。
似乎使用 gcc，函数copysignf()比 三元运算符。在下面的代码中，我重写了你的标量 带copysignf()的解决方案 而不是三元运算符。
即使使用相当旧的 gcc 4.9 编译器，代码也会矢量化，使用 选项gcc -std=c99 -O3 -m64 -Wall -march=haswell -fno-math-errno.sqrtf()函数矢量化为vsqrtps指令。Godbolt链接在这里。
#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 https://stackoverflow.com/a/49091530/2439725 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;
}