How to implement arctan(x) without using standard library?

[Gerald]

Recently, my program need to be run in embeded enviroment, and I cann't
use standard library. But I need to use arctan(x), so I implement it
like the following:

inline double pow(double x, size_t n) {
if (n == 0)
return 1;
else if (n % 2 == 0)
return pow(x * x, n >> 1);
else
return x * pow(x * x, n >> 1);
}

inline double atan_Gerald(double x) {
register double temp1 = x >= 0 ? x : -x;
register double temp2 = temp1 <= 1.0 ? temp1 : (temp1 - 1) / (temp1 +
1);
register double sum = temp2;

for (register size_t i = 1; i != 6; ++i)
sum += (i % 2 ? -1 : 1) * pow(temp2, (i << 1) + 1) / ((i << 1) + 1);

if (temp1 > 1.0) sum += 0.785398;
return x >= 0 ? sum : -sum;
}

But I found it doesn't work efficiently like atan in standard library,
Would you please
help me to optimize my code or give me another suggestion?

Thank you!

Regards Gerald

 

[Gerald]

I made a mistake: the 'register' before double is useless, I forget to
delete it.

 

[Dietmar Kuehl]

Recently, my program need to be run in embeded enviroment, and I cann't
use standard library.

I don't think that this is a reasonable argument for not using the
standard library! In particular, I know that there is at least one
standard library targeted also at embedded systems: the Dinkumware
library is available in a version tailored to embedded systems
(although I personally don't see any justification at all in EC++...).
See <http://www.dinkumware.com/>. The prices are definitely in a
range where it comes much cheaper to buy the library than provide a
quality implementation of even a single trigonometric function.

But I need to use arctan(x), so I implement it
like the following:

I haven't implemented these functions myself and thus I might be
entirely wrong about how I think these functions are implemented.
Still, my understanding is that they are implemented very different
from how you tried to implement them. The basic approach is something
like this:

1. Translate the function argument(s) to a minimal basic area (this
is almost certainly not the correct term but I don't know the
correct term). My maths is rather rusty but I think for arctan(x)
the basic area is the range [0, 1]: all other values can be
derived from an appropriate argument:
- x < 0: arctan(x) = -arctan(-x)
- 1 < x: arctan(1/x) = pi/2 - arctan(x)
(actually, you already do this part...)
2. For the values in basic area you look up the subrange your value
is located in. Each subrange is associated with an approximation
function which is just a suitably parameterized polynomial of a
relatively small degree which can be evaluated reasonably fast.
The real tricky thing is to determine suitable subranges and the
parameters for their approximating polynomials. I have no idea
at all how this done.

 

[Earl Purple]

I haven't implemented these functions myself and thus I might be
entirely wrong about how I think these functions are implemented.
Still, my understanding is that they are implemented very different
from how you tried to implement them. The basic approach is something
like this:

1. Translate the function argument(s) to a minimal basic area (this
is almost certainly not the correct term but I don't know the
correct term). My maths is rather rusty but I think for arctan(x)
the basic area is the range [0, 1]: all other values can be
derived from an appropriate argument:
- x < 0: arctan(x) = -arctan(-x)
- 1 < x: arctan(1/x) = pi/2 - arctan(x)
(actually, you already do this part...)
2. For the values in basic area you look up the subrange your value
is located in. Each subrange is associated with an approximation
function which is just a suitably parameterized polynomial of a
relatively small degree which can be evaluated reasonably fast.
The real tricky thing is to determine suitable subranges and the
parameters for their approximating polynomials. I have no idea
at all how this done.

For good accuracy use cubic splines. For optimal speed use a bigger
caching table and linear progression between the points. (A cubic
spline has 4 parameters in each range so using the linear version you
can store in theory 5 times as many points).

 

[Earl Purple]

Of course you don't have to store the cubic spline parameters as they
can be worked out
from the cached value, but that would give you even slower performance.

 

[Michiel.Salters]

Recently, my program need to be run in embeded enviroment, and I cann't
use standard library.
But I need to use arctan(x)

The basic approach is something like this:

1. Translate the function argument(s) to a minimal basic area (this
is almost certainly not the correct term but I don't know the
correct term). My maths is rather rusty but I think for arctan(x)
the basic area is the range [0, 1]: all other values can be
derived from an appropriate argument:
- x < 0: arctan(x) = -arctan(-x)
- 1 < x: arctan(1/x) = pi/2 - arctan(x)
(actually, you already do this part...)

There are two kinds of reductions. The first one is to keep the tables
finite
(This is essential for e.g. sin(x), as x may be any value) and the
second is
just to reduce the table size. The latter is something to worry about
later.
(It's a classic time/space tradeoff you'd need to profile).

2. For the values in basic area you look up the subrange your value
is located in. Each subrange is associated with an approximation
function which is just a suitably parameterized polynomial of a
relatively small degree which can be evaluated reasonably fast.
The real tricky thing is to determine suitable subranges and the
parameters for their approximating polynomials. I have no idea
at all how this done.

Polynomials are very efficient in finding the first few digits, but may
be terrible
in finding the next. E.g. when you already have arctan(x)=0.23....
you'd need
to look in table 23 to get the next coefficients, which implies 100
tables.
That's nog good for embedded systems. (Of course, real systems probably
use nibbles not digits, so you'd need 256 tables after the first two
nibbles)

However, it seems a real arctan algorithm is a lot more complex:
http://www.electronics.dit.ie/staff/aschwarzbacher/research/ecs99.pdf

I wouldn't invent that myselves. Implementing the basic CORDIC seems
easy, though.

HTH,
Michiel Salters