Approximating atan without a math library

Question

I've been googling around for a solution to this problem. I've seen a number of ways to calculate atan(theta) for any -1 <= theta <= 1, but I am not sure what to do when theta is bigger or smaller than those bounds.

I assume I need to add or subtract multiples of pi to offset it? Is this line of thinking correct?

Currently I have:

double my_atan(double x)
{
    return x - (x*x*x)/3 + (x*x*x*x*x)/5;
}

Which is using the taylor series.

And for the following code,

int x;
for (x=0; x<M_PI*2; x++)
{
    printf("Actual: %f\n",    atan(x));
    printf("Approx: %f\n", my_atan(x));
    printf("\n");
}

It quickly loses control (as expected, as it's out of range):

Actual: 0.000000
Approx: 0.000000

Actual: 0.785398
Approx: 0.866667

Actual: 1.107149
Approx: 5.733333

Actual: 1.249046
Approx: 42.600000

Actual: 1.325818
Approx: 187.466667

Actual: 1.373401
Approx: 588.333333

Actual: 1.405648
Approx: 1489.200000

Not pictured here, but the output is fairly accurate when theta is within the appropriate range.

So my question is what steps exactly need to be taken to the my_tan function to allow it to support wider bounds?

Been staring at this for a while and so any guidance that can be offered would be much appreciated

Lookup CORDIC for an efficient algorithm or implement the Taylor series as outline on the corresponding Wikipedia article (it only terminates for |x| < 1 so you need to scale the argument appropriately). — fuz, Feb 03 '16 at 17:41
Taylor is bad for this due to poor convergence rate and convergence radius. There are better techniques: best discovered by searching, but a Newton Raphson I believe works well. — Bathsheba, Feb 03 '16 at 18:11
keyword: cordic atan http://www.uio.no/studier/emner/matnat/ifi/INF5430/v12/undervisningsmateriale/dirk/Lecture_cordic.pdf — user3528438, Feb 03 '16 at 19:17
For a single-precision implementation, see [this question](http://stackoverflow.com/questions/26692859/best-machine-optimized-polynomial-minimax-approximation-to-arctangent-on-1-1) — njuffa, Feb 04 '16 at 01:52
the domain of arctan function is not the angles, its (-inf,inf). Get used to use the syntax 'atan(x) = theta' — crbah, Feb 15 '16 at 10:10

Gerhard Stein · Accepted Answer · 2016-02-03T20:21:13.090

Let me complete your example and talk about some things that might help:

#include <stdio.h>
#include <math.h>


double my_atan(double x)
{
    return x - (x*x*x)/3 + (x*x*x*x*x)/5 - (x*x*x*x*x*x*x)/7;
}

int main()
{
    double x;
    for (x=0.0; x<M_PI*2; x+=0.1)
    {
        printf("x: %f\n",    x);
        printf("Actual: %f\n",    atan(x));
        printf("Approx: %f\n", my_atan(x));
        printf("\n");
    }
}

The term int x is an integer and you are approximating large angles. try to use doubles here you get no conversioncast errors.

Now to your problem. The Taylor series you use only works if your |x| < 1.

Taylor Series get inaccurate the more you get away from a given point or in your case zero (0+x).

That series works well up to pi/4, even at that point it is very inaccurate but larger values get very bad. So for smaller angles it works well enough.

Shouldn't that be a 7 instead of 5? – Vicente Cunha Feb 03 '16 at 18:29 — Vicente Cunha, Feb 03 '16 at 18:29

Doug Currie · Answer 2 · 2016-02-03T18:36:38.823

3

The solution I used for a fixed point library was a minimax approximation using the Remez algorithm. Even there, I used a different set of coefficients for three ranges: 0 to 0.5; 0.5 to 0.75, and 0.75 to 1. With that breakdown I was able to get 1 ULP accuracy.

Then you need good argument reduction to get the argument in range. In my case I used a good reciprocal function for arguments over 1. Here are the identities:

atan(-x) == -atan(x)

atan(x) == pi/2 - atan(1/x) // for x > 1

There is a nice blog post on Taylor vs. Remez approximations here; at that site is also a Remez toolkit for finding the coefficients you need.

edited Feb 03 '16 at 18:36

answered Feb 03 '16 at 18:26

Doug Currie

40,708
1
95
119

Good thoughts about the issues +1. Not impressed with the blog post as it centers on absolute error over the interval and not relative error which is akin to `ULP` measurement. – chux - Reinstate Monica Feb 03 '16 at 22:06
@chux, about halfway through the lolremez toolkit tutorial there is a "switch to relative error" http://lolengine.net/wiki/doc/maths/remez/tutorial-relative-error – Doug Currie Feb 03 '16 at 23:21
@DougCurrie still missing one identity to maintain performance and accuracy as if `x>0.8` the Taylor series approach starts to converge more and more badly (the closer it gets to `1`). I am using `atan(x)-atan(y)=atan((x-y)/(1+x.y))` see [my C++ atan,asin](https://stackoverflow.com/a/50894477/2521214) – Spektre Jun 17 '18 at 08:18

Approximating atan without a math library

2 Answers2