atan2 approximation with x86 (SSE2 사용) 및 ARM (vfpv4 NEON 사용)

Question

11 비트 정확도의 빠른 atan2 (float)를 가수로 구현하려고합니다. atan2 구현은 이미지 처리에 사용됩니다. 따라서 SIMD 명령어 (x86 (SSE2 사용) & ARM (vpfv4 NEON 사용))를 사용하여 구현하는 것이 더 나을 수 있습니다.

지금은 Chebyshev 다항식 근사법 (https://jp.mathworks.com/help/fixedpoint/examples/calculate-fixed-point-arctangent.html)을 사용합니다.

나는 수동으로 벡터화 된 코드를 구현하고자합니다. SSE2 (또는 이상) & NEON 래퍼 라이브러리를 사용합니다. 나는 계획하거나

룩업 테이블 벡터화이

• 벡터화 다항식 근사
  • 체비 쇼프 다항식 (현재 구현) 구현 옵션
• 스칼라 룩업 테이블 (+ 선형 보간을) 시도 NEON에 익숙하지 않아 NEON의 VGATHER와 같은 지침을 알지 못합니다.
• 벡터화 된 CORDIC 방법 (11 비트 정확도를 얻기 위해 12 회 이상의 회전 반복이 필요하므로 느릴 수 있음)

기타 좋은 아이디어?

Answer 1

먼저 확인하고 싶은 것은 float 배열에 적용 할 때 컴파일러가 atan2f (y,x)을 벡터화 할 수 있는지 여부입니다. 대개 -O3과 같이 높은 최적화 수준이 필요하고 무 처리 및 NaN 등의 특수 입력 및 errno 처리, 비정규 및 특수 입력이 크게 무시되는 편안한 "빠른 연산"모드를 지정해야합니다. 이 방법을 사용하면 정확도가 필요한 것보다 훨씬 좋을 수 있지만 성능면에서 조심스럽게 조정 된 라이브러리 구현을 이길 수는 없습니다.

다음으로 시도해보십시오. 충분한 정확도로 간단한 스칼라 구현을 작성하고 컴파일러에서 벡터화하도록하는 것입니다. 일반적으로 이것은 if-conversion을 통해 분기없는 코드로 변환 될 수있는 매우 간단한 분기를 피하는 것을 의미합니다. 이러한 코드의 예는 아래에 표시된 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 flavor의 SSE 명령어를 사용하여 인라인되고 벡터화 된 것으로 나타납니다.

다른 모든 것이 실패하면 fast_atan2() 코드를 직접 플랫폼 고유의 SIMD 내장 함수로 변환 할 수 있습니다. 나는 그것을 빨리 할 충분한 경험이 없다.

나는이 코드에서 매우 간단한 알고리즘을 사용했다.이 알고리즘은 [0,1]에서 감소 된 인수를 생성하고, minimax 다항식 근사법을 따르는 단순한 인수 축소 및 최종 결과를 완전한 동그라미. 코어 근사는 Remez 알고리즘을 사용하여 생성되었으며 2 차 Horner 체계를 사용하여 평가됩니다. 사용 가능한 경우 Fused-multiply add (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