SSEで速度がアップすることがわかったので、さっそくatan2の計算に適用してみよう。
一つの関数で4つのatan2の計算を一度に行う。
場合わけの部分はそのまま残しておいて、有理関数の計算をSSEで行う（コードは下）。
環境が今までとちがっていて、PCが少し速くて、コンパイラはVC2008EE。あとから今までの環境でベンチマークしなおす。
SSE版は30nsが25nsとなった。
場合わけの部分はそのままであるし、代入が多いのであまり速くならない。
有理関数の計算部分で最も時間がかかっているのは割り算の部分だろう。
divps xmm4, xmm5 ; xmm4 /= xmm5
ここは逆数の近似を出すコマンドを使うことができる。
rcpps xmm6, xmm5 ; xmm6 = 1 / xmm5
mulps xmm4, xmm6 ; xmm4 *= xmm6
しかし、このコマンドは精度が悪い。1.5×2-12、すなわち3.7e-4程度である。
atan2の精度が5e-6程度だったのが、2e-4程度となった。
精度を上げるために、ここでニュートン法を使おう。
rcppsで得られた近似の逆数をx0、逆数を取る前の数値をaとする。
y = 1/xの(x0, 1/x0)における接線は、
これと、y = aの交点は、
より、
x0 = (1 + ε)/aとすると、
x = (1 - ε2)/aとなって、精度がアップする。
このアセンブラはこうなる。
rcpps xmm6, xmm5 ; x0 ≒ 1 / a
movups xmm1, xmm6
mulps xmm1, xmm1
mulps xmm1, xmm5 ; x2 = x0 * x0 * a
addps xmm6, xmm6 ; x0 += x0
subps xmm6, xmm1 ; x = 2x0 - x2
mulps xmm4, xmm6
これで、精度が元に戻り、24nsとなった。
double atan12(float *z, const float *y, const float *x) {
float rotate_angle[4];
float x2[4], y2[4], z2[4];for(int i = 0; i < 4; i++) {
if(x[i] + y[i] >= 0) {
if(x[i] >= y[i]) { // [-π/4,π/4]
x2[i] = x[i];
y2[i] = y[i];
rotate_angle[i] = 0;
}
else { // (π/4,3π/4)
x2[i] = y[i];
y2[i] = -x[i];
rotate_angle[i] = M_PI / 2;
}
}
else {
if(x[i] >= y[i]) { // (-3π/4,-π/4)
x2[i] = -y[i];
y2[i] = x[i];
rotate_angle[i] = -M_PI / 2;
}
else { // [-π,-3π/4], [3π/4,π]
y2[i] = -y[i];
x2[i] = -x[i];
rotate_angle[i] = 0 <= y2[i] ? -M_PI : M_PI;
}
}
}
float a0[4] = { 256.0f, 256.0f, 256.0f, 256.0f };
float a1[4] = { 5943.0f, 5943.0f, 5943.0f, 5943.0f };
float a2[4] = { 19250.0f, 19250.0f, 19250.0f, 19250.0f };
float a3[4] = { 15015.0f, 15015.0f, 15015.0f, 15015.0f };
float b0[4] = { 1225.0f, 1225.0f, 1225.0f, 1225.0f };
float b1[4] = { 11025.0f, 11025.0f, 11025.0f, 11025.0f };
float b2[4] = { 24255.0f, 24255.0f, 24255.0f, 24255.0f };
__asm {
movups xmm0, [y2]
movups xmm1, [x2]
movups xmm2, xmm0
movups xmm3, xmm1
mulps xmm2, xmm2 ; y^2
mulps xmm3, xmm3 ; x^2
movups xmm6, [a3]
movups xmm4, xmm3
mulps xmm4, xmm6 ; xmm4 = x^2 * 15015
movups xmm5, xmm4 ; xmm5 = xmm4
movups xmm6, [a2]
mulps xmm6, xmm2
addps xmm4, xmm6 ; xmm4 += a1 * y^2
mulps xmm4, xmm3 ; xmm4 *= x^2
movups xmm6, [b2]
mulps xmm6, xmm2
addps xmm5, xmm6 ; xmm5 += a1 * y^2
mulps xmm5, xmm3 ; xmm5 *= x^2
movups xmm6, [a1]
mulps xmm6, xmm2
mulps xmm6, xmm2
addps xmm4, xmm6
mulps xmm4, xmm3
movups xmm6, [b1]
mulps xmm6, xmm2
mulps xmm6, xmm2
addps xmm5, xmm6
mulps xmm5, xmm3
movups xmm6, [a0]
mulps xmm6, xmm2
mulps xmm6, xmm2
mulps xmm6, xmm2
addps xmm4, xmm6
mulps xmm4, xmm3
movups xmm6, [b0]
mulps xmm6, xmm2
mulps xmm6, xmm2
mulps xmm6, xmm2
addps xmm5, xmm6
mulps xmm5, xmm3
mulps xmm4, xmm0 ; xmm4 *= y
mulps xmm5, xmm1 ; xmm5 *= x
#if 1 // 逆数の近似計算
rcpps xmm6, xmm5 ; xmm6 = 1 / xmm5
#if 0 // 精度アップ
movups xmm1, xmm6
mulps xmm1, xmm1
mulps xmm1, xmm5 ; xmm1 = xmm6^2 * xmm5
addps xmm6, xmm6 ; xmm6 += xmm6
subps xmm6, xmm1 ; xmm6 -= xmm1
#endif
mulps xmm4, xmm6
#else
divps xmm4, xmm5 ; xmm4 /= xmm5
#endif
movups z2, xmm4 ; z = xmm4
}
z[0] = z2[0] + rotate_angle[0];
z[1] = z2[1] + rotate_angle[1];
z[2] = z2[2] + rotate_angle[2];
z[3] = z2[3] + rotate_angle[3];
}void bench_atan2() {
int t = timeGetTime();
double sum = 0;
#if 1
float x[4] = { -0.5f, -0.5f + dx, -0.5f + dx * 2, -0.5f + dx * 3 };
float y[4] = { 1.0f, 1.0f, 1.0f, 1.0f };
float z[4];
float max_diff = 0;
for(__int64 i = 0; i < N / 4; i++) {
atan12(z, y, x);
float diff = abs(z[0] - atan2(y[0], x[0]));
if(diff > max_diff)
max_diff = diff;
for(int k = 0; k < 4; k++) {
sum += z[k];
x[k] += dx * 4;
}
}
fprintf(stdout, "max diff : %e\n", max_diff);
#else
double x = -0.5;
double y = 1.0;
for(__int64 i = 0; i < N; i++) {
sum += atan10(y, x);
x += dx;
}
#endif
fprintf(stdout, "%.3fs\n", (timeGetTime() - t) * 1e-3);
fprintf(stdout, "%e\n", sum);
