尾数为11位的atan2近似(用SSE2)和ARM(用vfpv4霓虹灯)

Question

我试图用尾数实现11位精度的快速atan2(浮点数)。atan2实现将用于图像处理。因此，最好用SIMD指令(针对x86的impl (用SSE2)和ARM(用vpfv4霓虹灯)来实现)。

现在，我使用chebyshev多项式逼近(https://jp.mathworks.com/help/fixedpoint/examples/calculate-fixed-point-arctangent.html)。

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

• 矢量化多项式逼近
  • chebyshev多项式(现已实现)

​

• 标量查找表(+线性插值)
• 矢量化查找表(这可能吗？)我不熟悉霓虹灯，所以我不知道一些指令，比如霓虹灯中的VGATHER )
• 矢量化CORDIC方法(这可能很慢，因为它需要12+旋转迭代才能获得11位的精度)

还有个好主意？

Answer 1

首先要检查的是，当编译器应用于float数组时，编译器是否能够将float矢量化。这通常需要至少一个较高的优化级别，例如-O3，并可能指定一种放松的“快速数学”模式，在这种模式下，errno处理、取消以及诸如无穷大和NaNs等特殊输入在很大程度上被忽略。使用这种方法，精确性可能远远超出了所需的范围，但在性能方面，很难超过经过精心优化的库实现。

接下来要尝试的是编写一个简单的标量实现，并使其具有足够的准确性，并让编译器将其矢量化。通常，这意味着避免使用非常简单的分支，这些分支可以通过if-转换转换为无分支代码。这样的代码的一个例子是下面所示的fast_atan2f()。使用作为icl /O3 /fp:precise /Qvec_report=2 my_atan2f.c调用的Intel编译器，这将被成功地向量化：my_atan2f.c(67): (col. 9) remark: LOOP WAS VECTORIZED.通过反汇编对生成的代码进行双重检查，显示fast_atan2f()已被内联并矢量化，使用的是*ps风格的SSE指令。

如果所有这些都失败了，您可以手工将fast_atan2()的代码转换为特定于平台的SIMD本质，这不应该那么困难。不过，我没有足够的经验来迅速完成这项工作。

我在这段代码中使用了一个非常简单的算法，它是一个简单的参数约简，在0,1中生成一个约简参数，然后是一个极小极大多项式逼近，最后一步将结果映射回完整的圆。核近似用Remez算法生成，并采用二阶Horner格式进行计算.融合乘法加法(FMA)可以在可用的情况下使用.为了性能起见，不处理涉及无穷大、NaNs或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