Atan2 Faster Approximation 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
 A: Denote
$$
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.
A: 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

