在 x86（使用 SSE2）和 ARM（使用 vfpv4 NEON）上尾数为 11 位的 atan2 近似值

Posted 2023-02-16

【中文标题】在 x86（使用 SSE2）和 ARM（使用 vfpv4 NEON）上尾数为 11 位的 atan2 近似值

【英文标题】：atan2 approximation with 11bits in mantissa on x86(with SSE2) and ARM(with vfpv4 NEON)

【发布时间】：2017-09-14 04:46:32

【问题描述】：

我正在尝试在尾数中实现 11 位精度的快速 atan2(float)。 atan2 实现将用于图像处理。 所以最好用 SIMD 指令来实现（针对 x86（带 SSE2）和 ARM（带 vpfv4 NEON）的 impl）。

目前，我使用切比雪夫多项式逼近(https://jp.mathworks.com/help/fixedpoint/examples/calculate-fixed-point-arctangent.html)。

我愿意手动实现矢量化代码。 我将使用 SSE2（或更高版本）和 NEON 包装库。 我计划或尝试了这些实施选项

矢量化多项式逼近 切比雪夫多项式（现已实现） 标量查找表（+ 线性插值） 矢量化查找表（这可能吗？我对 NEON 不熟悉，所以我不知道 NEON 中的 VGATHER 之类的一些指令） 矢量化 CORDIC 方法（这可能会很慢，因为它需要 12 次以上的旋转迭代才能获得 11 位精度）

还有什么好主意？

【问题讨论】：

通常你会通过制作一个将向量作为输入的版本来向量化一个数学函数。这通常比尝试使用 SIMD 来加快对单个输入的评估更有效。

@PeterCordes 通常，您可以获得 gcc 或 clang 的自动向量化器，为您提供编写良好的近似值的向量化。由于缺少 ieee754 一致性（缺少非规范数字），只有 ARM NEON 将其搞砸了。

@EOF 这真的有效吗？标量实现通常会有点复杂，或者使用表查找。 LUT 可能仍然是一个胜利，但 gcc 不会自动矢量化解包/查找/重新打包。例如看这个 AVX2 __m256d log2 函数***.com/a/45898937/224132，它使用了一个收集指令并且只有一个小的多项式。

如果您基本上要手动矢量化，直接编写矢量化代码可能更容易，而不是自动矢量化器转换的标量代码。虽然您可能必须为 NEON 和 SSE2 内在函数编写两次，或者使用通用包装库或 GNU C“本机”向量语法，但您可能会发现 NEON 与 SSE2 的最佳策略不同。

@PeterCordes 我将使用 SSE2（或更高版本）和 NEON 包装库。我计划或尝试了这些实现选项 * 多项式近似 * chebyshev 多项式（现已实现） * 标量查找表（+ 线性插值） * 矢量化查找表（这可能吗？我不熟悉 NEON，所以我不知道一些指令，如 AVX2 中的 VGATHER）* CORDIC 方法（这可能很慢，因为它需要 12 次以上的旋转迭代才能获得 11 位精度）

【参考方案1】：

首先要检查的是，当您的编译器应用于floats 数组时，您的编译器是否能够向量化atan2f (y,x)。这通常需要至少一个高优化级别，例如-O3，并且可能指定一个宽松的“快速数学”模式，其中errno 处理、非正规和特殊输入（例如无穷大和NaN）在很大程度上被忽略。使用这种方法，准确度可能会远远超出要求，但在性能方面可能很难超越经过仔细调整的库实现。

接下来要尝试的是编写一个具有足够精度的简单标量实现，并让编译器对其进行向量化。通常，这意味着除了可以通过 if-conversion 转换为无分支代码的非常简单的分支之外，避免其他任何事情。此类代码的一个示例是如下所示的fast_atan2f()。使用 Intel 编译器，调用为 icl /O3 /fp:precise /Qvec_report=2 my_atan2f.c，这已成功矢量化：my_atan2f.c(67): (col. 9) remark: LOOP WAS VECTORIZED. 通过反汇编仔细检查生成的代码表明，fast_atan2f() 已使用 *ps 风格的 SSE 指令进行内联和矢量化。

如果一切都失败了，您可以手动将fast_atan2() 的代码转换为特定于平台的 SIMD 内部函数，这应该不难做到。不过，我没有足够的经验来快速完成。

我在这段代码中使用了一个非常简单的算法，它是一个简单的参数归约，以在 [0,1] 中产生一个归约参数，然后是一个极小极大多项式近似，最后一步将结果映射回完整的圆圈。核心近似是使用 Remez 算法生成的，并使用二阶霍纳方案进行评估。在可用的情况下，可以使用融合乘法 (FMA)。出于性能考虑，不处理涉及无穷大、NaN 或 0/0 的特殊情况。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* maximum relative error about 3.6e-5 */
float fast_atan2f (float y, float x)

    float a, r, s, t, c, q, ax, ay, mx, mn;
    ax = fabsf (x);
    ay = fabsf (y);
    mx = fmaxf (ay, ax);
    mn = fminf (ay, ax);
    a = mn / mx;
    /* Minimax polynomial approximation to atan(a) on [0,1] */
    s = a * a;
    c = s * a;
    q = s * s;
    r =  0.024840285f * q + 0.18681418f;
    t = -0.094097948f * q - 0.33213072f;
    r = r * s + t;
    r = r * c + a;
    /* Map to full circle */
    if (ay > ax) r = 1.57079637f - r;
    if (x <   0) r = 3.14159274f - r;
    if (y <   0) r = -r;
    return r;


/* Fixes via: Greg Rose, KISS: A Bit Too Simple. http://eprint.iacr.org/2011/007 */
static unsigned int z=362436069,w=521288629,jsr=362436069,jcong=123456789;
#define znew (z=36969*(z&0xffff)+(z>>16))
#define wnew (w=18000*(w&0xffff)+(w>>16))
#define MWC  ((znew<<16)+wnew)
#define SHR3 (jsr^=(jsr<<13),jsr^=(jsr>>17),jsr^=(jsr<<5)) /* 2^32-1 */
#define CONG (jcong=69069*jcong+13579)                     /* 2^32 */
#define KISS ((MWC^CONG)+SHR3)

float rand_float(void)

    volatile union 
        float f;
        unsigned int i;
     cvt;
    do 
        cvt.i = KISS;
     while (isnan(cvt.f) || isinf (cvt.f) || (fabsf (cvt.f) < powf (2.0f, -126)));
    return cvt.f;


int main (void)

    const int N = 10000;
    const int M = 100000;
    float ref, relerr, maxrelerr = 0.0f;
    float arga[N], argb[N], res[N];
    int i, j;

    printf ("testing atan2() with %d test vectors\n", N*M);

    for (j = 0; j < M; j++) 
        for (i = 0; i < N; i++) 
            arga[i] = rand_float();
            argb[i] = rand_float();
        

        // This loop should be vectorized
        for (i = 0; i < N; i++) 
            res[i] = fast_atan2f (argb[i], arga[i]);
        

        for (i = 0; i < N; i++) 
            ref = (float) atan2 ((double)argb[i], (double)arga[i]);
            relerr = (res[i] - ref) / ref;
            if ((fabsf (relerr) > maxrelerr) && 
                (fabsf (ref) >= powf (2.0f, -126)))  // result not denormal
                maxrelerr = fabsf (relerr);
            
        
    ;

    printf ("max rel err = % 15.8e\n\n", maxrelerr);

    printf ("atan2(1,0)  = % 15.8e\n", fast_atan2f(1,0));
    printf ("atan2(-1,0) = % 15.8e\n", fast_atan2f(-1,0));
    printf ("atan2(0,1)  = % 15.8e\n", fast_atan2f(0,1));
    printf ("atan2(0,-1) = % 15.8e\n", fast_atan2f(0,-1));
    return EXIT_SUCCESS;

上述程序的输出应该类似于以下内容：

testing atan2() with 1000000000 test vectors
max rel err =  3.53486939e-005

atan2(1,0)  =  1.57079637e+000
atan2(-1,0) = -1.57079637e+000
atan2(0,1)  =  0.00000000e+000
atan2(0,-1) =  3.14159274e+000

【讨论】：

Godbolt 链接显示fast_atan2f() 内部循环成功矢量化：ICC 17、gcc 7.2、clang 5.0。我不知道如何设置 ARM 编译器以使用 Neon 进行自动矢量化。