Algorithm:
```
  1)  atan(x) = sign(x) * atan(|x|)

  2)  If |x| > 1 + 1/32, atan(|x|) = pi/2 - atan(1/|x|)

  3)  For 1/16 < |x| < 1 + 1/32, we find k such that: | |x| - k/16 | <= 1/32.
      Let y = |x| - k/16, then using the angle summation formula, we have:
    atan(|x|) = atan(k/16) + atan( (|x| - k/16) / (1 + |x| * k/16) )
              = atan(k/16) + atan( y / (1 + (y + k/16) * k/16 )
              = atan(k/16) + atan( y / ((1 + k^2/256) + y * k/16) )

  4)  Let u = y / (1 + k^2/256), then we can rewritten the above as:
    atan(|x|) = atan(k/16) + atan( u / (1 + u * k/16) )
              ~ atan(k/16) + (u - k/16 * u^2 + (k^2/256 - 1/3) * u^3 +
                              + (k/16 - (k/16)^3) * u^4) + O(u^5)
```

With all the computations are done in single precision (float), the
total of approximation errors and rounding errors is bounded by 4 ULPs.