Even during 10 years in a past life working in a group that designed FIR filters, I never grappled with the Remez exchange algorithm hands-on. This week I finally found the excuse I needed.
A while ago, I needed a crude but very fast approximation to arctan. In fact it was atan2, but I’ll start with the easier problem first. With a release of the application imminent, I wanted to credit the original source I got the formula from, but I couldn’t find it any longer. So I set about recomputing the formula for myself.
The version I found a while ago looked like this:
atan(x) ≈ (0.97179803008 * x) - (0.19065470515 * x**3)
for x in the range [-1,+1], i.e. the result is in the range [-pi/4,+pi/4]. The maximum error in the result is just less than the equivalent of 0.3 degrees.
[If the person behind these coefficients gets in touch, I’ll be more than happy to give credit and link to the original work here.]
Having now got the Remez exchange algorithm working, I found a slightly improved version is
atan(x) ≈ (0.97239 * x) - (0.19195 * x**3)
which has a worst-case error of 0.2837 degrees. Even better, a plot shows the 6 extrema all have this absolute error, as expected.
Moving to the original atan2 problem, the requirement was rather unusual. Starting with the x and y arguments as integers, the answer is to be expressed as a 32-bit integer, where 2**32 represents 360 degrees. This allows subsequent rotations to be done using integer add and subtract, with the convenient property that multiples of 360 degrees are lost by the modulo wraparound at 2**32. The resulting function looks like this, with the coefficients as shown earlier:
#define K1 (0.97239)
#define K2 (0.19195)
#define K3 ((float)((double)(1<<29) * 4.0 * K1 / M_PI))
#define K4 ((float)((double)(1<<29) * 4.0 * K2 / M_PI))
int32_t atan2i_approx2(int32_t y, int32_t x)
{
int32_t u, v, s;
uint32_t w;
float f, f2, g;
int32_t r;
float fx, fy;
float tt, bb;
fx = (float) x;
fy = (float) y;
u = x + y;
v = x - y;
w = (u ^ v) >> 31;
tt = 
? -fx : fy;
bb = ? fy : fx;
f = tt / bb;
f2 = f * f;
g = f * (K3 - K4*f2);
s = (((u >> 30) & 3) ^ ((v >> 31) & 1)) << 30;
r = s + (int32_t) g;
return r;
}
https://rc0rc0.wordpress.com/2013/06/05/minimax-approximation-to-arctan-atan-atan2/
If a fast reciprocal instruction is available, rational approximations can be a good choice. The following code computes the arctangent normalized to the [0 1) interval with a maximum error of 0.1620º
normalized_atan(x) ~ (b x + x^2) / (1 + 2 b x + x^2)
where b = 0.596227
// Approximates atan2(y, x) normalized to the [0,4) range
// with a maximum error of 0.1620 degrees
float norm_atan2( float y, float x )
{
static const uint32_t sign_mask = 0x80000000;
static const float b = 0.596227f;
// Extract the sign bits
uint32_t ux_s = sign_mask & (uint32_t &)x;
uint32_t uy_s = sign_mask & (uint32_t &)y;
// Determine the quadrant offset
float q = (float)( ( ~ux_s & uy_s ) >> 29 | ux_s >> 30 );
// Calculate the arctangent in the first quadrant
float bxy_a = ::fabs( b * x * y );
float num = bxy_a + y * y;
float atan_1q = num / ( x * x + bxy_a + num );
// Translate it to the proper quadrant
uint32_t uatan_2q = (ux_s ^ uy_s) | (uint32_t &)atan_1q;
return q + (float &)uatan_2q;
}
Really nice, and thanks so much for the idea!
Here are some better coefficients that give a maximum error of 0.142781 degrees! Additionally, pi/4 gives the exact answer, so doing the pi/2-atan(1/x) trick gives continuous results over the entire function! Also it is a good idea to have the equation in the form below, for 1 less multiplication:
x*(0.983758618f-0.220587f*x*x)
I am using this in my code now.
~Rodol Phito
[…] the sine and cosine the web retrieved many interesting pages also for the arctan (e.g. here, here, here, here). The algorithms generally reduce the problem to calculating the arctan in the range . The […]