Calculate atan2 without std functions or C99

Question

I am calculating angles from a 3-axis accelerometer, but my compiler doesn't have a atan or atan2 function. It has a reserved memory slot, but it calls a function i can't find in any files.

My compiler is Keil µVision 4 running the ARMCC compiler. The compiles has the file math.h, but the function is extern and doesn't exist:

  extern _ARMABI double atan2(double /*y*/, double /*x*/);

Is there a lib or function I can include that has the function arctan implemented? Or is there an alternative function to calculate angles from accelerometer? I need full 3-axis calibration of the angles.

Edit: I was hoping to avoid a table full of pre-calculated values.

I can't verify this but if you are using the GNU ARM compiler with uVision then there should be an option in the IDE to use the math libraries as described [here](http://www.keil.com/support/man/docs/uv4/uv4_dg_armld.htm). You probably also need to be using the standard library instead of the micro library which excludes some floating point functions. — tinman, Aug 13 '12 at 09:38
This function definitely is present in the libraries that come with the ARMCC compiler, and should even work with microlib as the math functions are the same (the difference is in the low level software floating point functions). — Al Grant, Dec 31 '12 at 17:54

score 12 · Answer 1 · edited Apr 27 '22 at 11:06

The following code uses a rational approximation to get the arctangent normalized to the [0 1) interval (you can multiply the result by Pi/2 to get the real arctangent)

normalized_atan(x) ~ (b x + x^2) / (1 + 2 b x + x^2)

where b = 0.596227

The maximum error is 0.1620º

#include <stdint.h>
#include <math.h>

// Approximates atan(x) normalized to the [-1,1] range
// with a maximum error of 0.1620 degrees.

float normalized_atan( float x )
{
    static const uint32_t sign_mask = 0x80000000;
    static const float b = 0.596227f;

    // Extract the sign bit
    uint32_t ux_s  = sign_mask & (uint32_t &)x;

    // Calculate the arctangent in the first quadrant
    float bx_a = ::fabs( b * x );
    float num = bx_a + x * x;
    float atan_1q = num / ( 1.f + bx_a + num );

    // Restore the sign bit
    uint32_t atan_2q = ux_s | (uint32_t &)atan_1q;
    return (float &)atan_2q;
}

// Approximates atan2(y, x) normalized to the [0,4) range
// with a maximum error of 0.1620 degrees

float normalized_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;
}

In case you need more precision, there is a 3rd order rational function:

normalized_atan(x) ~ ( c x + x^2 + x^3) / ( 1 + (c + 1) x + (c + 1) x^2 + x^3)

where c = (1 + sqrt(17)) / 8

which has a maximum approximation error of 0.00811º

To fix an issue near (0,0) I changed one line to add .0001 which I presume is not going to cause a problem? `float atan_1q = num / (x * x + bxy_a + num + (float).0001); // .0001 to fix dbz and nan issue` — , Aug 06 '17 at 01:46
+1 These functions are well done, and even avoid branches. Do you have a reference on where you got the rational functions and the implementation? — Greg Allen, Dec 18 '19 at 20:06

score 9 · Accepted Answer · edited May 23 '17 at 12:18

9

Its not very difficult to implement your own arctan2. Convert arctan2 to arctan using this formula. And you can then calculate arctan using this infinite series. If you sum sufficient number of terms of this infinite series, you will get very close to what the library function arctan2 does.

Here is one similar implementation for exp() that you could use as a reference.

edited May 23 '17 at 12:18

Community

1
1

answered Aug 13 '12 at 08:51

Pavan Manjunath

27,404
12
99
125

5

Those infinite series are useful for exploring the mathematical properties of functions, but they are not suitable for computation. They are slow to converge, requiring many terms to get an accurate result, and are hard to calculate in floating-point without accumulating errors. Typically, math libraries compute these functions by using minimax polynomials, which are designed to minimize the maximum error using a few terms. – Eric Postpischil Aug 13 '12 at 11:19
2

Additionally, functions are partitioned into intervals for various reasons. arctan should be partitioned into |x| < 1 and |x| > 1 (|x| = 1 can be in either partition). For |x| < 1, arctan(x) is approximated with a minimax polynomial. For |x| > 1, arctan(x) may be approximated with a minimax polynomial using 1/x as input. – Eric Postpischil Aug 13 '12 at 11:22
Will that be an efficient implementation?. Which is the algorithm used in GCC libraries? – captain Aug 30 '12 at 05:00
@Jay Its not highly efficient but is very easy to implement. I am not sure about `libc` libraries but check Eric's answer for one of the open source implementation which extensively uses pre-defined constants to speed up calculations – Pavan Manjunath Aug 30 '12 at 06:13
@PavanManjunath Reason for not using Taylor series is that the error grows polynomially. A better approximation would be to use Remez exchange algorithm to find a polynomial (of same order) that gives bounded absolute error within the required interval, but no constraint outside the interval. – rwong Mar 13 '14 at 21:48

score 4 · Answer 3 · answered Aug 13 '12 at 11:25

4

There is an open source atan implementation here.

answered Aug 13 '12 at 11:25

Eric Postpischil

195,579
13
168
312

3

Apple has an open source portal? I never would have thought. +1 – andand Aug 13 '12 at 13:09

score 0 · Answer 4 · answered Aug 13 '12 at 08:33

0

The actual implementations of the math functions (or stubs to the HWFPU if one exists) should be in libm. With GCC this is indicated by passing -lm to the compiler, but I don't know how it is done with your specific tools.

answered Aug 13 '12 at 08:33

Ignacio Vazquez-Abrams

776,304
153
1,341
1,358

-lm isn't supported by Keil µVision. "Unrecognized option" – Jeffa Aug 13 '12 at 08:46
2

As far as I know Keil is a completely different compiler which has nothing to do with GCC whatsoever. – Lundin Aug 13 '12 at 08:59
@PavanManjunath: With GCC and similar compilers, “-l” does not mean a linker option. It means a library. “-l” means a library named “lib.” for appropriate suffixes, such as “.a” or “.dylib”, depending on the compiler version and the target platform. Libraries listed on the command line are passed to the linker, of course, by necessity. The GCC switch for passing a linker option is “-X”. – Eric Postpischil Aug 13 '12 at 12:40
@EricPostpischil. Yup you are specific and right! :) Thanks for the correction. I removed my ambiguous comment. – Pavan Manjunath Aug 13 '12 at 12:43

Calculate atan2 without std functions or C99

4 Answers4

Linked