Atan2 Faster Approximation

Question

I am using atan2(y, x) for finding the polar angle from the x-axis and a vector which contains the point (x,y) for converting Cartesian coordinates to polar coordinates. But, in my program which will be used for calculations to move a robot, atan2 is a very computationally expensive procedure. What are ways to approximate atan2's result with a fast calculation with error bounds of 1 degree.

Wikipedia reference for atan2: https://en.wikipedia.org/wiki/Atan2

Here is some sample output:

atan2(1, 0) --> pi/2 atan2(-1, 0) --> -pi/2 atan2(0, 1) --> 0 atan2(0, -1) --> pi

Answer 1

Denote $$ \operatorname{atan2}(y,x) = \arctan z, \quad \mbox{where} \quad z = \dfrac{y}{x}. $$

Note that (in degrees) if $z>1$, then $$ \arctan z = 90^\circ - \arctan \dfrac{1}{z}, $$ and that $$ \arctan(-z) = -\arctan ~ z. $$

So, it is enough to write approximation for $z\in [0,1]$.

This simple approximation $$ \arctan z \approx z\Bigl(45^\circ - 15.66^\circ(z-1)\Bigr)\tag{1} $$ gives us error $\pm0.22^\circ$.

---

This source gives us not bad approximation. I just translated it to degrees:

$$ \arctan z \approx z\Bigl(45^\circ- \bigl(z - 1\bigr)\bigl(14^\circ + 3.83^\circ z\bigr)\Bigr), \qquad 0\le z \le 1;\tag{2} $$ error is $\pm 0.09^\circ$.

Authors of the paper said that this formula is $3\times$ faster than standard one.

Answer 2

The pseudo code below computes r = atan2 (y,x) in radians. Using IEEE-754 single precision for all variables, the maximum relative error observed for 100 million trials was $6.14(10^{-4})$, thus complying with the requirement that the relative error be $\lt \frac{1}{360} \approx 2.77(10^{-3})$.

The algorithm does not deal with special cases such as $x=0, y=0$, nor does it consider special IEEE-754 floating-point operands such as infinities and NaN.

This is based on a simple argument reduction using the identities $arctan (z) = \frac{\pi}{2} - arctan (\frac{1}{z})$ if $z \gt 0$ and $arctan (-z) = -arctan(z)$ followed by a polynomial minimax approximation on $[0,1]$, with final result mapped back to $[-\pi,\pi]$ based on the octant of the original vector. The polynomial is evaluated using the Horner scheme. The minimax approximation was computed using the Remez algorithm.

a := min (|x|, |y|) / max (|x|, |y|)
s := a * a
r := ((-0.0464964749 * s + 0.15931422) * s - 0.327622764) * s * a + a
if |y| > |x| then r := 1.57079637 - r
if x < 0 then r := 3.14159274 - r
if y < 0 then r := -r