Approximating inverse trigonometric funcions

Approximating inverse trigonometric funcions

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.

Answer by amit for Approximating inverse trigonometric funcions

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)

Answer by aero for Approximating inverse trigonometric funcions

Maybe some kind of intelligent brute force like newton rapson.

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

Answer by nulvinge for Approximating inverse trigonometric funcions

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:

Answer by njuffa for Approximating inverse trigonometric funcions

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

Answer by Jason S for Approximating inverse trigonometric funcions

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?

Answer by Eugeny89 for Approximating inverse trigonometric funcions

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.

Answer by rcor for Approximating inverse trigonometric funcions

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   #include     // 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?

Answer by Fnord for Approximating inverse trigonometric funcions

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.

Answer by florin for Approximating inverse trigonometric funcions

For Trigonometric Functions, special for inverse trigonometric functions( polinomial aproximation for arcsin(x)=x+x^4/3 ( For x <91/100).), then you need to use Newton method for solve sin(x)=a.

If you know sin(x) or cos(x) or tan(x) you can find sin(x/2),cos(x/2),and tan(x/2).then you fiind sin(x/4),then sin(x/8)......then sin(x/2^n),cos(x/2^n),tan(x/2^n).but x/2^n tend fi zero when n tend at infinite.in practice n=20( max ) will be enought.

After this use Taylor expansion( Mac Laurin series):arcsin(x)=x+x^3/6+3/40x^5+5/112x^7+35/1152x^9.....arctan(x)=x-x^3/3+x^5/5-x^7/7+x^9/9....

Answer by florin for Approximating inverse trigonometric funcions

arcsin(x)=pi/2-sqrt(1-x) (1.5707288-0.2121144 x+0.074261 x^2-0.0187293 x^3)) max error 0.0000676

From Handbook of Mathematical Functions by M. Abramowitz and I.A. Stegun, Ed.

nVidia reference implementation:

• http://http.developer.nvidia.com/Cg/asin.html

Answer by florin for Approximating inverse trigonometric funcions

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).

Answer by florin for Approximating inverse trigonometric funcions

1.221253148 x^8-2.894350729 x^7+3.072144547 x^6-1.671552458 x^5+0.572199725 x^4+0.0606801186 x^3+0.01003401398 x^2+0.999636250403 x=arcsin(x).(max error 0.89 second of a degree).Aproximation is for x in range 0,sqrt(1/2).For x bigger that sqrt(1/2) use arcsin(x)=pi/2-arcsin(sqrt(1-x^2)).(where arcsin is approximed by given eight degree polynomial)

Fatal error: Call to a member function getElementsByTagName() on a non-object in D:\XAMPP INSTALLASTION\xampp\htdocs\endunpratama9i\www-stackoverflow-info-proses.php on line 72