atan2
とにかくこれを速くしたい。
場合分けはあとまわしにして、まずは、0 <= |y| <= xの範囲で考える。
前回のatanがそのまま使えて、
y /= x;
double y2 = y * y;
return y * (((256 * y2 + 5943) * y2 + 19250) * y2 + 15015)
/ (((1225 * y2 + 11025) * y2 + 24255) * y2 + 15015);
yを1に固定、xを1〜2に振ってみた。
そうすると、34nsと倍以上になった。
割り算が2回になったのがいけないのだろう。
演算の数は、加算6、乗算8、除算2である。これを[6,8,2]と表現しよう。
割り算を避けるには次のようにする。
これだと割り算が1回減って、高速化が期待できる。
実際のコードは、
double x2 = x * x;
double y2 = y * y;
return y * (15015 * x2 * x2 * x2 + 19250 * x2 * x2 * y2
+ 5943 * x2 * y2 * y2 + 256 * y2 * y2 * y2)
/ (x * (15015 * x2 * x2 * x2 + 24255 * x2 * x2 * y2
+ 11025 * x2 * y2 * y2 + 1225 * y2 * y2 * y2));
演算の回数は[6,28,1]で、15.9nsとなった。
掛け算がかなり多くなったが、次のようにすると減らせる。
double x2 = x * x;
double y2 = y * y;
return y * (((15015 * x2 + 19250 * y2) * x2
+ 5943 * y2 * y2) * x2 + 256 * y2 * y2 * y2)
/ (x * (((15015 * x2 + 24255 * y2) * x2
+ 11025 * y2 * y2) * x2 + 1225 * y2 * y2 * y2));
演算の回数は[6,22,1]で、14.6ns
さらにこうすると減らせる。
double x2 = x * x;
double y2 = y * y;
double y4 = y2 * y2;
double y6 = y2 * y4;
double a = 15015 * x2;
return y * (((a + 19250 * y2) * x2 + 5943 * y4) * x2 + 256 * y6)
/ (x * (((a + 24255 * y2) * x2 + 11025 * y4) * x2 + 1225 * y6));
演算の回数は[6,17,1]で、処理時間は上とほとんど変わらなかった。
場合分けを入れると、
double rotate_angle;
if(abs(y) <= abs(x)) {
if(0 < x) // [-π/4,π/4]
rotate_angle = 0;
else { // [-π,-3π/4], [3π/4,π]
y = -y;
x = -x;
rotate_angle = 0 <= y ? -M_PI : M_PI;
}
}
else {
double tmp = x;
if(0 < y) { // (π/4,3π/4)
x = y;
y = -tmp;
rotate_angle = M_PI / 2;
}
else { // (-3π/4,-π/4)
x = -y;
y = tmp;
rotate_angle = -M_PI / 2;
}
}
double x2 = x * x;
double y2 = y * y;
return y * (((15015 * x2 + 19250 * y2) * x2
+ 5943 * y2 * y2) * x2 + 256 * y2 * y2 * y2)
/ (x * (((15015 * x2 + 24255 * y2) * x2
+ 11025 * y2 * y2) * x2 + 1225 * y2 * y2 * y2))
+ rotate_angle;
これで35ns。
場合分けのほうが遅いのか。
