准确的矢量化acosf（）实现

Question

acosf()如果平台支持融合乘加（FMA），则相对于无限精确（数学）结果，简单的实现可以轻松实现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;
}

Answer 1

问题中的代码的无分支版本是可能的（几乎没有任何冗余工作，只是一些比较/混合来为 FMA 创建常量），但不知道编译器是否会自动向量化它。

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

的参数asinf_core是(r > -0.5625f) ? r : sqrtf (fmaf (0.5f, r, 0.5f)).

与此同时，您（或编译器）可以在输出上混合 FMA 的系数。

pi/2如果您通过将常数放入 1float而不是用 2 个常数被乘数创建它来牺牲常数的准确性fmaf，则可以

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

因此，您可以在两个常量或andps具有 SIMD 比较结果的常量之间进行选择，以有条件地将其归零（例如 x86 SSE）。

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

事实上，我们可以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|非 NaN 情况，我们可以将其组合成一个 2 位整数，并将其用作洗牌控制，从所有 4 个常量的向量中选择一个 FP 常量，例如使用AVX vpermilps（带可变控制的快速车道内洗牌）。即，不要混合 4 种不同的方式，而是使用 shuffle 作为 2 位 LUT！

如果我们这样做，我们也应该对乘法常数这样做，因为创建常数是主要成本。变量混合比 x86 上的 shuffle 更昂贵（通常是 2 uops 对 1 uops）。在 Skylake 上，变量混合（如vblendvps）可以使用任何端口（而随机播放仅在端口 5 上运行）。有足够的 ILP，这可能会成为整体 uop 吞吐量或整体 ALU 端口（而不是端口 5）的瓶颈。（Haswell 上的变量混合对于端口 5 来说是 2 uop，因此它比 更糟糕vpermilps ymm,ymm,ymm）。

我们将从 -1、1、-2 和 2 中进行选择。

具有三元运算符 的标量，使用 进行自动向量化（使用 8 vblendvps）gcc7.3 -O3 -march=skylake -ffast-math。自动向量化所需的快速数学：/不幸的是，gcc仍然使用rsqrtps+牛顿迭代（没有FMA？！？），即使有-mrecip=none，我认为应该禁用这个。

使用 clang5.0仅使用 5 个自动矢量化vblendvps（具有相同的选项）。请在 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 链接。

TODO：内在函数版本。