Approximating inverse trigonometric functions

问题:

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.

回答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.

回答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

回答3:

you might want to use approximation: use an infinite series until the solution is close enough for you.

for example:
arcsin(z) = Sigma((2n!)/((2^2n)*(n!)^2)*((z^(2n+1))/(2n+1))) where n in [0,infinity)

回答4:

Maybe some kind of intelligent brute force like newton rapson.

So for solving asin() you go with steepest descent on sin()

回答5:

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:

回答6:

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º

回答7:

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.

回答8:

Use a polynomial approximation. Least-squares fit is easiest (Microsoft Excel has it) and Chebyshev approximation is more accurate.

This question has been covered before: How do Trigonometric functions work?

回答9:

Only continous functions are approximable by polynomials. And arcsin(x) is discontinous in point x=1.same arccos(x).But a range reduction to interval 1,sqrt(1/2) in that case avoid this situation. We have arcsin(x)=pi/2- arccos(x),arccos(x)=pi/2-arcsin(x).you can use matlab for minimax approximation.Aproximate only in range [0,sqrt(1/2)](if angle for that arcsin is request is bigger that sqrt(1/2) find cos(x).arctangent function only for x<1.arctan(x)=pi/2-arctan(1/x).