Approximating inverse trigonometric functions

Question

I have to implement asin, acos and atan in environment where I have only following math tools:

• sine
• cosine
• elementary fixed point arithmetic (floating point numbers are not available)

I also already have reasonably good square root function.

Can I use those to implement reasonably efficient inverse trigonometric functions?

I don't need too big precision (the floating point numbers have very limited precision anyways), basic approximation will do.

I'm already half decided to go with table lookup, but I would like to know if there is some neater option (that doesn't need several hundred lines of code just to implement basic math).

EDIT:

To clear things up: I need to run the function hundreds of times per frame at 35 frames per second.

Solution 1

Do you need a large precision for arcsin(x) function? If no you may calculate arcsin in N nodes, and keep values in memory. I suggest using line aproximation. if x = A*x_(N) + (1-A)*x_(N+1) then x = A*arcsin(x_(N)) + (1-A)*arcsin(x_(N+1)) where arcsin(x_(N)) is known.

Solution 2

In a fixed-point environment (S15.16) I successfully used the CORDIC algorithm (see Wikipedia for a general description) to compute atan2(y,x), then derived asin() and acos() from that using well-known functional identities that involve the square root:

asin(x) = atan2 (x, sqrt ((1.0 + x) * (1.0 - x)))
acos(x) = atan2 (sqrt ((1.0 + x) * (1.0 - x)), x)

It turns out that finding a useful description of the CORDIC iteration for atan2() on the double is harder than I thought. The following website appears to contain a sufficiently detailed description, and also discusses two alternative approaches, polynomial approximation and lookup tables:

http://ch.mathworks.com/examples/matlab-fixed-point-designer/615-calculate-fixed-point-arctangent

Solution 3

http://en.wikipedia.org/wiki/Inverse_trigonometric_functions#Expression_as_definite_integrals

You could do that integration numerically with your square root function, approximating with an infinite series:

Solution 4

It should be easy to addapt the following code to fixed point. It employs a rational approximation to calculate the arctangent normalized to the [0 1) interval (you can multiply it by Pi/2 to get the real arctangent). Then, you can use well known identities to get the arcsin/arccos from the 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 norm_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 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;
} 

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º

Solution 5

Submitting here my answer from this other similar question.

nVidia has some great resources I've used for my own uses, few examples: acos asin atan2 etc etc...

These algorithms produce precise enough results. Here's a straight up Python example with their code copy pasted in:

import math
def nVidia_acos(x):
    negate = float(x<0)
    x=abs(x)
    ret = -0.0187293
    ret = ret * x
    ret = ret + 0.0742610
    ret = ret * x
    ret = ret - 0.2121144
    ret = ret * x
    ret = ret + 1.5707288
    ret = ret * math.sqrt(1.0-x)
    ret = ret - 2 * negate * ret
    return negate * 3.14159265358979 + ret

And here are the results for comparison:

nVidia_acos(0.5)  result: 1.0471513828611643
math.acos(0.5)    result: 1.0471975511965976

That's pretty close! Multiply by 57.29577951 to get results in degrees, which is also from their "degrees" formula.

View more solutions

Comments

[Karoly Horvath]

and you can choose the starting point from a small lookup-table to speed up the calculation.

[Matěj Zábský]

Yeah, that's the table lookup I was talking about in OP. don't see a reason why I would calculated that at runtime, I would simoly bake the values into the program, so the actual asin calculation would just be an interpolation between two values.

[petrelharp]

From wikipedia, CORDIC doesn't even use the trig functions (neat!); I imagine what you did was a search; given sin(), cos() it seems Newton-Raphson or some such would be better? (Require fewer iterations, although the cost of the iterations would be different.)

[njuffa]

The reason I suggested looking into CORDIC is because it only requires fixed-point arithmetic. The most common use of CORDIC is probably for implementing sin / cos, that is how I first learned about it (in 1987). But quite a few other functions can be computed with CORDIC as well, such as atan2. Since I do not have any code lying around for computing atan2 with CORDIC I tried to find a website with enough detail that someone could base an implementation on it. The link I posted above was the best page I could find via a search engine in the space of a few minutes.