Computing ArcSin, ArcCos, ArcTan
Kye
Posts: 2,200
Hey guys,
I'm working on integer math trig rountines for the prop. So far I've got the sin, cos, and tan working perfectly by using an angle and a·radius to make the calculations.
As in... sin(pheta) is never greater than 1. But if I scale the value that comes out of the sin function provided in the manual by a radius (polar form) I can increase the resolution so that the output swings back and forth between the +-radius.
Now, I'm trying to do the arcSin and arcCos, arctan functions. But I'm not really sure what is the best approach to this.
So far I have:
Now, this code most likely does not work, but I'm trying to figure out what to do.
Basically I'm taking the yLength of a right triangle and then scaling it down by the polar form radius. This then produces the 16 bit 2's complement value the sin function in the manual outputs.
Then if the value is negative I make it positive so that I can look it up in the sin table.
.... Now here's the tricky part that I'm not sure what to do. I have the value that the sin table would output essentially, but I want to find the index where that value belongs. Is there any other way to find that index without searching throught the whole table?
So, anyone else worked with this and know what to do?
I think I have correct·the parts after finding the index which would be the angle in the table.
Thanks for your help,
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Nyamekye,
I'm working on integer math trig rountines for the prop. So far I've got the sin, cos, and tan working perfectly by using an angle and a·radius to make the calculations.
As in... sin(pheta) is never greater than 1. But if I scale the value that comes out of the sin function provided in the manual by a radius (polar form) I can increase the resolution so that the output swings back and forth between the +-radius.
Now, I'm trying to do the arcSin and arcCos, arctan functions. But I'm not really sure what is the best approach to this.
So far I have:
PUB arcSin(yLength, radius) | mirror '' 6 Stack Longs
'' ┌──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
'' │ Computes the arc sin of an angle. │
'' │ │
'' │ Returns the angle in which ((radius) * sin(angle)) = yLength. │
'' │ │
'' │ yLength - The opposite of the angle length of a right triangle. │
'' │ Radius - The radius (hypotenuse) of the polar plane triangle. │
'' └──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
yLength := ((yLength << 16) / ((radius <# 65535) #> 0))
mirror := (yLength < 0)
if(mirror)
-yLength
repeat result from 0 to 2049
if(yLength == word[noparse][[/noparse](result | constant($E000 >> 1)) << 1])
quit
if(mirror)
-result
return ((result << 3) / 91)
Now, this code most likely does not work, but I'm trying to figure out what to do.
Basically I'm taking the yLength of a right triangle and then scaling it down by the polar form radius. This then produces the 16 bit 2's complement value the sin function in the manual outputs.
Then if the value is negative I make it positive so that I can look it up in the sin table.
.... Now here's the tricky part that I'm not sure what to do. I have the value that the sin table would output essentially, but I want to find the index where that value belongs. Is there any other way to find that index without searching throught the whole table?
So, anyone else worked with this and know what to do?
I think I have correct·the parts after finding the index which would be the angle in the table.
Thanks for your help,
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Nyamekye,
Comments
if your not crunched for mem. space, could you just create the inverse tables ahead of time and load them in for lookup later? As your code shows, the ROM Sine table's 2k words, the inverse would be a big chunk too... just kind of winging it here - the math stuff fascinates me.
cheers
- Howard
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
The problem right now is that I have to search through the table.
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Nyamekye,
-Phil
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Here is some pseudo code to test the idea using Quickbasic that uses a binary search successive approximation method.
Initially a Sine Table is built to resemble the Sine table within the Propeller.· From there the Angle is derived after 11 itterations.
I did end up writing Propeller code for this but can't seem to find it at the moment.
ARCSIN = 1 ' Range from -1 to 1 Pi# = ATN(1) * 4 DIM SineTable(2048) Index = 0 FOR Deg = 0 TO 90 STEP 90 / 2048 'Build Sine Table Rad = (Deg / 180) * Pi# SineTable(Index) = SIN(Rad) Index = Index + 1 NEXT Deg Index = Index - 1 RefHigh = Index 'Find ArcSin Value RefLow = 0 Pivot = RefHigh \ 2 Sign = SGN(ARCSIN) ARCSIN = ABS(ARCSIN) FOR Itterations = 1 TO 11 IF SineTable(Pivot) > ARCSIN THEN RefHigh = Pivot ELSE RefLow = Pivot END IF Pivot = (RefHigh + RefLow) \ 2 NEXT Itterations ARCSINDeg = Sign * ((Pivot * 90) / 2048) ARCSINRad = Sign * ((Pivot * (Pi# / 2)) / 2048) CLS PRINT ARCSINDeg 'Display ArcSin Value PRINT ARCSINRad▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Beau Schwabe
IC Layout Engineer
Parallax, Inc.
Post Edited (Beau Schwabe (Parallax)) : 7/26/2009 4:55:48 AM GMT
Seems like there should be a way to do that but I can't seems to think of anything.
Guess I'll need to do bsearch.
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Nyamekye,
I don't have the arc-x tables in front of me, but I bet if you stare at them long enough, you'll see patterns that might allow you to reduce/normalize them (e.g. what the calculus tricks would do for you.) But we have to realize that those tables are *already* reduced. Obviously, the more data you drop out, the more of a quantization-like error you introduce.
- H
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
I used the technique of starting with ATN.· For ATN, reduce the angle to < PI/12.· I just found again the concept I used on this page (towards the bottom, the FORTRAN code):
http://www.dsprelated.com/showmessage/21872/1.php
For arccosine, I did this:
acos=pi/2-atan(x/sqrt(1-x*x))
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
My Prop Info&Apps: ·http://www.rayslogic.com/propeller/propeller.htm
The most useful inverse function you could program would be a four-quadrant arctan: atan(x, y). As Rayman points out, the other functions can follow from arctan. Frankly, it's probably not worth your while to optimise arcsin and arccos beyond that, since they're very rarely used for anything anyway.
-Phil
Actually that's a better idea because I would need to use the same look up table for the arccos function since I can't just add 90 like in the cos function.
Okay, thanks.
Supporting only a few is not optional however as I'm making library drivers.
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Nyamekye,
As Phil and Rayman pointed out, the ATAN2(X, Y) routine can be a good starting point for your integer math trig routines. It can be done with the CORDIC algorithm simply and efficiently. You can find PASM code for that CORDIC in the 'H48C Tri-Axis Accelerometer' or in the 'HM55B Compass Module Asm' drivers of Beau Schwabe. This code does cartesian-to-polar conversion.
CORDIC with suitably chosen inputs·calculates a whole range of scientific functions including; SIN, COS, TAN, ATAN, ANGLE, ASIN, ACOS, SINH, COSH, TANH, ATANH, LOG, EXP, SQRT and even multiply and divide, although CORDIC itself does not use multiplication. It uses only addition, subtraction, bitshift and table lookup. So, in your future integer math library you·may have repeated and efficient uses of those fast CORDIC routines. They can be organized in linear, circular and hyperbolic kinds, and they can be built upon each other, resulting in a compact code.
Istvan··
I'm doing more or less the same for an integer only navigation library using planar projection.
I'm planning to post my planar projection an dinteger only math ideas on wikispaces, as soon as I'll find time. Coupled I have a modified atan2 (using an object from Beau Schwabe) obtaining both angle and distance.
At the moment it is rough on the edges, and I don't know if it will ever be worth obex, but I can post it if you are interested.
Massimo