最初にチェックしたいのは、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.生成されたコードを逆アセンブルでダブルチェックすると、*psフレーバーのSSE命令を使用してfast_atan2f()がインライン化およびベクトル化されます。
他のすべてが失敗した場合は、fast_atan2()のコードを手作業でプラットフォーム固有のSIMD組み込み関数に変換することができます。私はそれをすばやく行うには十分な経験がありません。
私は非常に単純なアルゴリズムをこのコードで使用しました。このアルゴリズムは、[0,1]に引数を減らし、その後にminimax多項式近似を行う単純な引数の削減と、結果を完全な円。コア近似はRemezアルゴリズムで生成され、2次ホーナー方式を使用して評価されます。 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
通常は、入力としてベクトル（複数可）かかりバージョンを作ることによって、数学関数をベクトル化すると思います。これは通常、SIMDを使用して単一入力の評価を高速化しようとするよりも効率的です。 –
@PeterCordes通常、gccまたはclangのオートベクトルを取得して、よく書かれた近似のベクトル化を行うことができます。 ieee754の適合性が失われてしまったために、この問題を解決するのはARM NEONだけです（デノーマル番号がありません）。 – EOF
@EOFこれは本当に機能しますか？しばしば、スカラーの実装はやや枝分かれするか、テーブルルックアップを使用します。 LUTはまだ勝つかもしれませんが、gccはunpack/lookup/repackを自動ベクトル化しません。たとえば、このAVX2 '__m256d' log2関数https://stackoverflow.com/a/45898937/224132を参照してください。この関数はギャザリング命令と小さな多項式のみを使用します。 –