2017-09-14 263 views
3

我试图以11位精度尾数实现快速atan2(浮点数)。atan2实现将用于图像处理。 所以用SIMD指令(impl瞄准x86(带有SSE2)& ARM(带有vpfv4 NEON))可能会更好。atan2近似11位尾数在x86(与SSE2)和ARM(与vfpv4 NEON)

现在,我使用切比雪夫多项式近似(https://jp.mathworks.com/help/fixedpoint/examples/calculate-fixed-point-arctangent.html)。

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

  • 矢量多项式逼近
    • 切比雪夫多项式(现已实现)
  • 标查找表(+线性插值)
  • 量化的查找表(这是可能的吗?我不熟悉NEON,所以我不知道NEON中的一些指令,如VGATHER)
  • 矢量化的CORDIC方法(这可能是缓慢的,因为它需要12次旋转迭代才能获得11位准确度)

其他好主意?

+2

通常你会通过制作版本,采用矢量(S)作为输入矢量化数学函数。这通常比尝试使用SIMD加速单个输入的评估更有效。 –

+0

@PeterCordes通常你可以得到gcc或clang的自动载入器给你一个写得很好的近似向量化。由于缺少ieee754一致性(缺少非正常数字),只有ARM NEON才将其拧紧。 – EOF

+0

@EOF虽然这确实有用吗?通常一个标量实现可能会稍微分散一些,或者使用表查找。 LUT可能仍然是一个胜利,但gcc不会自动矢量化解包/查找/重新包装。例如,请参阅此AVX2'__m256d' log2函数https://stackoverflow.com/a/45898937/224132,它使用收集指令和只有一个小的多项式。 –

回答

2

您想要检查的第一件事是您的编译器是否能够在应用于float的数组时向量化atan2f (y,x)。这通常需要至少一个较高的优化级别,如-O3,并且可能会指定一个宽松的“快速数学”模式,其中errno处理,非规范化和特殊输入(如infinities和NaN)在很大程度上被忽略。采用这种方法,准确度可能远远超出要求的范围,但在性能方面可能很难打败仔细调整过的库实现。

接下来要尝试的是编写一个具有足够准确性的简单标量实现,并让编译器对其进行矢量化。通常这意味着避免任何事情,但只能通过if转换转换为无分支代码的非常简单的分支。这种代码的一个例子是fast_atan2f(),如下所示。使用英特尔编译器,调用为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]中生成简化参数的简单参数约简,接着是最小多项式近似,最后一步将结果映射回完整的循环。 Core近似是用Remez算法生成的,并使用二阶霍纳方案进行评估。在可用的地方可以使用熔融多次加法(FMA)。为了表现,不处理涉及infinities,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 
+0

Godbolt链接显示'fast_atan2f()内循环:[ICC 17](https://godbolt.org/g/P7ztE3),[gcc 7.2](https://godbolt.org/g/bt8Xrs),[铿锵5.0](https:// godbolt .ORG /克/ 1WiEa8)。我不知道如何设置ARM编译器来实现与Neon的自动矢量化。 – njuffa