重力を計算するところ、
float dx = diff_x(sx[i], x);
float dy = diff_x(sy[i], y);
float dist_square = dx * dx + dy * dy;
float d = dist_square * sqrt(dist_square);
ax += dx / d;
ay += dy / d;
ここで、sqrtを使っているが、実は、
sqrt(x) = x / rsqrt(x)
を使って計算しているらしい。rsqrtは見ての通り平方根の逆数。そういえば、x86にもこういう命令があって、精度は悪いが速いらしい。これを同じ理屈だろうか。そういうわけで、ここは次のようにすると少し速くなる。
float dist_square = dx * dx + dy * dy;
float rec_d = rsqrt(dist_square) / dist_square;
ax += dx * rec_d;
ay += dy * rec_d;
これで、2.95秒で16%ほど速くなった。
