How to compute ATAN2 with integers? (and, of course, quickly)
cessnapilot
Posts: 182
Hi, All
I only show the main points in my vague translation of a Babylonian tablet into SPIN like pseudocode. They were not very sophisticated, so they called ATAN2 as ANGLE, X and Y meant steps forward/backward and left/rigth, and they measured the ANGLE in turn from right to the left. So, 10 steps forward, the turnig to the left, and another 10 steps gives 45 deg ANGLE. To calculate ANGLE from the X,Y coordinates (within +-1 degree) they used only integers.
(With a loupe I figured out some scratched lines, and I got totally confused, but the final version is maybe)
IF (Y==0)
· IF (X=>0)
··· RETURN 0
· RETURN 180·
··
IF (X==0)
· IF (Y>0)
··· RETURN 90
· RETURN -90
IF (ABS(Y)=<ABS(X))
· ANGLE := (3667 * X * Y) / (64 * X * X + 17 * Y * Y)
ELSE
· ANGLE := 90 - (3667 * X * Y) / (64 * Y * Y + 17 * X * X)
IF (X<0)
· IF (Y<0)
··· ANGLE := 180 + ANGLE
· ANGLE := 180 - ANGLE
ELSE
· IF (Y<0)
··· ANGLE := 360 -ANGLE
RETURN ANGLE
I might mixed up the quadrants, as there were scratches on the tablet. That must be Babylonian, as the ANGLE is in degree.· Although they know the 355/113 number. This or similar code for SINE, COSINE, or of the inverses, can save you 2 COGs sometimes.
This code·can be very fast in PASM and easy to program in SPIN. Note that 64x·and 17x multiplications. They can be done altogether with 2 shifts and·1 addition in assembly.
Cheers,
Istvan
I only show the main points in my vague translation of a Babylonian tablet into SPIN like pseudocode. They were not very sophisticated, so they called ATAN2 as ANGLE, X and Y meant steps forward/backward and left/rigth, and they measured the ANGLE in turn from right to the left. So, 10 steps forward, the turnig to the left, and another 10 steps gives 45 deg ANGLE. To calculate ANGLE from the X,Y coordinates (within +-1 degree) they used only integers.
(With a loupe I figured out some scratched lines, and I got totally confused, but the final version is maybe)
IF (Y==0)
· IF (X=>0)
··· RETURN 0
· RETURN 180·
··
IF (X==0)
· IF (Y>0)
··· RETURN 90
· RETURN -90
IF (ABS(Y)=<ABS(X))
· ANGLE := (3667 * X * Y) / (64 * X * X + 17 * Y * Y)
ELSE
· ANGLE := 90 - (3667 * X * Y) / (64 * Y * Y + 17 * X * X)
IF (X<0)
· IF (Y<0)
··· ANGLE := 180 + ANGLE
· ANGLE := 180 - ANGLE
ELSE
· IF (Y<0)
··· ANGLE := 360 -ANGLE
RETURN ANGLE
I might mixed up the quadrants, as there were scratches on the tablet. That must be Babylonian, as the ANGLE is in degree.· Although they know the 355/113 number. This or similar code for SINE, COSINE, or of the inverses, can save you 2 COGs sometimes.
This code·can be very fast in PASM and easy to program in SPIN. Note that 64x·and 17x multiplications. They can be done altogether with 2 shifts and·1 addition in assembly.
Cheers,
Istvan
Comments
Check http://forums.parallax.com/showthread.php?p=609053 for a couple of ideas.
I'm using a modified object from Beau Schwabe (I added the cordic gain rectification) and it is low footprint, leaving a lot of space in the cog for other operations. It would be nice to compare the performance of the two approaches.
Sine and cosine are already built in, so I guess if the resolution is enough it should be faster and smaller..
Thanks for directing me to that great thread. Two things came upon my mind while reading.· The 'undiscovered algorithm' for polar coordinate DSP is probably hiding in the deepness of the
e^(i*pi)+1=0,
equation, where all +* is done by the CORDIC algorithm. Maybe the four quarks (components) of quaternions are the 'elementary' particles of the atoms (rotations) of CORDIC? The second was, that nothing can prevent us to 'join' the CORDIC FFT with the WINOGRAD method. Hopefully, the advantages of both not just sum up, but multiply. Perhaps we can try, at least the second idea.
In the soft regime of computer methods, there are a lot of·new·ones claiming novelty and efficiency, and almost all of them are seeking for·applications.·CORDIC, however, is a nice example for me, where an application needed an efficient algorithm, and·found it. I like this sequence of events·much more·than the previous one.
That 'integer ANGLE' came from a back of the envelope Floating-point minimax approximation of ATAN. It's precision is about 0.005 radians, and that is remarkable with one parameter! It has only 2 multiplications and 2 divisions. (Something like Q=X/Y; ANGLE=Q/(1.0+0.24*Q*Q)). I only adapted it to integer arguments. With CORDIC, you would need about 12-16 iteration to reach that precision, I guess. Before any real comparison of running algorithms we should define the rules. The microcontroller and it's instruction set is given in our case, and we have to define the number format and the required precision, as well.
For· example, in the FPU or in the FloatMath lib, if you relax somewhat on the precision, that simple formula can be faster, than the built in ATAN2 (and CORDIC will beat both with Q numbers). But, if you are keen on precision, the simple formula needs a following Newton-Raphson to double the valid digits and library ATAN2 is competitive. At high precision, CORDIC suffers from rounding errors, huge tables and from many iterations (remember rounding errors, again).
As for the instruction set, if the microcontroller has a hardware multiplier, or at least CLZ (Count Leading Zeroes) and a CTZ (Count Trailing Zeroes) instructions, multiplication/division became quick, an so looses CORDIC it's main advantage. Anyway, I aggree with you and for the Propeller, there will (should) be a flexible CORDIC(IQF) library for the benefit of all of us, 'Propeller users'.
Cheers,
Istvan··
As usual a best all around solution is only a dream, and any choice is only a question of tradeoff.
I'm working on a small object to be used for integer only GPS navigation, using cartesian planar projection (there is a partial page on wikispace). Thanks to the 32 bits of the Propeller the system is fast and rather accurate.
So any integer approach is really interesting. After all the Cordic is an extremely old algorithm (it was used on the HP35 calculators, which I found on the internet were called the "slide rule killers").. In fact the first articles on Cordic are from HP people. Anyway with a different hardware, or a different level of precision required, a different approach could prove to be a better solution..
Massimo
http://perso.ens-lyon.fr/jean-michel.muller/BKM94.pdf
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Composite NTSC sprite driver: Forum
NTSC & PAL driver templates: ObEx Forum
OnePinTVText driver: ObEx Forum
If I were a software salesman, I'd offer you a best all-round solution,·trust me. However, I would only like to collect working and proven algorithms that will be actually written/optimized for the Propeller and each of them will be tested by the clock and precision.·
You wrote that you work on integer only GPS navigation object, that is small. I don't think so, that is great,·and I bet that you found interesting that ANGLE. If you have a better, quicker or more precise one, let us know, say within a month. If not, enjoy it. Other fast and precise integer trigonometric will be posted in the FFFC series. Concerning that, what method do you use to represent the outcome of sine or cosine with integers? Or for the angles? In other words, what scale factors are practical for GPS precision?
Cheers,
Istvan
Thats the mathematical positive direction!
Quadrants are numbered CCW.
Nick
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Never use force, just go for a bigger hammer!
The DIY Digital-Readout for mills, lathes etc.:
YADRO
That's because they used sundials those days. Well, that is only a convention. In math sometimes the angle is measured from the horizontal x-axis, upwards. In navigation, for example, its is just the opposite, as you correctly noticed. I tried to refer to this ambiquity writing
>I might mixed up the quadrants, as there were scratches on the tablet.
Insert one scratch more for a minus, and it is ok.
I have read a book about UFOs and ETs, where the author got to the conclusion that the people in Babylon knew wireless telegraphy (from the ETs), because not a single piece of wire was found in the ruins. If that is so, than they could even use time-machines to read current mathematical papers. (Or, that is not so.)
Cheers,
Istvan
I agree. One of the wisdoms of good programming practice says : "Know what you are doing.". To be a specialist in one field is nice and important, but it has its dangers. When you know more and more about less and less, finally you will know everything about nothing. Fortunately, good algorithms or solutions, like CORDIC, usually speak for themselves in layman's terms, too. If you are doomed or determined to program them, you will develop sharp senses to pick them out from a crowd of stupid patents and from the vast amount of technical literature written in greek. Google is your friend, as somebody wrote in this forum, but i was once told, it may happen that something is wrong on the internet, too.
(And as a practicing PPL pilot, I keep avoiding to be deep in any field, at least not too often.)
Cheers,
Istvan·
I'm trying to write a small page on wikispaces here: http://propeller.wikispaces.com/integer_navigation
At the moment there are only the basic concepts, but I'll get to the details..
Sine and cosine are scaled to a 32 bit long (-1 to +1 range). The sine cosine table in the propeller are already scaled this way. The multiply high is almost perfect, with the exception of the range (-0.5 +0.5). In the object I'm working on there is a modified multiply high that handles -1 +1 range, but you could use a spin multiply high and then shift left or multiply times 2.
The trig functions are using a 13bit resolution for the angles, less than 0.05 degrees.
For this system you only need sine cosine and atan2. I'm using fraction of prime to have a distance as a fraction of nautical mile. Using the 32 bits of the propeller I have a resolution of less than an inch, much better with respect to the GPS precision, so I have room for loosing something during the calculations.
Massimo