您想要檢查的第一件事是您的編譯器是否能夠在應用於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
通常你會通過製作版本,採用矢量(S)作爲輸入矢量化數學函數。這通常比嘗試使用SIMD加速單個輸入的評估更有效。 –
@PeterCordes通常你可以得到gcc或clang的自動載入器給你一個寫得很好的近似向量化。由於缺少ieee754一致性(缺少非正常數字),只有ARM NEON纔將其擰緊。 – EOF
@EOF雖然這確實有用嗎?通常一個標量實現可能會稍微分散一些,或者使用表查找。 LUT可能仍然是一個勝利,但gcc不會自動矢量化解包/查找/重新包裝。例如,請參閱此AVX2'__m256d' log2函數https://stackoverflow.com/a/45898937/224132,它使用收集指令和只有一個小的多項式。 –